R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("DAAG-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('DAAG') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "ACF1" > > ### * ACF1 > > flush(stderr()); flush(stdout()) > > ### Name: ACF1 > ### Title: Aberrant Crypt Foci in Rat Colons > ### Aliases: ACF1 > ### Keywords: datasets > > ### ** Examples > > sapply(split(ACF1$count,ACF1$endtime),var) 6 12 18 2.333333 4.410714 2.571429 > plot(count ~ endtime, data=ACF1, pch=16) > pause() Pause. Press to continue...print("Poisson Regression - Example 8.3") > ACF.glm0 <- glm(formula = count ~ endtime, family = poisson, data = ACF1) > summary(ACF.glm0) Call: glm(formula = count ~ endtime, family = poisson, data = ACF1) Deviance Residuals: Min 1Q Median 3Q Max -2.46204 -0.47851 -0.07943 0.38159 2.26332 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.32152 0.40046 -0.803 0.422 endtime 0.11920 0.02642 4.511 6.44e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 51.105 on 21 degrees of freedom Residual deviance: 28.369 on 20 degrees of freedom AIC: 92.21 Number of Fisher Scoring iterations: 5 > > # Is there a quadratic effect? > pause() Pause. Press to continue... > ACF.glm <- glm(formula = count ~ endtime + I(endtime^2), + family = poisson, data = ACF1) > summary(ACF.glm) Call: glm(formula = count ~ endtime + I(endtime^2), family = poisson, data = ACF1) Deviance Residuals: Min 1Q Median 3Q Max -2.0616 -0.7834 -0.2808 0.4510 2.1693 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.722364 1.092494 1.577 0.115 endtime -0.262356 0.199685 -1.314 0.189 I(endtime^2) 0.015137 0.007954 1.903 0.057 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 51.105 on 21 degrees of freedom Residual deviance: 24.515 on 19 degrees of freedom AIC: 90.354 Number of Fisher Scoring iterations: 5 > > # But is the data really Poisson? If not, try quasipoisson: > pause() Pause. Press to continue... > ACF.glm <- glm(formula = count ~ endtime + I(endtime^2), + family = quasipoisson, data = ACF1) > summary(ACF.glm) Call: glm(formula = count ~ endtime + I(endtime^2), family = quasipoisson, data = ACF1) Deviance Residuals: Min 1Q Median 3Q Max -2.0616 -0.7834 -0.2808 0.4510 2.1693 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.722364 1.223430 1.408 0.175 endtime -0.262356 0.223618 -1.173 0.255 I(endtime^2) 0.015137 0.008908 1.699 0.106 (Dispersion parameter for quasipoisson family taken to be 1.254065) Null deviance: 51.105 on 21 degrees of freedom Residual deviance: 24.515 on 19 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 > > > > cleanEx(); ..nameEx <- "Cars93.summary" > > ### * Cars93.summary > > flush(stderr()); flush(stdout()) > > ### Name: Cars93.summary > ### Title: A Summary of the Cars93 Data set > ### Aliases: Cars93.summary > ### Keywords: datasets > > ### ** Examples > > type <- Cars93.summary$abbrev > type <- Cars93.summary[,4] > type <- Cars93.summary[,"abbrev"] > type <- Cars93.summary[[4]] # Take the object that is stored > # in the fourth list element. > type [1] C L M Sm Sp V Levels: C L M Sm Sp V > pause() Pause. Press to continue... > attach(Cars93.summary) > # R can now access the columns of Cars93.summary directly > abbrev [1] C L M Sm Sp V Levels: C L M Sm Sp V > detach("Cars93.summary") > pause() Pause. Press to continue... > # To change the name of the \verb!abbrev! variable (the fourth column) > names(Cars93.summary)[4] <- "code" > pause() Pause. Press to continue... > # To change all of the names, try > names(Cars93.summary) <- c("minpass","maxpass","number","code") > > > > > cleanEx(); ..nameEx <- "Lottario" > > ### * Lottario > > flush(stderr()); flush(stdout()) > > ### Name: Lottario > ### Title: Ontario Lottery Data > ### Aliases: Lottario > ### Keywords: datasets > > ### ** Examples > > order(Lottario$Frequency)[33:39] # the 7 most frequently chosen numbers [1] 13 5 14 12 38 4 35 > > > > cleanEx(); ..nameEx <- "Manitoba.lakes" > > ### * Manitoba.lakes > > flush(stderr()); flush(stdout()) > > ### Name: Manitoba.lakes > ### Title: The Nine Largest Lakes in Manitoba > ### Aliases: Manitoba.lakes > ### Keywords: datasets > > ### ** Examples > > plot(Manitoba.lakes) > plot(Manitoba.lakes[-1,]) > > > > cleanEx(); ..nameEx <- "allbacks" > > ### * allbacks > > flush(stderr()); flush(stdout()) > > ### Name: allbacks > ### Title: Measurements on a Selection of Books > ### Aliases: allbacks > ### Keywords: datasets > > ### ** Examples > > print("Multiple Regression - Example 6.1") [1] "Multiple Regression - Example 6.1" > attach(allbacks) > volume.split <- split(volume, cover) > weight.split <- split(weight, cover) > plot(weight.split$hb ~ volume.split$hb, pch=16, xlim=range(volume), ylim=range(weight), + ylab="Weight (g)", xlab="Volume (cc)") > points(weight.split$pb ~ volume.split$pb, pch=16, col=2) > pause() Pause. Press to continue... > allbacks.lm <- lm(weight ~ volume+area) > summary(allbacks.lm) Call: lm(formula = weight ~ volume + area) Residuals: Min 1Q Median 3Q Max -104.06 -30.02 -15.46 16.76 212.30 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 22.41342 58.40247 0.384 0.707858 volume 0.70821 0.06107 11.597 7.07e-08 *** area 0.46843 0.10195 4.595 0.000616 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 77.66 on 12 degrees of freedom Multiple R-Squared: 0.9285, Adjusted R-squared: 0.9166 F-statistic: 77.89 on 2 and 12 DF, p-value: 1.339e-07 > detach(allbacks) > pause() Pause. Press to continue... > anova(allbacks.lm) Analysis of Variance Table Response: weight Df Sum Sq Mean Sq F value Pr(>F) volume 1 812132 812132 134.659 7.02e-08 *** area 1 127328 127328 21.112 0.0006165 *** Residuals 12 72373 6031 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > pause() Pause. Press to continue... > model.matrix(allbacks.lm) (Intercept) volume area 1 1 885 382 2 1 1016 468 3 1 1125 387 4 1 239 371 5 1 701 371 6 1 641 367 7 1 1228 396 8 1 412 0 9 1 953 0 10 1 929 0 11 1 1492 0 12 1 419 0 13 1 1010 0 14 1 595 0 15 1 1034 0 attr(,"assign") [1] 0 1 2 > pause() Pause. Press to continue... > print("Example 6.1.1") [1] "Example 6.1.1" > allbacks.lm0 <- lm(weight ~ -1+volume+area, data=allbacks); > summary(allbacks.lm0) Call: lm(formula = weight ~ -1 + volume + area, data = allbacks) Residuals: Min 1Q Median 3Q Max -112.53 -28.73 -10.52 24.62 213.80 Coefficients: Estimate Std. Error t value Pr(>|t|) volume 0.72891 0.02767 26.344 1.15e-12 *** area 0.48087 0.09344 5.146 0.000188 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 75.07 on 13 degrees of freedom Multiple R-Squared: 0.9914, Adjusted R-squared: 0.9901 F-statistic: 747.9 on 2 and 13 DF, p-value: 3.799e-14 > pause() Pause. Press to continue... > print("Example 6.1.2") [1] "Example 6.1.2" > oldpar <- par(mfrow=c(2,2)) > plot(allbacks.lm0) > par(oldpar) > allbacks.lm13 <- lm(weight ~ -1+volume+area, data=allbacks[-13,]) > summary(allbacks.lm13) Call: lm(formula = weight ~ -1 + volume + area, data = allbacks[-13, ]) Residuals: Min 1Q Median 3Q Max -61.721 -25.291 3.429 31.244 58.856 Coefficients: Estimate Std. Error t value Pr(>|t|) volume 0.69485 0.01629 42.65 1.79e-14 *** area 0.55390 0.05269 10.51 2.08e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 41.02 on 12 degrees of freedom Multiple R-Squared: 0.9973, Adjusted R-squared: 0.9969 F-statistic: 2252 on 2 and 12 DF, p-value: 3.521e-16 > pause() Pause. Press to continue... > print("Example 6.1.3") [1] "Example 6.1.3" > round(coef(allbacks.lm0),2) # Baseline for changes volume area 0.73 0.48 > round(lm.influence(allbacks.lm0)$coef,2) volume area 1 0.00 -0.01 2 0.00 -0.01 3 0.00 0.01 4 0.00 0.00 5 0.00 0.03 6 0.00 -0.02 7 0.00 0.00 8 0.00 0.01 9 0.00 0.00 10 0.00 0.01 11 -0.03 0.07 12 0.00 -0.01 13 0.03 -0.07 14 0.00 0.00 15 0.00 0.01 > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "anesthetic" > > ### * anesthetic > > flush(stderr()); flush(stdout()) > > ### Name: anesthetic > ### Title: Anesthetic Effectiveness > ### Aliases: anesthetic > ### Keywords: datasets > > ### ** Examples > > print("Logistic Regression - Example 8.1.4") [1] "Logistic Regression - Example 8.1.4" > > z <- table(anesthetic$nomove, anesthetic$conc) > tot <- apply(z, 2, sum) # totals at each concentration > prop <- z[2, ]/(tot) # proportions at each concentration > oprop <- sum(z[2, ])/sum(tot) # expected proportion moving if concentration had no effect > conc <- as.numeric(dimnames(z)[[2]]) > plot(conc, prop, xlab = "Concentration", ylab = "Proportion", xlim = c(.5,2.5), + ylim = c(0, 1), pch = 16) > chw <- par()$cxy[1] > text(conc - 0.75 * chw, prop, paste(tot), adj = 1) > abline(h = oprop, lty = 2) > > pause() Pause. Press to continue... > anes.logit <- glm(nomove ~ conc, family = binomial(link = logit), + data = anesthetic) > anova(anes.logit) Analysis of Deviance Table Model: binomial, link: logit Response: nomove Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 29 41.455 conc 1 13.701 28 27.754 > summary(anes.logit) Call: glm(formula = nomove ~ conc, family = binomial(link = logit), data = anesthetic) Deviance Residuals: Min 1Q Median 3Q Max -1.76666 -0.74407 0.03413 0.68666 2.06900 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -6.469 2.418 -2.675 0.00748 ** conc 5.567 2.044 2.724 0.00645 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 41.455 on 29 degrees of freedom Residual deviance: 27.754 on 28 degrees of freedom AIC: 31.754 Number of Fisher Scoring iterations: 5 > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "austpop" > > ### * austpop > > flush(stderr()); flush(stdout()) > > ### Name: austpop > ### Title: Population figures for Australian States and Territories > ### Aliases: austpop > ### Keywords: datasets > > ### ** Examples > > print("Looping - Example 1.7") [1] "Looping - Example 1.7" > > growth.rates <- numeric(8) > for (j in seq(2,9)) { + growth.rates[j-1] <- (austpop[9, j]-austpop[1, j])/austpop[1, j] } > growth.rates <- data.frame(growth.rates) > row.names(growth.rates) <- names(austpop[c(-1,-10)]) > # Note the use of row.names() to name the rows of the data frame > growth.rates growth.rates NSW 2.295168 Vic 2.268275 Qld 3.979502 SA 2.363636 WA 4.875817 Tas 1.455959 NT 36.400000 ACT 102.333333 > > pause() Pause. Press to continue...print("Avoiding Loops - Example 1.7b") > > sapply(austpop[,-c(1,10)], function(x){(x[9]-x[1])/x[1]}) NSW Vic Qld SA WA Tas NT 2.295168 2.268275 3.979502 2.363636 4.875817 1.455959 36.400000 ACT 102.333333 > > pause() Pause. Press to continue...print("Plot - Example 1.8a") > attach(austpop) > plot(year, ACT, type="l") # Join the points ("l" = "line") > detach(austpop) > > pause() Pause. Press to continue...print("Exerice 1.12.9") > attach(austpop) > oldpar <- par(mfrow=c(2,4)) > for (i in 2:9){ + plot(austpop[,1], log(austpop[, i]), xlab="Year", + ylab=names(austpop)[i], pch=16, ylim=c(0,10))} > par(oldpar) > detach(austpop) > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "bestset.noise" > > ### * bestset.noise > > flush(stderr()); flush(stdout()) > > ### Name: bestset.noise > ### Title: Best Subset Selection Applied to Noise > ### Aliases: bestset.noise > ### Keywords: models > > ### ** Examples > > leaps.out <- try(require(leaps, silent=TRUE)) Error in require(leaps, silent = TRUE) : unused argument(s) (silent ...) > leaps.out.log <- is.logical(leaps.out) > if ((leaps.out.log==TRUE)&(leaps.out==TRUE)) + bestset.noise(20,6) # `best' 3-variable regression for 20 simulated observations > # on 7 unrelated variables (including the response) > > > > cleanEx(); ..nameEx <- "bomsoi" > > ### * bomsoi > > flush(stderr()); flush(stdout()) > > ### Name: bomsoi > ### Title: Southern Oscillation Index Data > ### Aliases: bomsoi > ### Keywords: datasets > > ### ** Examples > > require(ts) Loading required package: ts Warning: package 'ts' has been merged into 'stats' [1] TRUE > plot(ts(bomsoi[, 15:14], start=1900), + panel=function(y,...)panel.smooth(1900:2001, y,...)) > > # Check for skewness by comparing the normal probability plots for > # different a, e.g. > par(mfrow = c(2,3)) > for (a in c(50, 100, 150, 200, 250, 300)) + qqnorm(log(bomsoi[, "avrain"] - a)) > # a = 250 leads to a nearly linear plot > par(mfrow = c(1,1)) > > plot(bomsoi$SOI, log(bomsoi$avrain - 250), xlab = "SOI", + ylab = "log(avrain = 250)") > lines(lowess(bomsoi$SOI)$y, lowess(log(bomsoi$avrain - 250))$y, lwd=2) > # NB: separate lowess fits against time > lines(lowess(bomsoi$SOI, log(bomsoi$avrain - 250))) > > detsoi <- data.frame( + detSOI = bomsoi[, "SOI"] - lowess(bomsoi[, "SOI"])$y, + detrain = log(bomsoi$avrain - 250) - lowess(log(bomsoi$avrain - 250))$y) > row.names(detsoi) <- paste(1900:2001) > > par(mfrow = c(1,2)) > plot(log(avrain-250) ~ SOI, data = bomsoi, ylab = + "log(Average rainfall - 250)") > lines(lowess(bomsoi$SOI, log(bomsoi$avrain-250))) > plot(detrain ~ detSOI, data = detsoi, + xlab="Detrended SOI", ylab = "Detrended log(Rainfall-250)") > lines(lowess(detsoi$detrain ~ detsoi$detSOI)) > par(mfrow = c(1,1)) > > require(nlme) Loading required package: nlme [1] TRUE > soi.gls <- gls(detrain ~ detSOI, data = detsoi, correlation = + corARMA(q=12)) > summary(soi.gls) Generalized least squares fit by REML Model: detrain ~ detSOI Data: detsoi AIC BIC logLik 70.41511 109.4927 -20.20756 Correlation Structure: ARMA(0,12) Formula: ~1 Parameter estimate(s): Theta1 Theta2 Theta3 Theta4 Theta5 Theta6 -0.06299578 -0.09703185 0.01605314 -0.13600311 -0.10565381 -0.12503400 Theta7 Theta8 Theta9 Theta10 Theta11 Theta12 -0.02863833 -0.03085423 0.10950782 0.23565425 -0.28234185 -0.47541693 Coefficients: Value Std.Error t-value p-value (Intercept) -0.01806652 0.006143586 -2.940713 0.0041 detSOI 0.03379488 0.003289733 10.272836 0.0000 Correlation: (Intr) detSOI -0.059 Standardized residuals: Min Q1 Med Q3 Max -3.1196396 -0.5932111 0.1264759 0.8330351 2.4050761 Residual standard error: 0.3133226 Degrees of freedom: 102 total; 100 residual > > soi1ML.gls <- update(soi.gls, method = "ML") > soi0ML.gls <- update(soi.gls, detrain ~ 1, method = "ML") > soi2ML.gls <- update(soi.gls, detrain ~ detSOI + detSOI^2, method = "ML") > anova(soi2ML.gls, soi1ML.gls) Model df AIC BIC logLik soi2ML.gls 1 15 52.42128 91.79588 -11.21064 soi1ML.gls 2 15 52.42128 91.79588 -11.21064 > > # compare with MA(11) and MA(13) > soi11.gls <- update(soi.gls, correlation=corARMA(q=11)) > soi13.gls <- update(soi.gls, correlation=corARMA(q=13)) > anova(soi11.gls, soi.gls, soi13.gls) Model df AIC BIC logLik Test L.Ratio p-value soi11.gls 1 14 81.63582 118.1082 -26.81791 soi.gls 2 15 70.41511 109.4927 -20.20756 1 vs 2 13.220708 0.0003 soi13.gls 3 16 71.55681 113.2395 -19.77841 2 vs 3 0.858299 0.3542 > > # compare with the white noise model > soi0.gls <- gls(detrain ~ detSOI, data=detsoi) > anova(soi0.gls, soi.gls) Model df AIC BIC logLik Test L.Ratio p-value soi0.gls 1 3 73.88653 81.70204 -33.94327 soi.gls 2 15 70.41511 109.49266 -20.20756 1 vs 2 27.47142 0.0066 > > # a Portmanteau test of whiteness for the white noise model residuals > > Box.test(resid(soi0.gls), lag=20, type="Ljung-Box") Box-Ljung test data: resid(soi0.gls) X-squared = 33.058, df = 20, p-value = 0.03325 > # check residual properties > acf(resid(soi.gls)) # (Correlated) residuals > acf(resid(soi.gls, type="normalized")) # Innovation estimates, uncorrelated > qqnorm(resid(soi.gls, type="normalized")) ## Examine normality > # Now extract the moving average parameters, and plot the > # the theoretical autocorrelation function that they imply. > beta <- summary(soi.gls$modelStruct)$corStruct > plot(ARMAacf(ma=beta,lag.max=20), type="h") > # Next, plot several simulated autocorrelation functions > # We can plot autocorrelation functions as though they were time series! > > plot.ts(ts(cbind( + "Series 1" = acf(arima.sim(list(ma=beta), n=83), plot=FALSE, + lag.max = 20)$acf, + "Series 2" = acf(arima.sim(list(ma=beta), n=83), plot=FALSE, + lag.max = 20)$acf, + "Series 3" = acf(arima.sim(list(ma=beta), n=83), plot=FALSE, + lag.max = 20)$acf), start=0), type="h", main = "", xlab = "Lag") > # Show confidence bounds for the MA parameters > intervals(soi.gls) Approximate 95% confidence intervals Coefficients: lower est. upper (Intercept) -0.03025522 -0.01806652 -0.005877821 detSOI 0.02726815 0.03379488 0.040321620 attr(,"label") [1] "Coefficients:" Correlation structure: lower est. upper Theta1 5.9031035 -0.06299578 6.4457770 Theta2 15.1634507 -0.09703185 19.2914786 Theta3 22.0080444 0.01605314 35.9501855 Theta4 19.4264295 -0.13600311 47.4265066 Theta5 9.4378374 -0.10565381 47.4379545 Theta6 -2.0691208 -0.12503400 36.2498085 Theta7 -12.9078515 -0.02863833 19.0675821 Theta8 -20.5364817 -0.03085423 3.7942682 Theta9 -20.1241744 0.10950782 -3.5666035 Theta10 -12.3207525 0.23565425 -3.6446151 Theta11 -4.3169596 -0.28234185 -1.4484117 Theta12 -0.6635251 -0.47541693 -0.2307739 attr(,"label") [1] "Correlation structure:" Residual standard error: lower est. upper 0.2665452 0.3133226 0.3683092 > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "carprice" > > ### * carprice > > flush(stderr()); flush(stdout()) > > ### Name: carprice > ### Title: US Car Price Data > ### Aliases: carprice > ### Keywords: datasets > > ### ** Examples > > print("Multicollinearity - Example 6.8") [1] "Multicollinearity - Example 6.8" > pairs(carprice[,-c(1,8,9)]) > > carprice1.lm <- lm(gpm100 ~ Type+Min.Price+Price+Max.Price+Range.Price, + data=carprice) > round(summary(carprice1.lm)$coef,3) Estimate Std. Error t value Pr(>|t|) (Intercept) 3.287 0.153 21.467 0.000 TypeLarge 0.324 0.173 1.867 0.069 TypeMidsize 0.185 0.166 1.114 0.272 TypeSmall -0.389 0.168 -2.317 0.026 TypeSporty 0.205 0.175 1.176 0.247 TypeVan 1.349 0.193 6.997 0.000 Min.Price 0.700 0.989 0.707 0.484 Price -1.377 1.983 -0.695 0.491 Max.Price 0.711 0.994 0.715 0.479 > pause() Pause. Press to continue... > alias(carprice1.lm) Model : gpm100 ~ Type + Min.Price + Price + Max.Price + Range.Price Complete : (Intercept) TypeLarge TypeMidsize TypeSmall TypeSporty TypeVan Range.Price Min.Price Price Max.Price Range.Price -1 1 > pause() Pause. Press to continue... > carprice2.lm <- lm(gpm100 ~ Type+Min.Price+Price+Max.Price+RoughRange, data=carprice) > round(summary(carprice2.lm)$coef, 2) Estimate Std. Error t value Pr(>|t|) (Intercept) 3.29 0.16 21.19 0.00 TypeLarge 0.33 0.18 1.84 0.07 TypeMidsize 0.18 0.17 1.09 0.28 TypeSmall -0.39 0.17 -2.29 0.03 TypeSporty 0.21 0.18 1.16 0.25 TypeVan 1.35 0.20 6.90 0.00 Min.Price 0.26 4.28 0.06 0.95 Price -1.35 2.02 -0.67 0.51 Max.Price 1.13 4.04 0.28 0.78 RoughRange -0.43 4.04 -0.11 0.92 > pause() Pause. Press to continue... > carprice.lm <- lm(gpm100 ~ Type + Price, data = carprice) > round(summary(carprice.lm)$coef,4) Estimate Std. Error t value Pr(>|t|) (Intercept) 3.3047 0.1478 22.3566 0.0000 TypeLarge 0.3288 0.1681 1.9553 0.0574 TypeMidsize 0.1957 0.1622 1.2066 0.2345 TypeSmall -0.3651 0.1620 -2.2541 0.0296 TypeSporty 0.2426 0.1629 1.4897 0.1440 TypeVan 1.3920 0.1805 7.7121 0.0000 Price 0.0330 0.0074 4.4820 0.0001 > pause() Pause. Press to continue... > summary(carprice1.lm)$sigma # residual standard error when fitting all 3 price variables [1] 0.3059430 > pause() Pause. Press to continue... > summary(carprice.lm)$sigma # residual standard error when only price is used [1] 0.3005687 > pause() Pause. Press to continue... > vif(lm(gpm100 ~ Price, data=carprice)) # Baseline Price Price 1 > pause() Pause. Press to continue... > vif(carprice1.lm) # includes Min.Price, Price & Max.Price TypeLarge TypeMidsize TypeSmall TypeSporty TypeVan Min.Price 2.7216e+00 2.3316e+00 1.8043e+00 2.1740e+00 1.7779e+00 2.9919e+04 Price Max.Price 1.2064e+05 3.3400e+04 > pause() Pause. Press to continue... > vif(carprice2.lm) # includes Min.Price, Price, Max.Price & RoughRange TypeLarge TypeMidsize TypeSmall TypeSporty TypeVan Min.Price 2.7800e+00 2.3466e+00 1.8237e+00 2.1757e+00 1.7845e+00 5.4494e+05 Price Max.Price RoughRange 1.2239e+05 5.3809e+05 9.6219e+04 > pause() Pause. Press to continue... > vif(carprice.lm) # Price alone TypeLarge TypeMidsize TypeSmall TypeSporty TypeVan Price 2.6534 2.3042 1.7362 1.9576 1.6152 1.7277 > > > > cleanEx(); ..nameEx <- "cfseal" > > ### * cfseal > > flush(stderr()); flush(stdout()) > > ### Name: cfseal > ### Title: Cape Fur Seal Data > ### Aliases: cfseal > ### Keywords: datasets > > ### ** Examples > > print("Allometric Growth - Example 5.7") [1] "Allometric Growth - Example 5.7" > > cfseal.lm <- lm(log(heart) ~ log(weight), data=cfseal); summary(cfseal.lm) Call: lm(formula = log(heart) ~ log(weight), data = cfseal) Residuals: Min 1Q Median 3Q Max -3.149e-01 -9.182e-02 2.475e-05 1.168e-01 3.205e-01 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.20434 0.21131 5.699 4.12e-06 *** log(weight) 1.12615 0.05467 20.597 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1796 on 28 degrees of freedom Multiple R-Squared: 0.9381, Adjusted R-squared: 0.9359 F-statistic: 424.3 on 1 and 28 DF, p-value: < 2.2e-16 > plot(log(heart) ~ log(weight), data = cfseal, pch=16, xlab = "Heart Weight (g, log scale)", + ylab = "Body weight (kg, log scale)", axes=FALSE) > heartaxis <- 100*(2^seq(0,3)) > bodyaxis <- c(20,40,60,100,180) > axis(1, at = log(bodyaxis), lab = bodyaxis) > axis(2, at = log(heartaxis), lab = heartaxis) > box() > abline(cfseal.lm) > > > > cleanEx(); ..nameEx <- "cities" > > ### * cities > > flush(stderr()); flush(stdout()) > > ### Name: cities > ### Title: Populations of Major Canadian Cities (1992-96) > ### Aliases: cities > ### Keywords: datasets > > ### ** Examples > > cities$have <- factor((cities$REGION=="ON")|(cities$REGION=="WEST")) > plot(POP1996~POP1992, data=cities, col=as.integer(cities$have)) > > > > cleanEx(); ..nameEx <- "component.residual" > > ### * component.residual > > flush(stderr()); flush(stdout()) > > ### Name: component.residual > ### Title: Component + Residual Plot > ### Aliases: component.residual > ### Keywords: models > > ### ** Examples > > mice12.lm <- lm(brainwt ~ bodywt + lsize, data=litters) > oldpar <- par(mfrow = c(1,2)) > component.residual(mice12.lm, 1, xlab = "Body weight", ylab= "t(Body weight) + e") > component.residual(mice12.lm, 2, xlab = "Litter size", ylab= "t(Litter size) + e") > par(oldpar) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "cuckoos" > > ### * cuckoos > > flush(stderr()); flush(stdout()) > > ### Name: cuckoos > ### Title: Cuckoo Eggs Data > ### Aliases: cuckoos > ### Keywords: datasets > > ### ** Examples > > print("Strip and Boxplots - Example 2.1.2") [1] "Strip and Boxplots - Example 2.1.2" > > attach(cuckoos) > oldpar <- par(las = 2) # labels at right angle to axis. > stripchart(length ~ species) > boxplot(split(cuckoos$length, cuckoos$species), + xlab="Length of egg", horizontal=TRUE) > detach(cuckoos) > par(oldpar) > pause() Pause. Press to continue... > print("Summaries - Example 2.2.2") [1] "Summaries - Example 2.2.2" > sapply(split(cuckoos$length, cuckoos$species), sd) hedge.sparrow meadow.pipit pied.wagtail robin tree.pipit 1.0494373 0.9195849 1.0722917 0.6821229 0.8800974 wren 0.7542262 > pause() Pause. Press to continue... > print("Example 4.1.4") [1] "Example 4.1.4" > wren <- split(cuckoos$length, cuckoos$species)$wren > median(wren) [1] 21 > n <- length(wren) > sqrt(pi/2)*sd(wren)/sqrt(n) # this s.e. computation assumes normality [1] 0.2440709 > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "cv.binary" > > ### * cv.binary > > flush(stderr()); flush(stdout()) > > ### Name: cv.binary > ### Title: Cross-Validation for Regression with a Binary Response > ### Aliases: cv.binary > ### Keywords: models > > ### ** Examples > > frogs.glm <- glm(pres.abs ~ log(distance) + log(NoOfPools), + family=binomial,data=frogs) > cv.binary(frogs.glm) Fold: 3 4 6 10 9 7 1 2 8 5 Internal estimate of accuracy = 0 Cross-validation estimate of accuracy = 0 > > > > cleanEx(); ..nameEx <- "cv.lm" > > ### * cv.lm > > flush(stderr()); flush(stdout()) > > ### Name: cv.lm > ### Title: Cross-Validation for Linear Regression > ### Aliases: cv.lm > ### Keywords: models > > ### ** Examples > > cv.lm() Analysis of Variance Table Response: yv Df Sum Sq Mean Sq F value Pr(>F) xv 1 18566 18566 8 0.014 * Residuals 13 30179 2321 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 fold 1 Observations in test set: 3 12 13 14 15 X1 X2 X3 X4 X5 area 802.0 696 771.0 1006.0 1191 Predicted 204.0 188 199.3 234.7 262 sale.price 215.0 255 260.0 293.0 375 Residual 11.0 67 60.7 58.3 113 Sum of squares = 24351 Mean square = 4870 n = 5 fold 2 Observations in test set: 2 5 6 9 10 X1 X2 X3 X4 X5 area 905 716 963.0 1018.00 887.00 Predicted 255 224 264.4 273.38 252.06 sale.price 215 113 185.0 276.00 260.00 Residual -40 -112 -79.4 2.62 7.94 Sum of squares = 20416 Mean square = 4083 n = 5 fold 3 Observations in test set: 1 4 7 8 11 X1 X2 X3 X4 X5 area 694.0 1366 821.00 714.0 790.00 Predicted 183.2 388 221.94 189.3 212.49 sale.price 192.0 274 212.00 220.0 221.50 Residual 8.8 -114 -9.94 30.7 9.01 Sum of squares = 14241 Mean square = 2848 n = 5 Overall ms 3934 > > > > cleanEx(); ..nameEx <- "dewpoint" > > ### * dewpoint > > flush(stderr()); flush(stdout()) > > ### Name: dewpoint > ### Title: Dewpoint Data > ### Aliases: dewpoint > ### Keywords: datasets > > ### ** Examples > > print("Additive Model - Example 7.5") [1] "Additive Model - Example 7.5" > require(splines) Loading required package: splines [1] TRUE > attach(dewpoint) > ds.lm <- lm(dewpt ~ bs(maxtemp,5) + bs(mintemp,5), data=dewpoint) > ds.fit <-predict(ds.lm, type="terms", se=TRUE) > oldpar <- par(mfrow=c(1,2)) > plot(maxtemp, ds.fit$fit[,1], xlab="Maximum temperature", + ylab="Change from dewpoint mean",type="n") > lines(maxtemp,ds.fit$fit[,1]) > lines(maxtemp,ds.fit$fit[,1]-2*ds.fit$se[,1],lty=2) > lines(maxtemp,ds.fit$fit[,1]+2*ds.fit$se[,1],lty=2) > plot(mintemp,ds.fit$fit[,2],xlab="Minimum temperature", + ylab="Change from dewpoint mean",type="n") > ord<-order(mintemp) > lines(mintemp[ord],ds.fit$fit[ord,2]) > lines(mintemp[ord],ds.fit$fit[ord,2]-2*ds.fit$se[ord,2],lty=2) > lines(mintemp[ord],ds.fit$fit[ord,2]+2*ds.fit$se[ord,2],lty=2) > detach(dewpoint) > par(oldpar) > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "elastic1" > > ### * elastic1 > > flush(stderr()); flush(stdout()) > > ### Name: elastic1 > ### Title: Elastic Band Data Replicated > ### Aliases: elastic1 > ### Keywords: datasets > > ### ** Examples > > plot(elastic1) > > print("Inline Functions - Example 12.2.2") [1] "Inline Functions - Example 12.2.2" > sapply(elastic1, mean) stretch distance 48 196 > pause() Pause. Press to continue... > sapply(elastic1, function(x)mean(x)) stretch distance 48 196 > pause() Pause. Press to continue... > sapply(elastic1, function(x)sum(log(x))) stretch distance 27.1 36.9 > pause() Pause. Press to continue... > print("Data Output - Example 12.3.2") [1] "Data Output - Example 12.3.2" > write.table(elastic1, file="bandsframe.txt") > > > > > cleanEx(); ..nameEx <- "elastic2" > > ### * elastic2 > > flush(stderr()); flush(stdout()) > > ### Name: elastic2 > ### Title: Elastic Band Data Replicated Again > ### Aliases: elastic2 > ### Keywords: datasets > > ### ** Examples > > plot(elastic2) > pause() Pause. Press to continue... > print("Chapter 5 Exercise") [1] "Chapter 5 Exercise" > > yrange <- range(c(elastic1$distance, elastic2$distance)) > xrange <- range(c(elastic1$stretch, elastic2$stretch)) > plot(distance ~ stretch, data = elastic1, pch = 16, ylim = yrange, xlim = + xrange) > points(distance ~ stretch, data = elastic2, pch = 15, col = 2) > legend(xrange[1], yrange[2], legend = c("Data set 1", "Data set 2"), pch = + c(16, 15), col = c(1, 2)) > > elastic1.lm <- lm(distance ~ stretch, data = elastic1) > elastic2.lm <- lm(distance ~ stretch, data = elastic2) > abline(elastic1.lm) > abline(elastic2.lm, col = 2) > summary(elastic1.lm) Call: lm(formula = distance ~ stretch, data = elastic1) Residuals: 1 2 3 4 5 6 7 -0.143 -18.714 -7.286 -1.429 8.000 -6.857 26.429 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -119.14 70.94 -1.68 0.1539 stretch 6.57 1.47 4.46 0.0066 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 15.6 on 5 degrees of freedom Multiple R-Squared: 0.799, Adjusted R-squared: 0.759 F-statistic: 19.9 on 1 and 5 DF, p-value: 0.00663 > summary(elastic2.lm) Call: lm(formula = distance ~ stretch, data = elastic2) Residuals: Min 1Q Median 3Q Max -10.083 -7.083 -0.583 5.167 20.167 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -100.917 15.610 -6.46 0.00035 *** stretch 5.950 0.315 18.90 2.9e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 10.4 on 7 degrees of freedom Multiple R-Squared: 0.981, Adjusted R-squared: 0.978 F-statistic: 357 on 1 and 7 DF, p-value: 2.89e-07 > pause() Pause. Press to continue... > predict(elastic1.lm, se.fit=TRUE) $fit 1 2 3 4 5 6 7 183 236 196 209 170 157 223 $se.fit [1] 6.59 10.62 5.89 6.59 8.33 10.62 8.33 $df [1] 5 $residual.scale [1] 15.6 > predict(elastic2.lm, se.fit=TRUE) $fit 1 2 3 4 5 6 7 8 9 77.6 196.6 137.1 166.8 256.1 226.3 107.3 226.3 285.8 $se.fit [1] 6.74 3.52 4.36 3.64 5.06 4.06 5.45 4.06 6.30 $df [1] 7 $residual.scale [1] 10.4 > > > > cleanEx(); ..nameEx <- "elasticband" > > ### * elasticband > > flush(stderr()); flush(stdout()) > > ### Name: elasticband > ### Title: Elastic Band Data > ### Aliases: elasticband > ### Keywords: datasets > > ### ** Examples > > print("Example 1.8.1") [1] "Example 1.8.1" > > attach(elasticband) # R now knows where to find stretch and distance > plot(stretch, distance) # Alternative: plot(distance ~ stretch) > detach(elasticband) > pause() Pause. Press to continue... > print("Output of Data Frames - Example 12.3.2") [1] "Output of Data Frames - Example 12.3.2" > > write(t(elasticband),file="bands.txt",ncol=2) > > sink("bands2.txt") > elasticband # NB: No output on screen > sink() > > print("Lists - Example 12.7") [1] "Lists - Example 12.7" > > elastic.lm <- lm(distance ~ stretch, data=elasticband) > names(elastic.lm) [1] "coefficients" "residuals" "effects" "rank" [5] "fitted.values" "assign" "qr" "df.residual" [9] "xlevels" "call" "terms" "model" > elastic.lm$coefficients (Intercept) stretch -63.57 4.55 > elastic.lm[["coefficients"]] (Intercept) stretch -63.57 4.55 > pause() Pause. Press to continue... > elastic.lm[[1]] (Intercept) stretch -63.57 4.55 > pause() Pause. Press to continue... > elastic.lm[1] $coefficients (Intercept) stretch -63.57 4.55 > pause() Pause. Press to continue... > options(digits=3) > elastic.lm$residuals 1 2 3 4 5 6 7 2.107 -0.321 18.000 1.893 -27.786 13.321 -7.214 > pause() Pause. Press to continue... > elastic.lm$call lm(formula = distance ~ stretch, data = elasticband) > pause() Pause. Press to continue... > mode(elastic.lm$call) [1] "call" > > > > > cleanEx(); ..nameEx <- "fossilfuel" > > ### * fossilfuel > > flush(stderr()); flush(stdout()) > > ### Name: fossilfuel > ### Title: Fossil Fuel Data > ### Aliases: fossilfuel > ### Keywords: datasets > > ### ** Examples > > plot(fossilfuel) > > > > cleanEx(); ..nameEx <- "frogs" > > ### * frogs > > flush(stderr()); flush(stdout()) > > ### Name: frogs > ### Title: Frogs Data > ### Aliases: frogs > ### Keywords: datasets > > ### ** Examples > > print("Multiple Logistic Regression - Example 8.2") [1] "Multiple Logistic Regression - Example 8.2" > > plot(northing ~ easting, data=frogs, pch=c(1,16)[frogs$pres.abs+1], + xlab="Meters east of reference point", ylab="Meters north") > > pause() Pause. Press to continue... > pairs(frogs[,4:10], oma=c(2,2,2,2), cex=0.5) > > pause() Pause. Press to continue... > oldpar <- par(mfrow=c(1,3)) > for(nam in c("distance","NoOfPools","NoOfSites")){ + y <- frogs[,nam] + plot(density(y),main="",xlab=nam) + par(oldpar) + } > > pause() Pause. Press to continue... > attach(frogs) > pairs(cbind(altitude,log(distance),log(NoOfPools),NoOfSites), + panel=panel.smooth, labels=c("altitude","log(distance)", + "log(NoOfPools)","NoOfSites")) > detach(frogs) > > frogs.glm0 <- glm(formula = pres.abs ~ altitude + log(distance) + + log(NoOfPools) + NoOfSites + avrain + meanmin + meanmax, + family = binomial, data = frogs) > summary(frogs.glm0) Call: glm(formula = pres.abs ~ altitude + log(distance) + log(NoOfPools) + NoOfSites + avrain + meanmin + meanmax, family = binomial, data = frogs) Deviance Residuals: Min 1Q Median 3Q Max -1.979 -0.719 -0.278 0.796 2.566 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 4.09e+01 1.33e+02 0.31 0.75785 altitude -6.65e-03 3.87e-02 -0.17 0.86347 log(distance) -7.59e-01 2.55e-01 -2.97 0.00295 ** log(NoOfPools) 5.73e-01 2.16e-01 2.65 0.00808 ** NoOfSites -8.98e-04 1.07e-01 -0.01 0.99333 avrain -6.79e-03 6.00e-02 -0.11 0.90985 meanmin 5.30e+00 1.54e+00 3.44 0.00058 *** meanmax -3.17e+00 4.84e+00 -0.66 0.51205 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 279.99 on 211 degrees of freedom Residual deviance: 197.62 on 204 degrees of freedom AIC: 213.6 Number of Fisher Scoring iterations: 5 > pause() Pause. Press to continue... > frogs.glm <- glm(formula = pres.abs ~ log(distance) + log(NoOfPools) + + meanmin + + meanmax, family = binomial, data = frogs) > oldpar <- par(mfrow=c(2,2)) > termplot(frogs.glm, data=frogs) > par(oldpar) > pause() Pause. Press to continue... > termplot(frogs.glm, data=frogs, partial.resid=TRUE) > > cv.binary(frogs.glm0) # All explanatory variables Fold: 3 4 6 10 9 7 1 2 8 5 Internal estimate of accuracy = 0 Cross-validation estimate of accuracy = 0 > pause() Pause. Press to continue... > cv.binary(frogs.glm) # Reduced set of explanatory variables Fold: 8 10 9 4 7 3 2 6 1 5 Internal estimate of accuracy = 0 Cross-validation estimate of accuracy = 0 > > pause() Pause. Press to continue... > for (j in 1:4){ + rand <- sample(1:10, 212, replace=TRUE) + all.acc <- cv.binary(frogs.glm0, rand=rand, print.details=FALSE)$acc.cv + reduced.acc <- cv.binary(frogs.glm, rand=rand, print.details=FALSE)$acc.cv + cat("\nAll:", round(all.acc,3), " Reduced:", round(reduced.acc,3)) + } All: 0 Reduced: 0 All: 0 Reduced: 0 All: 0 Reduced: 0 All: 0 Reduced: 0> > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "fruitohms" > > ### * fruitohms > > flush(stderr()); flush(stdout()) > > ### Name: fruitohms > ### Title: Electrical Resistance of Kiwi Fruit > ### Aliases: fruitohms > ### Keywords: datasets > > ### ** Examples > > plot(ohms ~ juice, xlab="Apparent juice content (%)",ylab="Resistance (ohms)", data=fruitohms) > lines(lowess(fruitohms$juice, fruitohms$ohms), lwd=2) > pause() Pause. Press to continue... > require(splines) Loading required package: splines [1] TRUE > attach(fruitohms) > plot(ohms ~ juice, cex=0.8, xlab="Apparent juice content (%)", + ylab="Resistance (ohms)", type="n") > fruit.lmb4 <- lm(ohms ~ bs(juice,4)) > ord <- order(juice) > lines(juice[ord], fitted(fruit.lmb4)[ord], lwd=2) > ci <- predict(fruit.lmb4, interval="confidence") > lines(juice[ord], ci[ord,"lwr"]) > lines(juice[ord], ci[ord,"upr"]) > > > > cleanEx(); ..nameEx <- "geophones" > > ### * geophones > > flush(stderr()); flush(stdout()) > > ### Name: geophones > ### Title: Seismic Timing Data > ### Aliases: geophones > ### Keywords: datasets > > ### ** Examples > > plot(geophones) > lines(lowess(geophones, f=.25)) > > > > cleanEx(); ..nameEx <- "hills" > > ### * hills > > flush(stderr()); flush(stdout()) > > ### Name: hills > ### Title: Scottish Hill Races Data > ### Aliases: hills > ### Keywords: datasets > > ### ** Examples > > print("Transformation - Example 6.4.3") [1] "Transformation - Example 6.4.3" > pairs(hills, labels=c("dist\n\n(miles)", "climb\n\n(feet)", + "time\n\n(hours)")) > pause() Pause. Press to continue... > pairs(log(hills), labels=c("dist\n\n(log(miles))", "climb\n\n(log(feet))", + "time\n\n(log(hours))")) > pause() Pause. Press to continue... > hills0.loglm <- lm(log(time) ~ log(dist) + log(climb), data = hills) > oldpar <- par(mfrow=c(2,2)) > plot(hills0.loglm) > pause() Pause. Press to continue... > hills.loglm <- lm(log(time) ~ log(dist) + log(climb), data = hills[-18,]) > summary(hills.loglm) Call: lm(formula = log(time) ~ log(dist) + log(climb), data = hills[-18, ]) Residuals: Min 1Q Median 3Q Max -0.51957 -0.07342 0.00735 0.07345 0.32893 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.8816 0.2826 -13.73 1.0e-14 *** log(dist) 0.9092 0.0650 13.99 6.2e-15 *** log(climb) 0.2601 0.0484 5.37 7.3e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.16 on 31 degrees of freedom Multiple R-Squared: 0.952, Adjusted R-squared: 0.949 F-statistic: 306 on 2 and 31 DF, p-value: <2e-16 > plot(hills.loglm) > pause() Pause. Press to continue... > hills2.loglm <- lm(log(time) ~ log(dist)+log(climb)+log(dist):log(climb), + data=hills[-18,]) > anova(hills.loglm, hills2.loglm) Analysis of Variance Table Model 1: log(time) ~ log(dist) + log(climb) Model 2: log(time) ~ log(dist) + log(climb) + log(dist):log(climb) Res.Df RSS Df Sum of Sq F Pr(>F) 1 31 0.798 2 30 0.741 1 0.057 2.31 0.14 > pause() Pause. Press to continue... > step(hills2.loglm) Start: AIC= -122 log(time) ~ log(dist) + log(climb) + log(dist):log(climb) Df Sum of Sq RSS AIC 0.7 -122.1 - log(dist):log(climb) 1 0.1 0.8 -121.6 Call: lm(formula = log(time) ~ log(dist) + log(climb) + log(dist):log(climb), data = hills[-18, ]) Coefficients: (Intercept) log(dist) log(climb) -2.4755 0.1685 0.0672 log(dist):log(climb) 0.0993 > pause() Pause. Press to continue... > summary(hills.loglm, corr=TRUE)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) -3.88 0.2826 -13.73 1.01e-14 log(dist) 0.91 0.0650 13.99 6.16e-15 log(climb) 0.26 0.0484 5.37 7.33e-06 > pause() Pause. Press to continue... > summary(hills2.loglm, corr=TRUE)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) -2.4755 0.9653 -2.565 0.0156 log(dist) 0.1685 0.4913 0.343 0.7339 log(climb) 0.0672 0.1354 0.497 0.6231 log(dist):log(climb) 0.0993 0.0653 1.520 0.1389 > par(oldpar) > pause() Pause. Press to continue... > print("Nonlinear - Example 6.9.4") [1] "Nonlinear - Example 6.9.4" > require(nls) Loading required package: nls Warning: package 'nls' has been merged into 'stats' [1] TRUE > hills.nls0 <- nls(time ~ (dist^alpha)*(climb^beta), start = + c(alpha = .909, beta = .260), data = hills[-18,]) > summary(hills.nls0) Formula: time ~ (dist^alpha) * (climb^beta) Parameters: Estimate Std. Error t value Pr(>|t|) alpha 1.0573 0.1261 8.39 1.4e-09 *** beta -0.2849 0.0436 -6.53 2.4e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.506 on 32 degrees of freedom Correlation of Parameter Estimates: alpha beta -0.976 > plot(residuals(hills.nls0) ~ predict(hills.nls0)) # residual plot > pause() Pause. Press to continue... > hills$climb.mi <- hills$climb/5280 > hills.nls <- nls(time ~ alpha + beta*dist + gamma*(climb.mi^delta), + start=c(alpha = 1, beta = 1, gamma = 1, delta = 1), data=hills[-18,]) > summary(hills.nls) Formula: time ~ alpha + beta * dist + gamma * (climb.mi^delta) Parameters: Estimate Std. Error t value Pr(>|t|) alpha -0.05310 0.02947 -1.80 0.082 . beta 0.11076 0.00398 27.82 < 2e-16 *** gamma 0.77277 0.07523 10.27 2.4e-11 *** delta 2.18414 0.22706 9.62 1.1e-10 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0953 on 30 degrees of freedom Correlation of Parameter Estimates: alpha beta gamma beta -0.5794 gamma -0.0592 -0.650 delta 0.2308 0.427 -0.773 > plot(residuals(hills.nls) ~ predict(hills.nls)) # residual plot > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "houseprices" > > ### * houseprices > > flush(stderr()); flush(stdout()) > > ### Name: houseprices > ### Title: Aranda House Prices > ### Aliases: houseprices > ### Keywords: datasets > > ### ** Examples > > plot(sale.price~area, data=houseprices) > pause() Pause. Press to continue... > coplot(sale.price~area|bedrooms, data=houseprices) > pause() Pause. Press to continue... > print("Cross-Validation - Example 5.5.2") [1] "Cross-Validation - Example 5.5.2" > > houseprices.lm <- lm(sale.price ~ area, data=houseprices) > summary(houseprices.lm)$sigma^2 [1] 2321 > pause() Pause. Press to continue... > cv.lm() Analysis of Variance Table Response: yv Df Sum Sq Mean Sq F value Pr(>F) xv 1 18566 18566 8 0.014 * Residuals 13 30179 2321 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 fold 1 Observations in test set: 3 12 13 14 15 X1 X2 X3 X4 X5 area 802.0 696 771.0 1006.0 1191 Predicted 204.0 188 199.3 234.7 262 sale.price 215.0 255 260.0 293.0 375 Residual 11.0 67 60.7 58.3 113 Sum of squares = 24351 Mean square = 4870 n = 5 fold 2 Observations in test set: 2 5 6 9 10 X1 X2 X3 X4 X5 area 905 716 963.0 1018.00 887.00 Predicted 255 224 264.4 273.38 252.06 sale.price 215 113 185.0 276.00 260.00 Residual -40 -112 -79.4 2.62 7.94 Sum of squares = 20416 Mean square = 4083 n = 5 fold 3 Observations in test set: 1 4 7 8 11 X1 X2 X3 X4 X5 area 694.0 1366 821.00 714.0 790.00 Predicted 183.2 388 221.94 189.3 212.49 sale.price 192.0 274 212.00 220.0 221.50 Residual 8.8 -114 -9.94 30.7 9.01 Sum of squares = 14241 Mean square = 2848 n = 5 Overall ms 3934 > pause() Pause. Press to continue... > print("Bootstrapping - Example 5.5.3") [1] "Bootstrapping - Example 5.5.3" > houseprices.fn <- function (houseprices, index){ + house.resample <- houseprices[index,] + house.lm <- lm(sale.price ~ area, data=house.resample) + coef(house.lm)[2] # slope estimate for resampled data + } > require(boot) # ensure that the boot package is loaded Loading required package: boot [1] TRUE > houseprices.boot <- boot(houseprices, R=999, statistic=houseprices.fn) > > houseprices1.fn <- function (houseprices, index){ + house.resample <- houseprices[index,] + house.lm <- lm(sale.price ~ area, data=house.resample) + predict(house.lm, newdata=data.frame(area=1200)) + } > > houseprices1.boot <- boot(houseprices, R=999, statistic=houseprices1.fn) > boot.ci(houseprices1.boot, type="perc") # "basic" is an alternative to "perc" BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 999 bootstrap replicates CALL : boot.ci(boot.out = houseprices1.boot, type = "perc") Intervals : Level Percentile 95% (244, 365 ) Calculations and Intervals on Original Scale > houseprices2.fn <- function (houseprices, index){ + house.resample <- houseprices[index,] + house.lm <- lm(sale.price ~ area, data=house.resample) + houseprices$sale.price-predict(house.lm, houseprices) # resampled prediction errors + } > > n <- length(houseprices$area) > R <- 200 > houseprices2.boot <- boot(houseprices, R=R, statistic=houseprices2.fn) > house.fac <- factor(rep(1:n, rep(R, n))) > plot(house.fac, as.vector(houseprices2.boot$t), ylab="Prediction Errors", + xlab="House") > pause() Pause. Press to continue... > plot(apply(houseprices2.boot$t,2, sd)/predict.lm(houseprices.lm, se.fit=TRUE)$se.fit, + ylab="Ratio of Bootstrap SE's to Model-Based SE's", xlab="House", pch=16) > abline(1,0) > > > > > cleanEx(); ..nameEx <- "ironslag" > > ### * ironslag > > flush(stderr()); flush(stdout()) > > ### Name: ironslag > ### Title: Iron Content Measurements > ### Aliases: ironslag > ### Keywords: datasets > > ### ** Examples > > iron.lm <- lm(chemical ~ magnetic, data = ironslag) > oldpar <- par(mfrow = c(2,2)) > plot(iron.lm) > par(oldpar) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "jobs" > > ### * jobs > > flush(stderr()); flush(stdout()) > > ### Name: jobs > ### Title: Canadian Labour Force Summary Data (1995-96) > ### Aliases: jobs > ### Keywords: datasets > > ### ** Examples > > print("Multiple Variables and Times - Example 2.1.4") [1] "Multiple Variables and Times - Example 2.1.4" > sapply(jobs, range) BC Alberta Prairies Ontario Quebec Atlantic Date [1,] 1737 1366 973 5212 3167 941 95 [2,] 1840 1436 999 5360 3257 968 97 > pause() Pause. Press to continue... > matplot(jobs[,7], jobs[,-7], type="l", xlim=c(95,97.1)) > # Notice that we have been able to use a data frame as the second argument to matplot(). > # For more information on matplot(), type help(matplot) > text(rep(jobs[24,7], 6), jobs[24,1:6], names(jobs)[1:6], adj=0) > pause() Pause. Press to continue... > sapply(log(jobs[,-7]), range) BC Alberta Prairies Ontario Quebec Atlantic [1,] 7.46 7.22 6.88 8.56 8.06 6.85 [2,] 7.52 7.27 6.91 8.59 8.09 6.88 > apply(sapply(log(jobs[,-7]), range), 2, diff) BC Alberta Prairies Ontario Quebec Atlantic 0.0576 0.0500 0.0264 0.0280 0.0280 0.0283 > pause() Pause. Press to continue... > oldpar <- par(mfrow=c(2,3)) > range.log <- sapply(log(jobs[,-7], 2), range) > maxdiff <- max(apply(range.log, 2, diff)) > range.log[2,] <- range.log[1,] + maxdiff > titles <- c("BC Jobs","Alberta Jobs","Prairie Jobs", + "Ontario Jobs", "Quebec Jobs", "Atlantic Jobs") > for (i in 1:6){ + plot(jobs$Date, log(jobs[,i], 2), type = "l", ylim = range.log[,i], + xlab = "Time", ylab = "Number of jobs", main = titles[i]) + } > par(oldpar) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "kiwishade" > > ### * kiwishade > > flush(stderr()); flush(stdout()) > > ### Name: kiwishade > ### Title: Kiwi Shading Data > ### Aliases: kiwishade > ### Keywords: datasets > > ### ** Examples > > print("Data Summary - Example 2.2.1") [1] "Data Summary - Example 2.2.1" > attach(kiwishade) > kiwimeans <- aggregate(yield, by=list(block, shade), mean) > names(kiwimeans) <- c("block","shade","meanyield") > > kiwimeans[1:4,] block shade meanyield 1 east none 99.0 2 north none 104.0 3 west none 97.6 4 east Aug2Dec 105.6 > pause() Pause. Press to continue... > print("Multilevel Design - Example 9.3") [1] "Multilevel Design - Example 9.3" > kiwishade.aov <- aov(yield ~ shade+Error(block/shade),data=kiwishade) > summary(kiwishade.aov) Error: block Df Sum Sq Mean Sq F value Pr(>F) Residuals 2 172.3 86.2 Error: block:shade Df Sum Sq Mean Sq F value Pr(>F) shade 3 1395 465 22.2 0.0012 ** Residuals 6 126 21 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 36 439 12 > pause() Pause. Press to continue... > sapply(split(yield, shade), mean) none Aug2Dec Dec2Feb Feb2May 100.2 103.2 89.9 92.8 > > pause() Pause. Press to continue... > kiwi.table <- t(sapply(split(yield, plot), as.vector)) > kiwi.means <- sapply(split(yield, plot), mean) > kiwi.means.table <- matrix(rep(kiwi.means,4), nrow=12, ncol=4) > kiwi.summary <- data.frame(kiwi.means, kiwi.table-kiwi.means.table) > names(kiwi.summary)<- c("Mean", "Vine 1", "Vine 2", "Vine 3", "Vine 4") > kiwi.summary Mean Vine 1 Vine 2 Vine 3 Vine 4 east.Aug2Dec 105.6 3.315 -2.415 2.765 -3.6650 east.Dec2Feb 88.9 -3.520 -7.490 5.020 5.9900 east.Feb2May 95.9 -3.430 1.690 1.270 0.4700 east.none 99.0 1.285 -2.025 -0.975 1.7150 north.Aug2Dec 103.6 -0.495 -1.075 0.425 1.1450 north.Dec2Feb 93.7 -1.528 0.112 1.062 0.3525 north.Feb2May 92.0 0.545 1.415 -1.745 -0.2150 north.none 104.0 -2.915 3.995 2.645 -3.7250 west.Aug2Dec 100.5 5.220 0.860 -2.950 -3.1300 west.Dec2Feb 87.1 -0.495 3.605 -0.085 -3.0250 west.Feb2May 90.4 -3.138 4.582 -1.347 -0.0975 west.none 97.6 -4.918 5.932 3.123 -4.1375 > mean(kiwi.means) # the grand mean (only for balanced design) [1] 96.5 > > require(nlme) Loading required package: nlme [1] TRUE > kiwishade.lme <- lme(fixed = yield ~ shade, random = ~ 1 | block/plot, + data=kiwishade) > res <- residuals(kiwishade.lme) > hat <- fitted(kiwishade.lme) # By default fitted(kiwishade.lme, level=2) > coplot(res ~ hat | kiwishade$block, pch=16, columns=3, + xlab= "Fitted", ylab="Residuals") > > res <- residuals(kiwishade.lme) > hat <- fitted(kiwishade.lme, level=0) # shade effects only > unique(hat) # There are just four distinct values, one per treatment [1] 100.2 103.2 89.9 92.8 > coplot(res ~ hat | kiwishade$block, pch=16, columns=3, + xlab="Fitted", ylab="Residuals") > > n.omit <- 2 > take <- rep(TRUE, 48) > take[sample(1:48,2)] <- FALSE > kiwishade.lme <- lme(yield ~ shade, data = kiwishade, + random = ~1 | block/plot, subset=take) > VarCorr(kiwishade.lme)[4, 1] # Plot component of variance [1] " 2.51" > VarCorr(kiwishade.lme)[4, 1] # Vine component of variance [1] " 2.51" > > detach(kiwishade) > > > > > cleanEx(); ..nameEx <- "leafshape17" > > ### * leafshape17 > > flush(stderr()); flush(stdout()) > > ### Name: leafshape17 > ### Title: Subset of Leaf Shape Data Set > ### Aliases: leafshape17 > ### Keywords: datasets > > ### ** Examples > > print("Discriminant Analysis - Example 11.2") [1] "Discriminant Analysis - Example 11.2" > > require(MASS) Loading required package: MASS Attaching package: 'MASS' The following object(s) are masked from package:DAAG : hills [1] TRUE > leaf17.lda <- lda(arch ~ logwid+loglen, data=leafshape17) > leaf17.hat <- predict(leaf17.lda) > leaf17.lda Call: lda(arch ~ logwid + loglen, data = leafshape17) Prior probabilities of groups: 0 1 0.672 0.328 Group means: logwid loglen 0 1.43 2.46 1 1.87 2.99 Coefficients of linear discriminants: LD1 logwid 0.156 loglen 3.066 > table(leafshape17$arch, leaf17.hat$class) 0 1 0 38 3 1 8 12 > pause() Pause. Press to continue... > tab <- table(leafshape17$arch, leaf17.hat$class) > sum(tab[row(tab)==col(tab)])/sum(tab) [1] 0.82 > leaf17cv.lda <- lda(arch ~ logwid+loglen, data=leafshape17, CV=TRUE) > tab <- table(leafshape17$arch, leaf17cv.lda$class) > pause() Pause. Press to continue... > leaf17.glm <- glm(arch ~ logwid + loglen, family=binomial, data=leafshape17) > options(digits=3) > summary(leaf17.glm)$coef Estimate Std. Error z value Pr(>|z|) (Intercept) -15.287 4.13 -3.703 0.000213 logwid 0.185 1.57 0.118 0.906224 loglen 5.269 1.96 2.684 0.007269 > pause() Pause. Press to continue... > leaf17.one <- cv.binary(leaf17.glm) Fold: 3 4 6 10 9 7 1 2 8 5 Internal estimate of accuracy = 0 Cross-validation estimate of accuracy = 0 > table(leafshape17$arch, round(leaf17.one$internal)) # Resubstitution 0 1 0 37 4 1 7 13 > pause() Pause. Press to continue... > table(leafshape17$arch, round(leaf17.one$cv)) # Cross-validation 0 1 0 36 5 1 9 11 > > > > cleanEx(); ..nameEx <- "leaftemp" > > ### * leaftemp > > flush(stderr()); flush(stdout()) > > ### Name: leaftemp > ### Title: Leaf and Air Temperature Data > ### Aliases: leaftemp > ### Keywords: datasets > > ### ** Examples > > print("Fitting Multiple Lines - Example 7.3") [1] "Fitting Multiple Lines - Example 7.3" > > leaf.lm1 <- lm(tempDiff ~ 1 , data = leaftemp) > leaf.lm2 <- lm(tempDiff ~ vapPress, data = leaftemp) > leaf.lm3 <- lm(tempDiff ~ CO2level + vapPress, data = leaftemp) > leaf.lm4 <- lm(tempDiff ~ CO2level + vapPress + vapPress:CO2level, + data = leaftemp) > > anova(leaf.lm1, leaf.lm2, leaf.lm3, leaf.lm4) Analysis of Variance Table Model 1: tempDiff ~ 1 Model 2: tempDiff ~ vapPress Model 3: tempDiff ~ CO2level + vapPress Model 4: tempDiff ~ CO2level + vapPress + vapPress:CO2level Res.Df RSS Df Sum of Sq F Pr(>F) 1 61 40.0 2 60 34.7 1 5.3 11.33 0.0014 ** 3 58 28.2 2 6.5 7.03 0.0019 ** 4 56 26.1 2 2.1 2.28 0.1112 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > summary(leaf.lm2) Call: lm(formula = tempDiff ~ vapPress, data = leaftemp) Residuals: Min 1Q Median 3Q Max -1.3680 -0.6119 0.0823 0.5040 1.7733 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.097 0.584 5.31 1.7e-06 *** vapPress -0.857 0.284 -3.02 0.0037 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.761 on 60 degrees of freedom Multiple R-Squared: 0.132, Adjusted R-squared: 0.117 F-statistic: 9.11 on 1 and 60 DF, p-value: 0.00373 > plot(leaf.lm2) > > > > > cleanEx(); ..nameEx <- "litters" > > ### * litters > > flush(stderr()); flush(stdout()) > > ### Name: litters > ### Title: Mouse Litters > ### Aliases: litters > ### Keywords: datasets > > ### ** Examples > > print("Multiple Regression - Example 6.2") [1] "Multiple Regression - Example 6.2" > > pairs(litters, labels=c("lsize\n\n(litter size)", "bodywt\n\n(Body Weight)", "brainwt\n\n(Brain Weight)")) > # pairs(litters) gives a scatterplot matrix with less adequate labeling > > mice1.lm <- lm(brainwt ~ lsize, data = litters) # Regress on lsize > mice2.lm <- lm(brainwt ~ bodywt, data = litters) #Regress on bodywt > mice12.lm <- lm(brainwt ~ lsize + bodywt, data = litters) # Regress on lsize & bodywt > > summary(mice1.lm)$coef # Similarly for other coefficients. Estimate Std. Error t value Pr(>|t|) (Intercept) 0.44700 0.00962 46.44 3.39e-20 lsize -0.00403 0.00120 -3.37 3.44e-03 > # results are consistent with the biological concept of brain sparing > > pause() Pause. Press to continue... > hat(model.matrix(mice12.lm)) # hat diagonal [1] 0.1991 0.1734 0.1303 0.1574 0.0879 0.2732 0.0677 0.0693 0.1555 0.0947 [11] 0.1279 0.0757 0.1397 0.0733 0.1251 0.0879 0.4326 0.1259 0.2042 0.1993 > pause() Pause. Press to continue... > plot(lm.influence(mice12.lm)$hat, residuals(mice12.lm)) > > print("Diagnostics - Example 6.3") [1] "Diagnostics - Example 6.3" > > mice12.lm <- lm(brainwt ~ bodywt+lsize, data=litters) > oldpar <-par(mfrow = c(1,2)) > bx <- mice12.lm$coef[2]; bz <- mice12.lm$coef[3] > res <- residuals(mice12.lm) > plot(litters$bodywt, bx*litters$bodywt+res, xlab="Body weight", + ylab="Component + Residual") > panel.smooth(litters$bodywt, bx*litters$bodywt+res) # Overlay > plot(litters$lsize, bz*litters$lsize+res, xlab="Litter size", + ylab="Component + Residual") > panel.smooth(litters$lsize, bz*litters$lsize+res) > par(oldpar) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "mifem" > > ### * mifem > > flush(stderr()); flush(stdout()) > > ### Name: mifem > ### Title: Mortality Outcomes for Females Suffering Myocardial Infarction > ### Aliases: mifem > ### Keywords: datasets > > ### ** Examples > > print("CART - Example 10.7") [1] "CART - Example 10.7" > summary(mifem) outcome age yronset premi smstat diabetes highbp live:974 Min. :35.0 Min. :85.0 y :311 c :390 y :248 y :813 dead:321 1st Qu.:57.0 1st Qu.:87.0 n :928 x :280 n :978 n :406 Median :63.0 Median :89.0 nk: 56 n :522 nk: 69 nk: 76 Mean :60.9 Mean :88.8 nk:103 3rd Qu.:66.0 3rd Qu.:91.0 Max. :69.0 Max. :93.0 hichol angina stroke y :452 y :472 y : 153 n :655 n :724 n :1063 nk:188 nk: 99 nk: 79 > pause() Pause. Press to continue... > require(rpart) Loading required package: rpart [1] TRUE > mifem.rpart <- rpart(outcome ~ ., data = mifem, cp = 0.0025) > plotcp(mifem.rpart) > printcp(mifem.rpart) Classification tree: rpart(formula = outcome ~ ., data = mifem, cp = 0.0025) Variables actually used in tree construction: [1] age angina diabetes hichol premi smstat stroke yronset Root node error: 321/1295 = 0.2 n= 1295 CP nsplit rel error xerror xstd 1 0.202 0 1.0 1.0 0.05 2 0.006 1 0.8 0.8 0.04 3 0.005 13 0.7 0.8 0.05 4 0.003 17 0.7 0.8 0.05 5 0.003 18 0.7 0.8 0.05 > pause() Pause. Press to continue... > mifemb.rpart <- prune(mifem.rpart, cp=0.006) > print(mifemb.rpart) n= 1295 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 1295 321 live (0.752 0.248) 2) angina=y,n 1196 239 live (0.800 0.200) * 3) angina=nk 99 17 dead (0.172 0.828) * > > > > cleanEx(); ..nameEx <- "mignonette" > > ### * mignonette > > flush(stderr()); flush(stdout()) > > ### Name: mignonette > ### Title: Darwin's Wild Mignonette Data > ### Aliases: mignonette > ### Keywords: datasets > > ### ** Examples > > print("Is Pairing Helpful? - Example 4.3.1") [1] "Is Pairing Helpful? - Example 4.3.1" > > attach(mignonette) > plot(cross ~ self, pch=rep(c(4,1), c(3,12))); abline(0,1) > abline(mean(cross-self), 1, lty=2) > detach(mignonette) > > > > cleanEx(); ..nameEx <- "milk" > > ### * milk > > flush(stderr()); flush(stdout()) > > ### Name: milk > ### Title: Milk Sweetness Study > ### Aliases: milk > ### Keywords: datasets > > ### ** Examples > > print("Rug Plot - Example 1.8.1") [1] "Rug Plot - Example 1.8.1" > xyrange <- range(milk) > plot(four ~ one, data = milk, xlim = xyrange, ylim = xyrange, pch = 16) > rug(milk$one) > rug(milk$four, side = 2) > abline(0, 1) > > > > cleanEx(); ..nameEx <- "modelcars" > > ### * modelcars > > flush(stderr()); flush(stdout()) > > ### Name: modelcars > ### Title: Model Car Data > ### Aliases: modelcars > ### Keywords: datasets > > ### ** Examples > > plot(modelcars) > modelcars.lm <- lm(distance.traveled ~ starting.point, data=modelcars) > aov(modelcars.lm) Call: aov(formula = modelcars.lm) Terms: starting.point Residuals Sum of Squares 548 23 Deg. of Freedom 1 10 Residual standard error: 1.52 Estimated effects may be unbalanced > pause() Pause. Press to continue... > print("Response Curves - Example 4.6") [1] "Response Curves - Example 4.6" > attach(modelcars) > stripchart(distance.traveled ~ starting.point, vertical=TRUE, pch=15, xlab = "Distance up ramp", ylab="Distance traveled") > detach(modelcars) > > > > > cleanEx(); ..nameEx <- "monica" > > ### * monica > > flush(stderr()); flush(stdout()) > > ### Name: monica > ### Title: WHO Monica Data > ### Aliases: monica > ### Keywords: datasets > > ### ** Examples > > print("CART - Example 10.7") [1] "CART - Example 10.7" > summary(monica) outcome sex age yronset premi smstat live:3525 m:4605 Min. :35.0 Min. :85.0 y :1511 c :2051 dead:2842 f:1762 1st Qu.:55.0 1st Qu.:87.0 n :4122 x :1938 Median :61.0 Median :89.0 nk: 734 n :1460 Mean :59.4 Mean :88.7 nk: 918 3rd Qu.:66.0 3rd Qu.:91.0 Max. :69.0 Max. :93.0 diabetes highbp hichol angina stroke hosp y : 818 y :2877 y :1840 y :1919 y : 560 y:4442 n :4664 n :2542 n :3294 n :3473 n :4881 n:1925 nk: 885 nk: 948 nk:1233 nk: 975 nk: 926 > pause() Pause. Press to continue... > require(rpart) Loading required package: rpart [1] TRUE > monica.rpart <- rpart(outcome ~ ., data = monica, cp = 0.0025) > plotcp(monica.rpart) > printcp(monica.rpart) Classification tree: rpart(formula = outcome ~ ., data = monica, cp = 0.0025) Variables actually used in tree construction: [1] highbp hosp Root node error: 2842/6367 = 0.4 n= 6367 CP nsplit rel error xerror xstd 1 0.675 0 1.0 1.0 0.01 2 0.052 1 0.3 0.3 0.01 3 0.003 2 0.3 0.3 0.01 > pause() Pause. Press to continue... > monicab.rpart <- prune(monica.rpart, cp=0.006) > print(monicab.rpart) n= 6367 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 6367 2840 live (0.55364 0.44636) 2) hosp=y 4442 920 live (0.79289 0.20711) 4) highbp=y,n 4199 724 live (0.82758 0.17242) * 5) highbp=nk 243 47 dead (0.19342 0.80658) * 3) hosp=n 1925 3 dead (0.00156 0.99844) * > > > > cleanEx(); ..nameEx <- "moths" > > ### * moths > > flush(stderr()); flush(stdout()) > > ### Name: moths > ### Title: Moths Data > ### Aliases: moths > ### Keywords: datasets > > ### ** Examples > > print("Quasi Poisson Regression - Example 8.3") [1] "Quasi Poisson Regression - Example 8.3" > rbind(table(moths[,4]), sapply(split(moths[,-4], moths$habitat), apply,2, + sum)) Bank Disturbed Lowerside NEsoak NWsoak SEsoak SWsoak Upperside 1 7 9 6 3 7 3 5 meters 21 49 191 254 65 193 116 952 A 0 8 41 14 71 37 20 28 P 4 33 17 14 19 6 48 8 > A.glm <- glm(formula = A ~ log(meters) + factor(habitat), family = + quasipoisson, data = moths) > summary(A.glm) Call: glm(formula = A ~ log(meters) + factor(habitat), family = quasipoisson, data = moths) Deviance Residuals: Min 1Q Median 3Q Max -3.556190 -1.463482 -0.000672 0.924806 3.082152 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -15.696 2096.588 -0.01 0.99 log(meters) 0.129 0.153 0.85 0.40 factor(habitat)Disturbed 15.622 2096.588 0.01 0.99 factor(habitat)Lowerside 16.906 2096.588 0.01 0.99 factor(habitat)NEsoak 16.084 2096.588 0.01 0.99 factor(habitat)NWsoak 18.468 2096.588 0.01 0.99 factor(habitat)SEsoak 16.968 2096.588 0.01 0.99 factor(habitat)SWsoak 17.137 2096.588 0.01 0.99 factor(habitat)Upperside 16.743 2096.588 0.01 0.99 (Dispersion parameter for quasipoisson family taken to be 2.7) Null deviance: 257.108 on 40 degrees of freedom Residual deviance: 93.991 on 32 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 13 > moths$habitat <- relevel(moths$habitat, ref="Lowerside") > A.glm <- glm(A ~ habitat + log(meters), family=quasipoisson, data=moths) > summary(A.glm)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 1.2098 0.455 2.65718 1.22e-02 habitatBank -16.9057 2096.588 -0.00806 9.94e-01 habitatDisturbed -1.2832 0.647 -1.98223 5.61e-02 habitatNEsoak -0.8217 0.537 -1.53016 1.36e-01 habitatNWsoak 1.5624 0.334 4.67424 5.11e-05 habitatSEsoak 0.0621 0.385 0.16114 8.73e-01 habitatSWsoak 0.2314 0.478 0.48410 6.32e-01 habitatUpperside -0.1623 0.584 -0.27784 7.83e-01 log(meters) 0.1292 0.153 0.84544 4.04e-01 > > > > cleanEx(); ..nameEx <- "nsw74psid1" > > ### * nsw74psid1 > > flush(stderr()); flush(stdout()) > > ### Name: nsw74psid1 > ### Title: Labour Training Evaluation Data > ### Aliases: nsw74psid1 > ### Keywords: datasets > > ### ** Examples > > print("Interpretation of Regression Coefficients - Example 6.6") [1] "Interpretation of Regression Coefficients - Example 6.6" > > nsw74psid1.lm <- lm(re78~ trt+ (age + educ + re74 + re75) + + (black + hisp + marr + nodeg), data = nsw74psid1) > summary(nsw74psid1.lm)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) -129.743 1.69e+03 -0.0768 9.39e-01 trt 751.946 9.15e+02 0.8216 4.11e-01 age -83.566 2.08e+01 -4.0149 6.11e-05 educ 592.610 1.03e+02 5.7366 1.07e-08 re74 0.278 2.79e-02 9.9598 5.73e-23 re75 0.568 2.76e-02 20.6130 1.02e-87 black -570.928 4.95e+02 -1.1530 2.49e-01 hisp 2163.281 1.09e+03 1.9805 4.78e-02 marr 1240.520 5.86e+02 2.1160 3.44e-02 nodeg 590.467 6.47e+02 0.9129 3.61e-01 > options(digits=4) > sapply(nsw74psid1[, c(2,3,8,9,10)], quantile, prob=c(.25,.5,.75,.95,1)) age educ re74 re75 re78 25% 25.0 10 8817 7605 9243 50% 32.0 12 17437 17008 19432 75% 43.5 14 25470 25584 28816 95% 53.0 17 41145 41177 45810 100% 55.0 17 137149 156653 121174 > attach(nsw74psid1) > sapply(nsw74psid1[trt==1, c(2,3,8,9,10)], quantile, + prob=c(.25,.5,.75,.95,1)) age educ re74 re75 re78 25% 20 9 0 0 485.2 50% 25 11 0 0 4232.3 75% 29 12 1291 1817 9643.0 95% 42 13 12342 6855 18774.7 100% 48 16 35040 25142 60307.9 > pause() Pause. Press to continue... > here <- age <= 40 & re74<=5000 & re75 <= 5000 & re78 < 30000 > nsw74psidA <- nsw74psid1[here, ] > detach(nsw74psid1) > table(nsw74psidA$trt) 0 1 119 133 > pause() Pause. Press to continue... > A1.lm <- lm(re78 ~ trt + (age + educ + re74 + re75) + (black + + hisp + marr + nodeg), data = nsw74psidA) > summary(A1.lm)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) -2093.4414 4253.8344 -0.49213 0.62307 trt 2942.9665 1169.9433 2.51548 0.01254 age 71.2005 77.1393 0.92301 0.35692 educ 245.9885 250.6877 0.98125 0.32745 re74 0.3140 0.4002 0.78455 0.43348 re75 0.8717 0.3826 2.27799 0.02360 black 905.5948 1135.1582 0.79777 0.42579 hisp 3723.0383 2211.8204 1.68325 0.09362 marr 24.7218 1104.6055 0.02238 0.98216 nodeg -883.4520 1294.2084 -0.68262 0.49550 > pause() Pause. Press to continue... > A2.lm <- lm(re78 ~ trt + (age + educ + re74 + re75) * (black + + hisp + marr + nodeg), data = nsw74psidA) > anova(A1.lm, A2.lm) Analysis of Variance Table Model 1: re78 ~ trt + (age + educ + re74 + re75) + (black + hisp + marr + nodeg) Model 2: re78 ~ trt + (age + educ + re74 + re75) * (black + hisp + marr + nodeg) Res.Df RSS Df Sum of Sq F Pr(>F) 1 242 9.72e+09 2 226 9.08e+09 16 6.47e+08 1.01 0.45 > > > > > cleanEx(); ..nameEx <- "nsw74psid3" > > ### * nsw74psid3 > > flush(stderr()); flush(stdout()) > > ### Name: nsw74psid3 > ### Title: Labour Training Evaluation Data > ### Aliases: nsw74psid3 > ### Keywords: datasets > > ### ** Examples > > print("Contingency Tables - Example 4.4") [1] "Contingency Tables - Example 4.4" > table(nsw74psid3$trt, nsw74psid3$nodeg) 0 1 0 63 65 1 54 131 > chisq.test(table(nsw74psid3$trt,nsw74psid3$nodeg)) Pearson's Chi-squared test with Yates' continuity correction data: table(nsw74psid3$trt, nsw74psid3$nodeg) X-squared = 12.12, df = 1, p-value = 0.0004975 > > > > cleanEx(); ..nameEx <- "nsw74psidA" > > ### * nsw74psidA > > flush(stderr()); flush(stdout()) > > ### Name: nsw74psidA > ### Title: A Subset of the nsw74psid1 Data Set > ### Aliases: nsw74psidA > ### Keywords: datasets > > ### ** Examples > > table(nsw74psidA$trt) 0 1 119 133 > > A1.lm <- lm(re78 ~ trt + (age + educ + re74 + re75) + (black + + hisp + marr + nodeg), data = nsw74psidA) > summary(A1.lm)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) -2093.4414 4253.8344 -0.49213 0.62307 trt 2942.9665 1169.9433 2.51548 0.01254 age 71.2005 77.1393 0.92301 0.35692 educ 245.9885 250.6877 0.98125 0.32745 re74 0.3140 0.4002 0.78455 0.43348 re75 0.8717 0.3826 2.27799 0.02360 black 905.5948 1135.1582 0.79777 0.42579 hisp 3723.0383 2211.8204 1.68325 0.09362 marr 24.7218 1104.6055 0.02238 0.98216 nodeg -883.4520 1294.2084 -0.68262 0.49550 > > discA.glm <- glm(formula = trt ~ age + educ + black + hisp + + marr + nodeg + re74 + re75, family = binomial, data = nsw74psidA) > A.scores <- predict(discA.glm) > > options(digits=4) > overlap <- A.scores > -3.5 & A.scores < 3.8 > A.lm <- lm(re78 ~ trt + A.scores, data=nsw74psidA, subset = overlap) > summary(A.lm)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 3588.0 739.8 4.850 2.286e-06 trt 2890.2 1204.6 2.399 1.723e-02 A.scores -333.4 270.7 -1.232 2.193e-01 > > > > > cleanEx(); ..nameEx <- "onesamp" > > ### * onesamp > > flush(stderr()); flush(stdout()) > > ### Name: onesamp > ### Title: Paired Sample t-test > ### Aliases: onesamp > ### Keywords: models > > ### ** Examples > > onesamp(dset = pair65, x = "ambient", y = "heated", xlab = + "Amount of stretch (ambient)", ylab = + "Amount of stretch (heated)") heated 14.60 16.28 6.103 One Sample t-test data: d t = 3.113, df = 8, p-value = 0.01438 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 1.642 11.025 sample estimates: mean of x 6.333 > > > > cleanEx(); ..nameEx <- "onet.permutation" > > ### * onet.permutation > > flush(stderr()); flush(stdout()) > > ### Name: onet.permutation > ### Title: One Sample Permutation t-test > ### Aliases: onet.permutation > ### Keywords: models > > ### ** Examples > > onet.permutation() [1] 0.014 > > > > cleanEx(); ..nameEx <- "oneway.plot" > > ### * oneway.plot > > flush(stderr()); flush(stdout()) > > ### Name: oneway.plot > ### Title: Display of One Way Analysis Results > ### Aliases: oneway.plot > ### Keywords: models > > ### ** Examples > > rice.aov <- aov(ShootDryMass ~ trt, data=rice) > oneway.plot(obj=rice.aov) [1] 1 1 1 1 > > > > cleanEx(); ..nameEx <- "orings" > > ### * orings > > flush(stderr()); flush(stdout()) > > ### Name: orings > ### Title: Challenger O-rings Data > ### Aliases: orings > ### Keywords: datasets > > ### ** Examples > > oldpar <- par(mfrow=c(1,2)) > plot(Total~Temperature, data = orings[c(1,2,4,11,13,18),]) # the > # observations included in the pre-launch charts > plot(Total~Temperature, data = orings) > par(oldpar) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "overlap.density" > > ### * overlap.density > > flush(stderr()); flush(stdout()) > > ### Name: overlap.density > ### Title: Overlapping Density Plots > ### Aliases: overlap.density > ### Keywords: models > > ### ** Examples > > attach(two65) > overlap.density(ambient,heated) [1] 228.9 269.9 > t.test(ambient,heated) Welch Two Sample t-test data: ambient and heated t = -1.990, df = 18.92, p-value = 0.0613 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -19.3109 0.4927 sample estimates: mean of x mean of y 244.1 253.5 > > > > cleanEx(); ..nameEx <- "ozone" > > ### * ozone > > flush(stderr()); flush(stdout()) > > ### Name: ozone > ### Title: Ozone Data > ### Aliases: ozone > ### Keywords: datasets > > ### ** Examples > > AnnualOzone <- ts(ozone$Annual, start=1956) > plot(AnnualOzone) > > > > cleanEx(); ..nameEx <- "pair65" > > ### * pair65 > > flush(stderr()); flush(stdout()) > > ### Name: pair65 > ### Title: Heated Elastic Bands > ### Aliases: pair65 > ### Keywords: datasets > > ### ** Examples > > mean(pair65$heated - pair65$ambient) [1] 6.333 > sd(pair65$heated - pair65$ambient) [1] 6.103 > > > > cleanEx(); ..nameEx <- "panel.corr" > > ### * panel.corr > > flush(stderr()); flush(stdout()) > > ### Name: panel.corr > ### Title: Scatterplot Panel > ### Aliases: panel.corr > ### Keywords: models > > ### ** Examples > > > # correlation between body and brain weights for 20 mice: > > weights <- litters[,-1] > names(weights) <- c("x","y") > weights <- list(weights) > weights[[1]]$xlim <- range(litters[,2]) > weights[[1]]$ylim <- range(litters[,3]) > panelplot(weights, panel.corr, totrows=1, totcols=1) > > > > cleanEx(); ..nameEx <- "panelplot" > > ### * panelplot > > flush(stderr()); flush(stdout()) > > ### Name: panelplot > ### Title: Panel Plot > ### Aliases: panelplot > ### Keywords: models > > ### ** Examples > > x1 <- x2 <- x3 <- (11:30)/5 > y1 <- x1 + rnorm(20)/2 > y2 <- 2 - 0.05 * x1 + 0.1 * ((x1 - 1.75))^4 + 1.25 * rnorm(20) > r <- round(cor(x1, y2), 3) > rho <- round(cor(rank(x1), rank(y2)), 3) > y3 <- (x1 - 3.85)^2 + 0.015 + rnorm(20)/4 > theta <- ((2 * pi) * (1:20))/20 > x4 <- 10 + 4 * cos(theta) > y4 <- 10 + 4 * sin(theta) + (0.5 * rnorm(20)) > r1 <- cor(x1, y1) > xy <- data.frame(x = c(rep(x1, 3), x4), y = c(y1, y2, y3, y4), gp = rep( 1:4, rep(20, 4))) > xy <- split(xy,xy$gp) > xlimdf<-lapply(list(x1,x2,x3,x4),range) > ylimdf<-lapply(list(y1,y2,y3,y4),range) > xy <- lapply(1:4,function(i,u,v,w){u[[i]]$xlim<-v[[i]]; + u[[i]]$ylim<-w[[i]]; u[[i]]},u=xy,v=xlimdf,w=ylimdf) > panelplot(xy,panel=panel.corr,totrows=2,totcols=2, oma=rep(1,4)) > > > > cleanEx(); ..nameEx <- "possum" > > ### * possum > > flush(stderr()); flush(stdout()) > > ### Name: possum > ### Title: Possum Measurements > ### Aliases: possum > ### Keywords: datasets > > ### ** Examples > > boxplot(earconch~sex, data=possum) > pause() Pause. Press to continue... > sex <- as.integer(possum$sex) > pairs(possum[, c(9:11)], oma=c(2,4,5,4), pch=c(0,2:7), col=c("red","blue"), + labels=c("tail\nlength","foot\nlength","ear conch\nlength")) > chh <- par()$cxy[2]; xleg <- 0.05; yleg <- 1.04 > oldpar <- par(xpd=TRUE) > legend(xleg, yleg, c("Cambarville", "Bellbird", "Whian Whian ", + "Byrangery", "Conondale ","Allyn River", "Bulburin"), pch=c(0,2:7), + x.intersp=1, y.intersp=0.75, cex=0.8, xjust=0, bty="n", ncol=4) > text(x=0.2, y=yleg - 2.25*chh, "female", col="red", cex=0.8, bty="n") > text(x=0.75, y=yleg - 2.25*chh, "male", col="blue", cex=0.8, bty="n") > par(oldpar) > pause() Pause. Press to continue... > sapply(possum[,6:14], function(x)max(x,na.rm=TRUE)/min(x,na.rm=TRUE)) hdlngth skullw totlngth taill footlgth earconch eye chest 1.250 1.372 1.287 1.344 1.292 1.395 1.391 1.455 belly 1.600 > pause() Pause. Press to continue... > here <- na.omit(possum$footlgth) > possum.prc <- princomp(possum[here, 6:14]) > pause() Pause. Press to continue... > plot(possum.prc$scores[,1] ~ possum.prc$scores[,2], + col=c("red","blue")[as.numeric(possum$sex[here])], + pch=c(0,2:7)[possum$site[here]], xlab = "PC1", ylab = "PC2") > # NB: We have abbreviated the axis titles > chh <- par()$cxy[2]; xleg <- -15; yleg <- 20.5 > oldpar <- par(xpd=TRUE) > legend(xleg, yleg, c("Cambarville", "Bellbird", "Whian Whian ", + "Byrangery", "Conondale ","Allyn River", "Bulburin"), pch=c(0,2:7), + x.intersp=1, y.intersp=0.75, cex=0.8, xjust=0, bty="n", ncol=4) > text(x=-9, y=yleg - 2.25*chh, "female", col="red", cex=0.8, bty="n") > summary(possum.prc, loadings=TRUE, digits=2) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Standard deviation 5.8341 2.8883 2.25072 1.41491 1.15307 1.03547 0.75705 Proportion of Variance 0.6432 0.1577 0.09574 0.03783 0.02513 0.02026 0.01083 Cumulative Proportion 0.6432 0.8009 0.89664 0.93448 0.95961 0.97987 0.99070 Comp.8 Comp.9 Standard deviation 0.544904 0.441774 Proportion of Variance 0.005611 0.003688 Cumulative Proportion 0.996312 1.000000 Loadings: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 hdlngth 0.659 -0.306 0.479 -0.290 -0.205 0.213 0.104 0.236 skullw 0.322 -0.249 -0.241 -0.297 -0.792 -0.151 0.155 totlngth 0.504 -0.217 -0.455 0.316 0.277 0.154 0.334 -0.425 taill 0.165 -0.223 0.431 0.173 0.270 -0.533 0.584 footlgth 0.124 0.266 -0.546 -0.118 -0.734 0.144 0.106 0.146 earconch -0.362 0.233 0.582 -0.547 -0.232 -0.233 -0.117 -0.230 eye -0.261 0.812 -0.197 -0.420 -0.221 chest 0.197 0.230 -0.262 -0.344 0.269 -0.609 -0.523 belly 0.352 0.725 0.309 0.383 0.227 -0.234 > par(oldpar) > pause() Pause. Press to continue... > require(MASS) Loading required package: MASS Attaching package: 'MASS' The following object(s) are masked from package:DAAG : hills [1] TRUE > here <- !is.na(possum$footlgth) > possum.lda <- lda(site ~ hdlngth+skullw+totlngth+ taill+footlgth+ + earconch+eye+chest+belly, data=possum, subset=here) > options(digits=4) > possum.lda$svd # Examine the singular values [1] 15.7578 3.9372 3.1860 1.5078 1.1420 0.7772 > plot(possum.lda, dimen=3) > # Scatterplot matrix - scores on 1st 3 canonical variates (Figure 11.4) > possum.lda Call: lda(site ~ hdlngth + skullw + totlngth + taill + footlgth + earconch + eye + chest + belly, data = possum, subset = here) Prior probabilities of groups: 1 2 3 4 5 6 7 0.32039 0.11650 0.06796 0.06796 0.12621 0.12621 0.17476 Group means: hdlngth skullw totlngth taill footlgth earconch eye chest belly 1 93.72 57.20 89.71 36.38 73.00 52.58 15.02 27.88 33.27 2 89.85 55.13 81.67 34.67 70.75 52.11 14.37 26.29 31.17 3 94.57 58.90 88.07 37.21 66.60 45.26 16.07 27.57 34.86 4 97.61 61.69 92.24 39.71 68.93 45.84 15.47 29.64 34.57 5 92.18 56.23 86.92 37.65 64.73 43.87 15.38 26.65 32.04 6 89.25 54.19 84.54 37.65 63.07 43.97 15.34 25.23 31.50 7 92.63 57.23 85.69 37.69 65.74 45.86 14.47 26.14 31.92 Coefficients of linear discriminants: LD1 LD2 LD3 LD4 LD5 LD6 hdlngth -0.150534 0.05833 -0.25569 -0.01238 -0.08195 -0.18713 skullw -0.026530 0.03985 -0.24977 0.12448 -0.13648 0.14378 totlngth 0.106957 0.27924 0.30519 -0.18493 -0.13904 -0.08918 taill -0.450631 -0.08960 -0.44579 -0.17304 0.32416 0.49507 footlgth 0.301903 -0.03601 -0.03909 0.07558 0.11910 -0.12606 earconch 0.586270 -0.04357 -0.07310 -0.08819 -0.06376 0.28015 eye -0.056141 0.08921 0.78453 0.46439 0.28481 0.29723 chest 0.090620 0.10418 0.04985 0.12749 0.64746 -0.07871 belly 0.009966 -0.05166 0.09358 0.16715 -0.29033 0.19391 Proportion of trace: LD1 LD2 LD3 LD4 LD5 LD6 0.8927 0.0557 0.0365 0.0082 0.0047 0.0022 > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "possumsites" > > ### * possumsites > > flush(stderr()); flush(stdout()) > > ### Name: possumsites > ### Title: Possum Sites > ### Aliases: possumsites > ### Keywords: datasets > > ### ** Examples > > require(oz) Loading required package: oz [1] TRUE > oz(sections=c(3:5, 11:16)) > attach(possumsites) > points(latitude, longitude, pch=16, col=2) > chw <- par()$cxy[1] > chh <- par()$cxy[2] > posval <- c(2,4,2,2,4,2,2) > text(latitude+(3-posval)*chw/4, longitude, row.names(possumsites), pos=posval) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "powerplot" > > ### * powerplot > > flush(stderr()); flush(stdout()) > > ### Name: powerplot > ### Title: Plot of Power Functions > ### Aliases: powerplot > ### Keywords: models > > ### ** Examples > > oldpar <- par(mfrow = c(2, 3), mar = par()$mar - c( + 1, 1, 1.0, 1), mgp = c(1.5, 0.5, 0), oma=c(0,1,0,1)) > # on.exit(par(oldpar)) > powerplot(expr="sqrt(x)", xlab="") > powerplot(expr="x^0.25", xlab="", ylab="") > powerplot(expr="log(x)", xlab="", ylab="") > powerplot(expr="x^2") > powerplot(expr="x^4", ylab="") > powerplot(expr="exp(x)", ylab="") > par(oldpar) > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "primates" > > ### * primates > > flush(stderr()); flush(stdout()) > > ### Name: primates > ### Title: Primate Body and Brain Weights > ### Aliases: primates > ### Keywords: datasets > > ### ** Examples > > attach(primates) > plot(x=Bodywt, y=Brainwt, pch=16, + xlab="Body weight (kg)", ylab="Brain weight (g)", + xlim=c(5,300), ylim=c(0,1500)) > chw <- par()$cxy[1] > chh <- par()$cxy[2] > text(x=Bodywt+chw, y=Brainwt+c(-.1,0,0,.1,0)*chh, + labels=row.names(primates), adj=0) > detach(primates) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "qreference" > > ### * qreference > > flush(stderr()); flush(stdout()) > > ### Name: qreference > ### Title: Normal QQ Reference Plot > ### Aliases: qreference > ### Keywords: models > > ### ** Examples > > qreference(rt(180,1)) > > qreference(rt(180,1), distribution=function(x) qt(x, df=1)) > > qreference(rexp(180), nrep = 4) > > toycars.lm <- lm(distance ~ angle + factor(car), data = toycars) > qreference(residuals(toycars.lm), nrep = 9) > > > > cleanEx(); ..nameEx <- "races2000" > > ### * races2000 > > flush(stderr()); flush(stdout()) > > ### Name: races2000 > ### Title: Scottish Hill Races Data - 2000 > ### Aliases: races2000 > ### Keywords: datasets > > ### ** Examples > > pairs(races2000[,-11]) > > > > cleanEx(); ..nameEx <- "rainforest" > > ### * rainforest > > flush(stderr()); flush(stdout()) > > ### Name: rainforest > ### Title: Rainforest Data > ### Aliases: rainforest > ### Keywords: datasets > > ### ** Examples > > table(rainforest$species) Acacia mabellae C. fraseri Acmena smithii B. myrtifolia 16 12 26 11 > > > > cleanEx(); ..nameEx <- "rareplants" > > ### * rareplants > > flush(stderr()); flush(stdout()) > > ### Name: rareplants > ### Title: Rare and Endangered Plant Species > ### Aliases: rareplants > ### Keywords: datasets datasets > > ### ** Examples > > chisq.test(rareplants) Pearson's Chi-squared test data: rareplants X-squared = 34.99, df = 6, p-value = 4.336e-06 > > > > cleanEx(); ..nameEx <- "rice" > > ### * rice > > flush(stderr()); flush(stdout()) > > ### Name: rice > ### Title: Genetically Modified and Wild Type Rice Data > ### Aliases: rice > ### Keywords: datasets > > ### ** Examples > > print("One and Two-Way Comparisons - Example 4.5") [1] "One and Two-Way Comparisons - Example 4.5" > attach(rice) > oldpar <- par(las = 2) > stripchart(ShootDryMass ~ trt, pch=1, cex=1, xlab="Level of factor 1") > detach(rice) > pause() Pause. Press to continue... > rice.aov <- aov(ShootDryMass ~ trt, data=rice); anova(rice.aov) Analysis of Variance Table Response: ShootDryMass Df Sum Sq Mean Sq F value Pr(>F) trt 5 68326 13665 36.7 <2e-16 *** Residuals 66 24562 372 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(rice.aov) Analysis of Variance Table Response: ShootDryMass Df Sum Sq Mean Sq F value Pr(>F) trt 5 68326 13665 36.7 <2e-16 *** Residuals 66 24562 372 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > pause() Pause. Press to continue... > summary.lm(rice.aov)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 108.33 5.569 19.453 4.919e-29 trtNH4Cl -58.08 7.876 -7.375 3.475e-10 trtNH4NO3 -35.00 7.876 -4.444 3.454e-05 trtF10 +ANU843 -101.00 7.876 -12.824 1.381e-19 trtNH4Cl +ANU843 -61.75 7.876 -7.841 5.107e-11 trtNH4NO3 +ANU843 -36.83 7.876 -4.677 1.488e-05 > pause() Pause. Press to continue... > rice$trt <- relevel(rice$trt, ref="NH4Cl") > # Set NH4Cl as the baseline > > fac1 <- factor(sapply(strsplit(as.character(rice$trt)," \\+"), function(x)x[1])) > anu843 <- sapply(strsplit(as.character(rice$trt), "\\+"), + function(x)c("wt","ANU843")[length(x)]) > anu843 <- factor(anu843, levels=c("wt", "ANU843")) > attach(rice) > interaction.plot(fac1, anu843, ShootDryMass) > detach(rice) > par(oldpar) > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "roller" > > ### * roller > > flush(stderr()); flush(stdout()) > > ### Name: roller > ### Title: Lawn Roller Data > ### Aliases: roller > ### Keywords: datasets > > ### ** Examples > > plot(roller) > roller.lm <- lm(depression ~ weight, data = roller) > plot(roller.lm, which = 4) > > > > cleanEx(); ..nameEx <- "science" > > ### * science > > flush(stderr()); flush(stdout()) > > ### Name: science > ### Title: School Science Survey Data > ### Aliases: science > ### Keywords: datasets > > ### ** Examples > > attach(science) > classmeans <- aggregate(like, by=list(PrivPub, Class), mean) > names(classmeans) <- c("PrivPub","Class","like") > dim(classmeans) [1] 66 3 > > attach(classmeans) > boxplot(split(like, PrivPub), ylab = "Class average of attitude to science score", boxwex = 0.4) > rug(like[PrivPub == "private"], side = 2) > rug(like[PrivPub == "public"], side = 4) > detach(classmeans) > > require(nlme) Loading required package: nlme [1] TRUE > science.lme <- lme(fixed = like ~ sex + PrivPub, + data = science, random = ~ 1 | school/Class, na.action=na.omit) > summary(science.lme)$tTable # Print coefficients. Value Std.Error DF t-value p-value (Intercept) 4.7220 0.1625 1316 29.066 6.785e-144 sexm 0.1824 0.0982 1316 1.857 6.347e-02 PrivPubpublic 0.4118 0.1857 39 2.217 3.253e-02 > > science1.lme <- lme(fixed = like ~ sex + PrivPub, data = science, + random = ~ 1 | Class, na.action=na.exclude) > summary(science1.lme)$tTable # Table of coefficients Value Std.Error DF t-value p-value (Intercept) 4.7218 0.1624 1316 29.072 6.202e-144 sexm 0.1823 0.0982 1316 1.857 6.356e-02 PrivPubpublic 0.4117 0.1857 64 2.217 3.018e-02 > > intervals(science1.lme, which="var-cov")[[1]]$Class^2 lower est. upper sd((Intercept)) 0.1904 0.3206 0.5398 > intervals(science1.lme, which="var-cov")[[2]]^2 lower est. upper 2.828 3.052 3.294 attr(,"label") [1] "Within-group standard error:" > > science.lme <- lme(fixed = like ~ sex + PrivPub, + data = science, random = ~ 1 | Class/school, na.action=na.exclude) > res <- residuals(science.lme) > hat <- fitted(science.lme) > coplot(res ~ hat|science$PrivPub[is.na(science$sex)!=TRUE], + xlab="Fitted values", ylab="Residuals") > detach(science) > > > > cleanEx(); ..nameEx <- "seedrates" > > ### * seedrates > > flush(stderr()); flush(stdout()) > > ### Name: seedrates > ### Title: Barley Seeding Rate Data > ### Aliases: seedrates > ### Keywords: datasets > > ### ** Examples > > plot(grain~rate,data=seedrates,xlim=c(50,180),ylim=c(15.5,22),axes=FALSE) > new.df<-data.frame(rate=(2:8)*25) > seedrates.lm1<-lm(grain~rate,data=seedrates) > seedrates.lm2<-lm(grain~rate+I(rate^2),data=seedrates) > hat1<-predict(seedrates.lm1,newdata=new.df,interval="confidence") > hat2<-predict(seedrates.lm2,newdata=new.df,interval="confidence") > axis(1,at=new.df$rate); axis(2); box() > z1<-spline(new.df$rate, hat1[,"fit"]); z2<-spline(new.df$rate, + hat2[,"fit"]) > rate<-new.df$rate; lines(z1$x,z1$y) > lines(spline(rate,hat1[,"lwr"]),lty=1,col=3) > lines(spline(rate,hat1[,"upr"]),lty=1,col=3) > lines(z2$x,z2$y,lty=4) > lines(spline(rate,hat2[,"lwr"]),lty=4,col=3) > lines(spline(rate,hat2[,"upr"]),lty=4,col=3) > > > cleanEx(); ..nameEx <- "show.colors" > > ### * show.colors > > flush(stderr()); flush(stdout()) > > ### Name: show.colors > ### Title: Show R's Colors > ### Aliases: show.colors > ### Keywords: models > > ### ** Examples > > require(MASS) Loading required package: MASS Attaching package: 'MASS' The following object(s) are masked from package:DAAG : hills [1] TRUE > show.colors() Loading required package: mva Warning: package 'mva' has been merged into 'stats' [1] 49 390 > > > > cleanEx(); ..nameEx <- "simulate.linear" > > ### * simulate.linear > > flush(stderr()); flush(stdout()) > > ### Name: simulate.linear > ### Title: Simulation of Linear Models for ANOVA vs. Regression Comparison > ### Aliases: simulate.linear > ### Keywords: models > > ### ** Examples > > simulate.linear() Proportion of datasets where linear p-value < aov p-value = 0.87 > > > > cleanEx(); ..nameEx <- "socsupport" > > ### * socsupport > > flush(stderr()); flush(stdout()) > > ### Name: socsupport > ### Title: Social Support Data > ### Aliases: socsupport > ### Keywords: datasets > > ### ** Examples > > attach(socsupport) > > not.na <- apply(socsupport[,9:19], 1, function(x)!any(is.na(x))) > ss.pr1 <- princomp(as.matrix(socsupport[not.na, 9:19]), cor=TRUE) > pairs(ss.pr1$scores[,1:3]) > sort(-ss.pr1$scores[,1]) # Minus the largest value appears first 36 30 62 66 10 32 81 73 -9.89867 -5.61459 -4.70607 -4.52908 -4.20789 -3.90331 -3.90226 -3.70374 2 75 78 5 59 12 14 63 -3.49741 -3.29753 -3.02037 -2.94128 -2.82118 -2.71366 -2.66089 -2.44295 85 83 44 68 31 89 42 22 -2.41774 -2.07015 -1.76393 -1.73443 -1.70018 -1.65712 -1.54761 -1.45063 18 92 49 93 79 74 21 38 -1.44271 -1.21288 -1.18136 -1.01838 -1.00279 -0.93347 -0.70089 -0.65869 11 45 52 69 80 33 15 8 -0.58748 -0.50429 -0.33703 -0.32970 -0.11435 -0.07399 -0.04064 -0.03503 6 17 60 39 41 50 90 72 -0.03264 0.16242 0.39149 0.51309 0.51694 0.56755 0.56803 0.61637 40 88 94 57 95 77 28 56 0.63450 0.69124 0.69680 0.70580 0.75363 0.78963 0.81077 0.93662 46 24 86 26 27 70 4 43 0.96231 1.10715 1.11452 1.12373 1.12384 1.12652 1.14337 1.21303 3 20 54 65 1 84 29 76 1.26994 1.27087 1.27507 1.56123 1.70041 1.86789 2.15292 2.21968 64 7 16 91 71 82 37 35 2.50528 2.67983 2.77527 2.79633 2.92374 2.95060 2.98158 2.99505 87 23 61 53 34 13 58 48 3.01844 3.02578 3.16507 3.25874 3.33247 3.41662 3.53631 3.60287 47 25 3.79593 4.06173 > pause() Pause. Press to continue... > not.na[36] <- FALSE > ss.pr <- princomp(as.matrix(socsupport[not.na, 9:19]), cor=TRUE) > summary(ss.pr) # Examine the contribution of the components Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Standard deviation 2.3937 1.2190 1.1366 0.84477 0.75447 0.69536 0.49726 Proportion of Variance 0.5209 0.1351 0.1174 0.06488 0.05175 0.04396 0.02248 Cumulative Proportion 0.5209 0.6560 0.7734 0.83833 0.89007 0.93403 0.95651 Comp.8 Comp.9 Comp.10 Comp.11 Standard deviation 0.45610 0.35952 0.295553 0.231892 Proportion of Variance 0.01891 0.01175 0.007941 0.004889 Cumulative Proportion 0.97542 0.98717 0.995111 1.000000 > pause() Pause. Press to continue... > # We now regress BDI on the first six principal components: > ss.lm <- lm(BDI[not.na] ~ ss.pr$scores[, 1:6], data=socsupport) > summary(ss.lm)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 10.4607 0.8934 11.7088 3.489e-19 ss.pr$scores[, 1:6]Comp.1 1.3113 0.3732 3.5134 7.229e-04 ss.pr$scores[, 1:6]Comp.2 -0.3959 0.7329 -0.5402 5.905e-01 ss.pr$scores[, 1:6]Comp.3 0.6036 0.7860 0.7679 4.447e-01 ss.pr$scores[, 1:6]Comp.4 1.4248 1.0576 1.3473 1.816e-01 ss.pr$scores[, 1:6]Comp.5 2.1459 1.1841 1.8122 7.362e-02 ss.pr$scores[, 1:6]Comp.6 1.2882 1.2848 1.0027 3.190e-01 > pause() Pause. Press to continue... > ss.pr$loadings[,1] emotional emotionalsat tangible tangiblesat affect affectsat -0.3201 -0.2979 -0.2467 -0.2887 -0.3069 -0.2883 psi psisat esupport psupport supsources -0.3628 -0.3322 -0.2886 -0.2846 -0.2846 > plot(BDI[not.na] ~ ss.pr$scores[ ,1], col=as.numeric(gender), + pch=as.numeric(gender), xlab ="1st principal component", ylab="BDI") > topleft <- par()$usr[c(1,4)] > legend(topleft[1], topleft[2], col=1:2, pch=1:2, legend=levels(gender)) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "softbacks" > > ### * softbacks > > flush(stderr()); flush(stdout()) > > ### Name: softbacks > ### Title: Measurements on a Selection of Paperback Books > ### Aliases: softbacks > ### Keywords: datasets > > ### ** Examples > > print("Outliers in Simple Regression - Example 5.2") [1] "Outliers in Simple Regression - Example 5.2" > paperback.lm <- lm(weight ~ volume, data=softbacks) > summary(paperback.lm) Call: lm(formula = weight ~ volume, data = softbacks) Residuals: Min 1Q Median 3Q Max -89.67 -39.89 -25.00 9.07 215.91 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 41.372 97.559 0.42 0.68629 volume 0.686 0.106 6.47 0.00064 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 102 on 6 degrees of freedom Multiple R-Squared: 0.875, Adjusted R-squared: 0.854 F-statistic: 41.9 on 1 and 6 DF, p-value: 0.000644 > plot(paperback.lm) > > > > cleanEx(); ..nameEx <- "spam7" > > ### * spam7 > > flush(stderr()); flush(stdout()) > > ### Name: spam7 > ### Title: Spam E-mail Data > ### Aliases: spam7 > ### Keywords: datasets > > ### ** Examples > > require(rpart) Loading required package: rpart [1] TRUE > spam.rpart <- rpart(formula = yesno ~ crl.tot + dollar + bang + + money + n000 + make, data=spam7) > plot(spam.rpart) > text(spam.rpart) > > > > > cleanEx(); ..nameEx <- "sugar" > > ### * sugar > > flush(stderr()); flush(stdout()) > > ### Name: sugar > ### Title: Sugar Data > ### Aliases: sugar > ### Keywords: datasets > > ### ** Examples > > sugar.aov <- aov(weight ~ trt, data=sugar) > fitted.values(sugar.aov) 1 2 3 4 5 6 7 8 9 10 11 12 83.23 83.23 83.23 61.83 61.83 61.83 67.50 67.50 67.50 48.90 48.90 48.90 > summary.lm(sugar.aov) Call: aov(formula = weight ~ trt, data = sugar) Residuals: Min 1Q Median 3Q Max -13.333 -2.783 -0.617 2.175 14.567 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 83.23 4.47 18.61 7.2e-08 *** trtA -21.40 6.33 -3.38 0.00960 ** trtB -15.73 6.33 -2.49 0.03768 * trtC -34.33 6.33 -5.43 0.00062 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 7.75 on 8 degrees of freedom Multiple R-Squared: 0.791, Adjusted R-squared: 0.713 F-statistic: 10.1 on 3 and 8 DF, p-value: 0.00425 > sugar.aov <- aov(formula = weight ~ trt, data = sugar) > summary.lm(sugar.aov) Call: aov(formula = weight ~ trt, data = sugar) Residuals: Min 1Q Median 3Q Max -13.333 -2.783 -0.617 2.175 14.567 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 83.23 4.47 18.61 7.2e-08 *** trtA -21.40 6.33 -3.38 0.00960 ** trtB -15.73 6.33 -2.49 0.03768 * trtC -34.33 6.33 -5.43 0.00062 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 7.75 on 8 degrees of freedom Multiple R-Squared: 0.791, Adjusted R-squared: 0.713 F-statistic: 10.1 on 3 and 8 DF, p-value: 0.00425 > > > > cleanEx(); ..nameEx <- "tinting" > > ### * tinting > > flush(stderr()); flush(stdout()) > > ### Name: tinting > ### Title: Car Window Tinting Experiment Data > ### Aliases: tinting > ### Keywords: datasets > > ### ** Examples > > require(lattice) Loading required package: lattice [1] TRUE > levels(tinting$agegp) <- toupper.initial(levels(tinting$agegp)) > xyplot(csoa ~ it | sex * agegp, data=tinting) # Simple use of xyplot() > pause() Pause. Press to continue... > xyplot(csoa ~ it|sex*agegp, data=tinting, panel=panel.superpose, groups=target) > pause() Pause. Press to continue... > xyplot(csoa ~ it|sex*agegp, data=tinting, panel=panel.superpose, col=1:2, + groups=target, key=list(x=0.14, y=0.84, points=list(pch=rep(1,2), + col=1:2), text=list(levels(tinting$target), col=1:2), border=TRUE)) > pause() Pause. Press to continue... > xyplot(csoa ~ it|sex*agegp, data=tinting, panel=panel.superpose, + groups=tint, type=c("p","smooth"), span=0.8, col=1:3, + key=list(x=0.14, y=0.84, points=list(pch=rep(1,2), col=1:3), + text=list(levels(tinting$tint), col=1:3), border=TRUE)) > > > > cleanEx(); ..nameEx <- "toupper.initial" > > ### * toupper.initial > > flush(stderr()); flush(stdout()) > > ### Name: toupper.initial > ### Title: Converts initial character of a string to upper case > ### Aliases: toupper.initial > ### Keywords: models > > ### ** Examples > > toupper.initial(names(tinting)[c(3,4)]) [1] "Age" "Sex" > > library(lattice) > levels(tinting$agegp) <- toupper.initial(levels(tinting$agegp)) > xyplot(csoa ~ it | sex * agegp, data=tinting) > > > > cleanEx(); ..nameEx <- "toycars" > > ### * toycars > > flush(stderr()); flush(stdout()) > > ### Name: toycars > ### Title: Toy Cars Data > ### Aliases: toycars > ### Keywords: datasets > > ### ** Examples > > toycars.lm <- lm(distance ~ angle + factor(car), data=toycars) > summary(toycars.lm) Call: lm(formula = distance ~ angle + factor(car), data = toycars) Residuals: Min 1Q Median 3Q Max -0.09811 -0.04240 -0.00669 0.01741 0.17251 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.09252 0.03467 2.67 0.0137 * angle 0.18854 0.00995 18.96 1.5e-15 *** factor(car)2 0.11111 0.03195 3.48 0.0020 ** factor(car)3 -0.08222 0.03195 -2.57 0.0170 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0678 on 23 degrees of freedom Multiple R-Squared: 0.945, Adjusted R-squared: 0.938 F-statistic: 132 on 3 and 23 DF, p-value: 1.22e-14 > > > > cleanEx(); ..nameEx <- "two65" > > ### * two65 > > flush(stderr()); flush(stdout()) > > ### Name: two65 > ### Title: Unpaired Heated Elastic Bands > ### Aliases: two65 > ### Keywords: datasets > > ### ** Examples > > twot.permutation(two65$ambient,two65$heated) # two sample permutation test [1] 0.0545 > > > > cleanEx(); ..nameEx <- "twot.permutation" > > ### * twot.permutation > > flush(stderr()); flush(stdout()) > > ### Name: twot.permutation > ### Title: Two Sample Permutation Test > ### Aliases: twot.permutation > ### Keywords: models > > ### ** Examples > > twot.permutation() [1] 0.0545 > > > > cleanEx(); ..nameEx <- "vif" > > ### * vif > > flush(stderr()); flush(stdout()) > > ### Name: vif > ### Title: Variance Inflation Factors > ### Aliases: vif > ### Keywords: models > > ### ** Examples > > litters.lm <- lm(brainwt ~ bodywt + lsize, data = litters) > vif(litters.lm) bodywt lsize 11.33 11.33 > > carprice1.lm <- lm(gpm100 ~ Type+Min.Price+Price+Max.Price+Range.Price, + data=carprice) > vif(carprice1.lm) TypeLarge TypeMidsize TypeSmall TypeSporty TypeVan Min.Price 2.722e+00 2.332e+00 1.804e+00 2.174e+00 1.778e+00 2.992e+04 Price Max.Price 1.206e+05 3.340e+04 > > carprice.lm <- lm(gpm100 ~ Type + Price, data = carprice) > vif(carprice1.lm) TypeLarge TypeMidsize TypeSmall TypeSporty TypeVan Min.Price 2.722e+00 2.332e+00 1.804e+00 2.174e+00 1.778e+00 2.992e+04 Price Max.Price 1.206e+05 3.340e+04 > > > > > cleanEx(); ..nameEx <- "wages1833" > > ### * wages1833 > > flush(stderr()); flush(stdout()) > > ### Name: wages1833 > ### Title: Wages of Lancashire Cotton Factory Workers in 1833 > ### Aliases: wages1833 > ### Keywords: datasets > > ### ** Examples > > attach(wages1833) > plot(mwage~age,ylim=range(c(mwage,fwage[fwage>0]))) > points(fwage[fwage>0]~age[fwage>0],pch=15,col="red") > lines(lowess(age,mwage)) > lines(lowess(age[fwage>0],fwage[fwage>0]),col="red") > > > > ### *