R Under development (unstable) (2019-03-04 r76193) -- "Unsuffered Consequences" Copyright (C) 2019 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) 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. Natural language support but running in an English locale 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 an HTML browser interface to help. Type 'q()' to quit R. > pkgname <- "datasets" > source(file.path(R.home("share"), "R", "examples-header.R")) > options(warn = 1) > library('datasets') > > base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') > base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv') > cleanEx() > nameEx("AirPassengers") > ### * AirPassengers > > flush(stderr()); flush(stdout()) > > ### Name: AirPassengers > ### Title: Monthly Airline Passenger Numbers 1949-1960 > ### Aliases: AirPassengers > ### Keywords: datasets > > ### ** Examples > > ## Not run: > ##D ## These are quite slow and so not run by example(AirPassengers) > ##D > ##D ## The classic 'airline model', by full ML > ##D (fit <- arima(log10(AirPassengers), c(0, 1, 1), > ##D seasonal = list(order = c(0, 1, 1), period = 12))) > ##D update(fit, method = "CSS") > ##D update(fit, x = window(log10(AirPassengers), start = 1954)) > ##D pred <- predict(fit, n.ahead = 24) > ##D tl <- pred$pred - 1.96 * pred$se > ##D tu <- pred$pred + 1.96 * pred$se > ##D ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty = c(1, 2, 2)) > ##D > ##D ## full ML fit is the same if the series is reversed, CSS fit is not > ##D ap0 <- rev(log10(AirPassengers)) > ##D attributes(ap0) <- attributes(AirPassengers) > ##D arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12)) > ##D arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12), > ##D method = "CSS") > ##D > ##D ## Structural Time Series > ##D ap <- log10(AirPassengers) - 2 > ##D (fit <- StructTS(ap, type = "BSM")) > ##D par(mfrow = c(1, 2)) > ##D plot(cbind(ap, fitted(fit)), plot.type = "single") > ##D plot(cbind(ap, tsSmooth(fit)), plot.type = "single") > ## End(Not run) > > > cleanEx() > nameEx("BOD") > ### * BOD > > flush(stderr()); flush(stdout()) > > ### Name: BOD > ### Title: Biochemical Oxygen Demand > ### Aliases: BOD > ### Keywords: datasets > > ### ** Examples > > ## Don't show: > options(show.nls.convergence=FALSE) > old <- options(digits = 5) > ## End(Don't show) > require(stats) > # simplest form of fitting a first-order model to these data > fm1 <- nls(demand ~ A*(1-exp(-exp(lrc)*Time)), data = BOD, + start = c(A = 20, lrc = log(.35))) > coef(fm1) A lrc 19.14258 -0.63282 > fm1 Nonlinear regression model model: demand ~ A * (1 - exp(-exp(lrc) * Time)) data: BOD A lrc 19.143 -0.633 residual sum-of-squares: 26 > # using the plinear algorithm > fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD, + start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE) 32.946 : -1.0498 22.1260 25.992 : -0.62572 19.10319 25.99 : -0.6327 19.1419 25.99 : -0.63282 19.14256 > # using a self-starting model > fm3 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD) > summary(fm3) Formula: demand ~ SSasympOrig(Time, A, lrc) Parameters: Estimate Std. Error t value Pr(>|t|) A 19.143 2.496 7.67 0.0016 ** lrc -0.633 0.382 -1.65 0.1733 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.55 on 4 degrees of freedom > ## Don't show: > options(old) > ## End(Don't show) > > > > cleanEx() > nameEx("ChickWeight") > ### * ChickWeight > > flush(stderr()); flush(stdout()) > > ### Name: ChickWeight > ### Title: Weight versus age of chicks on different diets > ### Aliases: ChickWeight > ### Keywords: datasets > > ### ** Examples > > > cleanEx() > nameEx("DNase") > ### * DNase > > flush(stderr()); flush(stdout()) > > ### Name: DNase > ### Title: Elisa assay of DNase > ### Aliases: DNase > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## Don't show: > options(show.nls.convergence=FALSE) > ## End(Don't show) > coplot(density ~ conc | Run, data = DNase, + show.given = FALSE, type = "b") > coplot(density ~ log(conc) | Run, data = DNase, + show.given = FALSE, type = "b") > ## fit a representative run > fm1 <- nls(density ~ SSlogis( log(conc), Asym, xmid, scal ), + data = DNase, subset = Run == 1) > ## compare with a four-parameter logistic > fm2 <- nls(density ~ SSfpl( log(conc), A, B, xmid, scal ), + data = DNase, subset = Run == 1) > summary(fm2) Formula: density ~ SSfpl(log(conc), A, B, xmid, scal) Parameters: Estimate Std. Error t value Pr(>|t|) A -0.007897 0.017200 -0.459 0.654 B 2.377239 0.109516 21.707 5.35e-11 *** xmid 1.507403 0.102080 14.767 4.65e-09 *** scal 1.062579 0.056996 18.643 3.16e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.01981 on 12 degrees of freedom > anova(fm1, fm2) Analysis of Variance Table Model 1: density ~ SSlogis(log(conc), Asym, xmid, scal) Model 2: density ~ SSfpl(log(conc), A, B, xmid, scal) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 13 0.0047896 2 12 0.0047073 1 8.2314e-05 0.2098 0.6551 > > > > cleanEx() > nameEx("Formaldehyde") > ### * Formaldehyde > > flush(stderr()); flush(stdout()) > > ### Name: Formaldehyde > ### Title: Determination of Formaldehyde > ### Aliases: Formaldehyde > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > plot(optden ~ carb, data = Formaldehyde, + xlab = "Carbohydrate (ml)", ylab = "Optical Density", + main = "Formaldehyde data", col = 4, las = 1) > abline(fm1 <- lm(optden ~ carb, data = Formaldehyde)) > summary(fm1) Call: lm(formula = optden ~ carb, data = Formaldehyde) Residuals: 1 2 3 4 5 6 -0.006714 0.001029 0.002771 0.007143 0.007514 -0.011743 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.005086 0.007834 0.649 0.552 carb 0.876286 0.013535 64.744 3.41e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.008649 on 4 degrees of freedom Multiple R-squared: 0.999, Adjusted R-squared: 0.9988 F-statistic: 4192 on 1 and 4 DF, p-value: 3.409e-07 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0)) > plot(fm1) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("HairEyeColor") > ### * HairEyeColor > > flush(stderr()); flush(stdout()) > > ### Name: HairEyeColor > ### Title: Hair and Eye Color of Statistics Students > ### Aliases: HairEyeColor > ### Keywords: datasets > > ### ** Examples > > require(graphics) > ## Full mosaic > mosaicplot(HairEyeColor) > ## Aggregate over sex (as in Snee's original data) > x <- apply(HairEyeColor, c(1, 2), sum) > x Eye Hair Brown Blue Hazel Green Black 68 20 15 5 Brown 119 84 54 29 Red 26 17 14 14 Blond 7 94 10 16 > mosaicplot(x, main = "Relation between hair and eye color") > > > > cleanEx() > nameEx("Harman23.cor") > ### * Harman23.cor > > flush(stderr()); flush(stdout()) > > ### Name: Harman23.cor > ### Title: Harman Example 2.3 > ### Aliases: Harman23.cor > ### Keywords: datasets > > ### ** Examples > > require(stats) > (Harman23.FA <- factanal(factors = 1, covmat = Harman23.cor)) Call: factanal(factors = 1, covmat = Harman23.cor) Uniquenesses: height arm.span forearm lower.leg weight 0.158 0.135 0.190 0.187 0.760 bitro.diameter chest.girth chest.width 0.829 0.877 0.801 Loadings: Factor1 height 0.918 arm.span 0.930 forearm 0.900 lower.leg 0.902 weight 0.490 bitro.diameter 0.413 chest.girth 0.351 chest.width 0.446 Factor1 SS loadings 4.064 Proportion Var 0.508 Test of the hypothesis that 1 factor is sufficient. The chi square statistic is 611.44 on 20 degrees of freedom. The p-value is 1.12e-116 > for(factors in 2:4) print(update(Harman23.FA, factors = factors)) Call: factanal(factors = factors, covmat = Harman23.cor) Uniquenesses: height arm.span forearm lower.leg weight 0.170 0.107 0.166 0.199 0.089 bitro.diameter chest.girth chest.width 0.364 0.416 0.537 Loadings: Factor1 Factor2 height 0.865 0.287 arm.span 0.927 0.181 forearm 0.895 0.179 lower.leg 0.859 0.252 weight 0.233 0.925 bitro.diameter 0.194 0.774 chest.girth 0.134 0.752 chest.width 0.278 0.621 Factor1 Factor2 SS loadings 3.335 2.617 Proportion Var 0.417 0.327 Cumulative Var 0.417 0.744 Test of the hypothesis that 2 factors are sufficient. The chi square statistic is 75.74 on 13 degrees of freedom. The p-value is 6.94e-11 Call: factanal(factors = factors, covmat = Harman23.cor) Uniquenesses: height arm.span forearm lower.leg weight 0.127 0.005 0.193 0.157 0.090 bitro.diameter chest.girth chest.width 0.359 0.411 0.490 Loadings: Factor1 Factor2 Factor3 height 0.886 0.267 -0.130 arm.span 0.937 0.195 0.280 forearm 0.874 0.188 lower.leg 0.877 0.230 -0.145 weight 0.242 0.916 -0.106 bitro.diameter 0.193 0.777 chest.girth 0.137 0.755 chest.width 0.261 0.646 0.159 Factor1 Factor2 Factor3 SS loadings 3.379 2.628 0.162 Proportion Var 0.422 0.329 0.020 Cumulative Var 0.422 0.751 0.771 Test of the hypothesis that 3 factors are sufficient. The chi square statistic is 22.81 on 7 degrees of freedom. The p-value is 0.00184 Call: factanal(factors = factors, covmat = Harman23.cor) Uniquenesses: height arm.span forearm lower.leg weight 0.137 0.005 0.191 0.116 0.138 bitro.diameter chest.girth chest.width 0.283 0.178 0.488 Loadings: Factor1 Factor2 Factor3 Factor4 height 0.879 0.277 -0.115 arm.span 0.937 0.194 0.277 forearm 0.875 0.191 lower.leg 0.887 0.209 0.135 -0.188 weight 0.246 0.882 0.111 -0.109 bitro.diameter 0.187 0.822 chest.girth 0.117 0.729 0.526 chest.width 0.263 0.644 0.141 Factor1 Factor2 Factor3 Factor4 SS loadings 3.382 2.595 0.323 0.165 Proportion Var 0.423 0.324 0.040 0.021 Cumulative Var 0.423 0.747 0.787 0.808 Test of the hypothesis that 4 factors are sufficient. The chi square statistic is 4.63 on 2 degrees of freedom. The p-value is 0.0988 > > > > cleanEx() > nameEx("Harman74.cor") > ### * Harman74.cor > > flush(stderr()); flush(stdout()) > > ### Name: Harman74.cor > ### Title: Harman Example 7.4 > ### Aliases: Harman74.cor > ### Keywords: datasets > > ### ** Examples > > require(stats) > (Harman74.FA <- factanal(factors = 1, covmat = Harman74.cor)) Call: factanal(factors = 1, covmat = Harman74.cor) Uniquenesses: VisualPerception Cubes PaperFormBoard 0.677 0.866 0.830 Flags GeneralInformation PargraphComprehension 0.768 0.487 0.491 SentenceCompletion WordClassification WordMeaning 0.500 0.514 0.474 Addition Code CountingDots 0.818 0.731 0.824 StraightCurvedCapitals WordRecognition NumberRecognition 0.681 0.833 0.863 FigureRecognition ObjectNumber NumberFigure 0.775 0.812 0.778 FigureWord Deduction NumericalPuzzles 0.816 0.612 0.676 ProblemReasoning SeriesCompletion ArithmeticProblems 0.619 0.524 0.593 Loadings: Factor1 VisualPerception 0.569 Cubes 0.366 PaperFormBoard 0.412 Flags 0.482 GeneralInformation 0.716 PargraphComprehension 0.713 SentenceCompletion 0.707 WordClassification 0.697 WordMeaning 0.725 Addition 0.426 Code 0.519 CountingDots 0.419 StraightCurvedCapitals 0.565 WordRecognition 0.408 NumberRecognition 0.370 FigureRecognition 0.474 ObjectNumber 0.434 NumberFigure 0.471 FigureWord 0.429 Deduction 0.623 NumericalPuzzles 0.569 ProblemReasoning 0.617 SeriesCompletion 0.690 ArithmeticProblems 0.638 Factor1 SS loadings 7.438 Proportion Var 0.310 Test of the hypothesis that 1 factor is sufficient. The chi square statistic is 622.91 on 252 degrees of freedom. The p-value is 2.28e-33 > for(factors in 2:5) print(update(Harman74.FA, factors = factors)) Call: factanal(factors = factors, covmat = Harman74.cor) Uniquenesses: VisualPerception Cubes PaperFormBoard 0.650 0.864 0.844 Flags GeneralInformation PargraphComprehension 0.778 0.375 0.316 SentenceCompletion WordClassification WordMeaning 0.319 0.503 0.258 Addition Code CountingDots 0.670 0.608 0.581 StraightCurvedCapitals WordRecognition NumberRecognition 0.567 0.832 0.850 FigureRecognition ObjectNumber NumberFigure 0.743 0.770 0.625 FigureWord Deduction NumericalPuzzles 0.792 0.629 0.579 ProblemReasoning SeriesCompletion ArithmeticProblems 0.634 0.539 0.553 Loadings: Factor1 Factor2 VisualPerception 0.506 0.306 Cubes 0.304 0.209 PaperFormBoard 0.297 0.260 Flags 0.327 0.339 GeneralInformation 0.240 0.753 PargraphComprehension 0.171 0.809 SentenceCompletion 0.163 0.809 WordClassification 0.344 0.615 WordMeaning 0.148 0.849 Addition 0.563 0.115 Code 0.591 0.207 CountingDots 0.647 StraightCurvedCapitals 0.612 0.241 WordRecognition 0.315 0.263 NumberRecognition 0.328 0.205 FigureRecognition 0.457 0.218 ObjectNumber 0.431 0.209 NumberFigure 0.601 0.116 FigureWord 0.399 0.222 Deduction 0.379 0.477 NumericalPuzzles 0.604 0.237 ProblemReasoning 0.390 0.462 SeriesCompletion 0.486 0.474 ArithmeticProblems 0.544 0.389 Factor1 Factor2 SS loadings 4.573 4.548 Proportion Var 0.191 0.190 Cumulative Var 0.191 0.380 Test of the hypothesis that 2 factors are sufficient. The chi square statistic is 420.24 on 229 degrees of freedom. The p-value is 2.01e-13 Call: factanal(factors = factors, covmat = Harman74.cor) Uniquenesses: VisualPerception Cubes PaperFormBoard 0.500 0.793 0.662 Flags GeneralInformation PargraphComprehension 0.694 0.352 0.316 SentenceCompletion WordClassification WordMeaning 0.300 0.502 0.256 Addition Code CountingDots 0.200 0.586 0.494 StraightCurvedCapitals WordRecognition NumberRecognition 0.569 0.838 0.848 FigureRecognition ObjectNumber NumberFigure 0.643 0.780 0.635 FigureWord Deduction NumericalPuzzles 0.788 0.590 0.580 ProblemReasoning SeriesCompletion ArithmeticProblems 0.597 0.498 0.500 Loadings: Factor1 Factor2 Factor3 VisualPerception 0.176 0.656 0.198 Cubes 0.122 0.428 PaperFormBoard 0.145 0.563 Flags 0.239 0.487 0.107 GeneralInformation 0.745 0.191 0.237 PargraphComprehension 0.780 0.249 0.118 SentenceCompletion 0.802 0.175 0.160 WordClassification 0.571 0.327 0.256 WordMeaning 0.821 0.248 Addition 0.162 -0.118 0.871 Code 0.198 0.219 0.572 CountingDots 0.179 0.688 StraightCurvedCapitals 0.190 0.381 0.499 WordRecognition 0.231 0.253 0.210 NumberRecognition 0.158 0.299 0.195 FigureRecognition 0.108 0.557 0.186 ObjectNumber 0.178 0.267 0.342 NumberFigure 0.427 0.424 FigureWord 0.167 0.355 0.240 Deduction 0.392 0.472 0.181 NumericalPuzzles 0.178 0.406 0.473 ProblemReasoning 0.382 0.473 0.182 SeriesCompletion 0.379 0.528 0.283 ArithmeticProblems 0.377 0.226 0.554 Factor1 Factor2 Factor3 SS loadings 3.802 3.488 3.186 Proportion Var 0.158 0.145 0.133 Cumulative Var 0.158 0.304 0.436 Test of the hypothesis that 3 factors are sufficient. The chi square statistic is 295.59 on 207 degrees of freedom. The p-value is 5.12e-05 Call: factanal(factors = factors, covmat = Harman74.cor) Uniquenesses: VisualPerception Cubes PaperFormBoard 0.438 0.780 0.644 Flags GeneralInformation PargraphComprehension 0.651 0.352 0.312 SentenceCompletion WordClassification WordMeaning 0.283 0.485 0.257 Addition Code CountingDots 0.240 0.551 0.435 StraightCurvedCapitals WordRecognition NumberRecognition 0.491 0.646 0.696 FigureRecognition ObjectNumber NumberFigure 0.549 0.598 0.593 FigureWord Deduction NumericalPuzzles 0.762 0.592 0.583 ProblemReasoning SeriesCompletion ArithmeticProblems 0.601 0.497 0.500 Loadings: Factor1 Factor2 Factor3 Factor4 VisualPerception 0.160 0.689 0.187 0.160 Cubes 0.117 0.436 PaperFormBoard 0.137 0.570 0.110 Flags 0.233 0.527 GeneralInformation 0.739 0.185 0.213 0.150 PargraphComprehension 0.767 0.205 0.233 SentenceCompletion 0.806 0.197 0.153 WordClassification 0.569 0.339 0.242 0.132 WordMeaning 0.806 0.201 0.227 Addition 0.167 -0.118 0.831 0.166 Code 0.180 0.120 0.512 0.374 CountingDots 0.210 0.716 StraightCurvedCapitals 0.188 0.438 0.525 WordRecognition 0.197 0.553 NumberRecognition 0.122 0.116 0.520 FigureRecognition 0.408 0.525 ObjectNumber 0.142 0.219 0.574 NumberFigure 0.293 0.336 0.456 FigureWord 0.148 0.239 0.161 0.365 Deduction 0.378 0.402 0.118 0.301 NumericalPuzzles 0.175 0.381 0.438 0.223 ProblemReasoning 0.366 0.399 0.123 0.301 SeriesCompletion 0.369 0.500 0.244 0.239 ArithmeticProblems 0.370 0.158 0.496 0.304 Factor1 Factor2 Factor3 Factor4 SS loadings 3.647 2.872 2.657 2.290 Proportion Var 0.152 0.120 0.111 0.095 Cumulative Var 0.152 0.272 0.382 0.478 Test of the hypothesis that 4 factors are sufficient. The chi square statistic is 226.68 on 186 degrees of freedom. The p-value is 0.0224 Call: factanal(factors = factors, covmat = Harman74.cor) Uniquenesses: VisualPerception Cubes PaperFormBoard 0.450 0.781 0.639 Flags GeneralInformation PargraphComprehension 0.649 0.357 0.288 SentenceCompletion WordClassification WordMeaning 0.277 0.485 0.262 Addition Code CountingDots 0.215 0.386 0.444 StraightCurvedCapitals WordRecognition NumberRecognition 0.256 0.639 0.706 FigureRecognition ObjectNumber NumberFigure 0.550 0.614 0.596 FigureWord Deduction NumericalPuzzles 0.764 0.521 0.564 ProblemReasoning SeriesCompletion ArithmeticProblems 0.580 0.442 0.478 Loadings: Factor1 Factor2 Factor3 Factor4 Factor5 VisualPerception 0.161 0.658 0.136 0.182 0.199 Cubes 0.113 0.435 0.107 PaperFormBoard 0.135 0.562 0.107 0.116 Flags 0.231 0.533 GeneralInformation 0.736 0.188 0.192 0.162 PargraphComprehension 0.775 0.187 0.251 0.113 SentenceCompletion 0.809 0.208 0.136 WordClassification 0.568 0.348 0.223 0.131 WordMeaning 0.800 0.215 0.224 Addition 0.175 -0.100 0.844 0.176 Code 0.185 0.438 0.451 0.426 CountingDots 0.222 0.690 0.101 0.140 StraightCurvedCapitals 0.186 0.425 0.458 0.559 WordRecognition 0.197 0.557 NumberRecognition 0.121 0.130 0.508 FigureRecognition 0.400 0.529 ObjectNumber 0.145 0.208 0.562 NumberFigure 0.306 0.325 0.452 FigureWord 0.147 0.242 0.145 0.364 Deduction 0.370 0.452 0.139 0.287 -0.190 NumericalPuzzles 0.170 0.402 0.439 0.230 ProblemReasoning 0.358 0.423 0.126 0.302 SeriesCompletion 0.360 0.549 0.256 0.223 -0.107 ArithmeticProblems 0.371 0.185 0.502 0.307 Factor1 Factor2 Factor3 Factor4 Factor5 SS loadings 3.632 2.964 2.456 2.345 0.663 Proportion Var 0.151 0.124 0.102 0.098 0.028 Cumulative Var 0.151 0.275 0.377 0.475 0.503 Test of the hypothesis that 5 factors are sufficient. The chi square statistic is 186.82 on 166 degrees of freedom. The p-value is 0.128 > Harman74.FA <- factanal(factors = 5, covmat = Harman74.cor, + rotation = "promax") > print(Harman74.FA$loadings, sort = TRUE) Loadings: Factor1 Factor2 Factor3 Factor4 Factor5 VisualPerception 0.831 -0.127 0.230 Cubes 0.534 PaperFormBoard 0.736 -0.290 0.136 Flags 0.647 -0.104 SeriesCompletion 0.555 0.126 0.127 GeneralInformation 0.764 PargraphComprehension 0.845 -0.140 0.140 SentenceCompletion 0.872 -0.140 WordClassification 0.277 0.505 0.104 WordMeaning 0.846 -0.108 Addition -0.334 1.012 CountingDots 0.206 -0.200 0.722 0.185 ArithmeticProblems 0.197 0.500 0.139 WordRecognition -0.126 0.127 -0.103 0.657 NumberRecognition 0.568 FigureRecognition 0.399 -0.142 -0.207 0.562 ObjectNumber -0.108 0.107 0.613 StraightCurvedCapitals 0.542 0.247 0.618 Code 0.112 0.288 0.486 0.424 NumberFigure 0.255 -0.230 0.211 0.413 FigureWord 0.187 0.347 Deduction 0.404 0.169 0.117 -0.203 NumericalPuzzles 0.393 0.368 ProblemReasoning 0.381 0.188 0.169 Factor1 Factor2 Factor3 Factor4 Factor5 SS loadings 3.529 3.311 2.367 2.109 0.762 Proportion Var 0.147 0.138 0.099 0.088 0.032 Cumulative Var 0.147 0.285 0.384 0.471 0.503 > > > > cleanEx() > nameEx("InsectSprays") > ### * InsectSprays > > flush(stderr()); flush(stdout()) > > ### Name: InsectSprays > ### Title: Effectiveness of Insect Sprays > ### Aliases: InsectSprays > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > boxplot(count ~ spray, data = InsectSprays, + xlab = "Type of spray", ylab = "Insect count", + main = "InsectSprays data", varwidth = TRUE, col = "lightgray") > fm1 <- aov(count ~ spray, data = InsectSprays) > summary(fm1) Df Sum Sq Mean Sq F value Pr(>F) spray 5 2669 533.8 34.7 <2e-16 *** Residuals 66 1015 15.4 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0)) > plot(fm1) > fm2 <- aov(sqrt(count) ~ spray, data = InsectSprays) > summary(fm2) Df Sum Sq Mean Sq F value Pr(>F) spray 5 88.44 17.688 44.8 <2e-16 *** Residuals 66 26.06 0.395 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plot(fm2) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("JohnsonJohnson") > ### * JohnsonJohnson > > flush(stderr()); flush(stdout()) > > ### Name: JohnsonJohnson > ### Title: Quarterly Earnings per Johnson & Johnson Share > ### Aliases: JohnsonJohnson > ### Keywords: datasets > > ### ** Examples > > > cleanEx() > nameEx("LifeCycleSavings") > ### * LifeCycleSavings > > flush(stderr()); flush(stdout()) > > ### Name: LifeCycleSavings > ### Title: Intercountry Life-Cycle Savings Data > ### Aliases: LifeCycleSavings > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > pairs(LifeCycleSavings, panel = panel.smooth, + main = "LifeCycleSavings data") > fm1 <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings) > summary(fm1) Call: lm(formula = sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings) Residuals: Min 1Q Median 3Q Max -8.2422 -2.6857 -0.2488 2.4280 9.7509 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 28.5660865 7.3545161 3.884 0.000334 *** pop15 -0.4611931 0.1446422 -3.189 0.002603 ** pop75 -1.6914977 1.0835989 -1.561 0.125530 dpi -0.0003369 0.0009311 -0.362 0.719173 ddpi 0.4096949 0.1961971 2.088 0.042471 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.803 on 45 degrees of freedom Multiple R-squared: 0.3385, Adjusted R-squared: 0.2797 F-statistic: 5.756 on 4 and 45 DF, p-value: 0.0007904 > > > > cleanEx() > nameEx("Loblolly") > ### * Loblolly > > flush(stderr()); flush(stdout()) > > ### Name: Loblolly > ### Title: Growth of Loblolly pine trees > ### Aliases: Loblolly > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > plot(height ~ age, data = Loblolly, subset = Seed == 329, + xlab = "Tree age (yr)", las = 1, + ylab = "Tree height (ft)", + main = "Loblolly data and fitted curve (Seed 329 only)") > fm1 <- nls(height ~ SSasymp(age, Asym, R0, lrc), + data = Loblolly, subset = Seed == 329) > age <- seq(0, 30, length.out = 101) > lines(age, predict(fm1, list(age = age))) > > > > cleanEx() > nameEx("Nile") > ### * Nile > > flush(stderr()); flush(stdout()) > > ### Name: Nile > ### Title: Flow of the River Nile > ### Aliases: Nile > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > par(mfrow = c(2, 2)) > plot(Nile) > acf(Nile) > pacf(Nile) > ar(Nile) # selects order 2 Call: ar(x = Nile) Coefficients: 1 2 0.4081 0.1812 Order selected 2 sigma^2 estimated as 21247 > cpgram(ar(Nile)$resid) > par(mfrow = c(1, 1)) > arima(Nile, c(2, 0, 0)) Call: arima(x = Nile, order = c(2, 0, 0)) Coefficients: ar1 ar2 intercept 0.4096 0.1987 919.8397 s.e. 0.0974 0.0990 35.6410 sigma^2 estimated as 20291: log likelihood = -637.98, aic = 1283.96 > > ## Now consider missing values, following Durbin & Koopman > NileNA <- Nile > NileNA[c(21:40, 61:80)] <- NA > arima(NileNA, c(2, 0, 0)) Call: arima(x = NileNA, order = c(2, 0, 0)) Coefficients: ar1 ar2 intercept 0.3622 0.1678 918.3103 s.e. 0.1273 0.1323 39.5037 sigma^2 estimated as 23676: log likelihood = -387.7, aic = 783.41 > plot(NileNA) > pred <- + predict(arima(window(NileNA, 1871, 1890), c(2, 0, 0)), n.ahead = 20) > lines(pred$pred, lty = 3, col = "red") > lines(pred$pred + 2*pred$se, lty = 2, col = "blue") > lines(pred$pred - 2*pred$se, lty = 2, col = "blue") > pred <- + predict(arima(window(NileNA, 1871, 1930), c(2, 0, 0)), n.ahead = 20) > lines(pred$pred, lty = 3, col = "red") > lines(pred$pred + 2*pred$se, lty = 2, col = "blue") > lines(pred$pred - 2*pred$se, lty = 2, col = "blue") > > ## Structural time series models > par(mfrow = c(3, 1)) > plot(Nile) > ## local level model > (fit <- StructTS(Nile, type = "level")) Call: StructTS(x = Nile, type = "level") Variances: level epsilon 1469 15099 > lines(fitted(fit), lty = 2) # contemporaneous smoothing > lines(tsSmooth(fit), lty = 2, col = 4) # fixed-interval smoothing > plot(residuals(fit)); abline(h = 0, lty = 3) > ## local trend model > (fit2 <- StructTS(Nile, type = "trend")) ## constant trend fitted Call: StructTS(x = Nile, type = "trend") Variances: level slope epsilon 1427 0 15047 > pred <- predict(fit, n.ahead = 30) > ## with 50% confidence interval > ts.plot(Nile, pred$pred, + pred$pred + 0.67*pred$se, pred$pred -0.67*pred$se) > > ## Now consider missing values > plot(NileNA) > (fit3 <- StructTS(NileNA, type = "level")) Call: StructTS(x = NileNA, type = "level") Variances: level epsilon 685.8 17899.8 > lines(fitted(fit3), lty = 2) > lines(tsSmooth(fit3), lty = 3) > plot(residuals(fit3)); abline(h = 0, lty = 3) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("Orange") > ### * Orange > > flush(stderr()); flush(stdout()) > > ### Name: Orange > ### Title: Growth of Orange Trees > ### Aliases: Orange > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > coplot(circumference ~ age | Tree, data = Orange, show.given = FALSE) > fm1 <- nls(circumference ~ SSlogis(age, Asym, xmid, scal), + data = Orange, subset = Tree == 3) > plot(circumference ~ age, data = Orange, subset = Tree == 3, + xlab = "Tree age (days since 1968/12/31)", + ylab = "Tree circumference (mm)", las = 1, + main = "Orange tree data and fitted model (Tree 3 only)") > age <- seq(0, 1600, length.out = 101) > lines(age, predict(fm1, list(age = age))) > > > > cleanEx() > nameEx("OrchardSprays") > ### * OrchardSprays > > flush(stderr()); flush(stdout()) > > ### Name: OrchardSprays > ### Title: Potency of Orchard Sprays > ### Aliases: OrchardSprays > ### Keywords: datasets > > ### ** Examples > > require(graphics) > pairs(OrchardSprays, main = "OrchardSprays data") > > > > cleanEx() > nameEx("PlantGrowth") > ### * PlantGrowth > > flush(stderr()); flush(stdout()) > > ### Name: PlantGrowth > ### Title: Results from an Experiment on Plant Growth > ### Aliases: PlantGrowth > ### Keywords: datasets > > ### ** Examples > > ## One factor ANOVA example from Dobson's book, cf. Table 7.4: > require(stats); require(graphics) > boxplot(weight ~ group, data = PlantGrowth, main = "PlantGrowth data", + ylab = "Dried weight of plants", col = "lightgray", + notch = TRUE, varwidth = TRUE) Warning in bxp(list(stats = c(4.17, 4.53, 5.155, 5.33, 6.11, 3.59, 4.17, : some notches went outside hinges ('box'): maybe set notch=FALSE > anova(lm(weight ~ group, data = PlantGrowth)) Analysis of Variance Table Response: weight Df Sum Sq Mean Sq F value Pr(>F) group 2 3.7663 1.8832 4.8461 0.01591 * Residuals 27 10.4921 0.3886 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > > > cleanEx() > nameEx("Puromycin") > ### * Puromycin > > flush(stderr()); flush(stdout()) > > ### Name: Puromycin > ### Title: Reaction Velocity of an Enzymatic Reaction > ### Aliases: Puromycin > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## Don't show: > options(show.nls.convergence=FALSE) > ## End(Don't show) > plot(rate ~ conc, data = Puromycin, las = 1, + xlab = "Substrate concentration (ppm)", + ylab = "Reaction velocity (counts/min/min)", + pch = as.integer(Puromycin$state), + col = as.integer(Puromycin$state), + main = "Puromycin data and fitted Michaelis-Menten curves") > ## simplest form of fitting the Michaelis-Menten model to these data > fm1 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin, + subset = state == "treated", + start = c(Vm = 200, K = 0.05)) > fm2 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin, + subset = state == "untreated", + start = c(Vm = 160, K = 0.05)) > summary(fm1) Formula: rate ~ Vm * conc/(K + conc) Parameters: Estimate Std. Error t value Pr(>|t|) Vm 2.127e+02 6.947e+00 30.615 3.24e-11 *** K 6.412e-02 8.281e-03 7.743 1.57e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.93 on 10 degrees of freedom > summary(fm2) Formula: rate ~ Vm * conc/(K + conc) Parameters: Estimate Std. Error t value Pr(>|t|) Vm 1.603e+02 6.480e+00 24.734 1.38e-09 *** K 4.771e-02 7.782e-03 6.131 0.000173 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 9.773 on 9 degrees of freedom > ## add fitted lines to the plot > conc <- seq(0, 1.2, length.out = 101) > lines(conc, predict(fm1, list(conc = conc)), lty = 1, col = 1) > lines(conc, predict(fm2, list(conc = conc)), lty = 2, col = 2) > legend(0.8, 120, levels(Puromycin$state), + col = 1:2, lty = 1:2, pch = 1:2) > > ## using partial linearity > fm3 <- nls(rate ~ conc/(K + conc), data = Puromycin, + subset = state == "treated", start = c(K = 0.05), + algorithm = "plinear") > > > > cleanEx() > nameEx("Theoph") > ### * Theoph > > flush(stderr()); flush(stdout()) > > ### Name: Theoph > ### Title: Pharmacokinetics of Theophylline > ### Aliases: Theoph > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## Don't show: > options(show.nls.convergence=FALSE) > ## End(Don't show) > coplot(conc ~ Time | Subject, data = Theoph, show.given = FALSE) > Theoph.4 <- subset(Theoph, Subject == 4) > fm1 <- nls(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), + data = Theoph.4) > summary(fm1) Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) Parameters: Estimate Std. Error t value Pr(>|t|) lKe -2.4365 0.2257 -10.797 4.77e-06 *** lKa 0.1583 0.2297 0.689 0.51 lCl -3.2861 0.1448 -22.695 1.51e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.8465 on 8 degrees of freedom > plot(conc ~ Time, data = Theoph.4, + xlab = "Time since drug administration (hr)", + ylab = "Theophylline concentration (mg/L)", + main = "Observed concentrations and fitted model", + sub = "Theophylline data - Subject 4 only", + las = 1, col = 4) > xvals <- seq(0, par("usr")[2], length.out = 55) > lines(xvals, predict(fm1, newdata = list(Time = xvals)), + col = 4) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("Titanic") > ### * Titanic > > flush(stderr()); flush(stdout()) > > ### Name: Titanic > ### Title: Survival of passengers on the Titanic > ### Aliases: Titanic > ### Keywords: datasets > > ### ** Examples > > require(graphics) > mosaicplot(Titanic, main = "Survival on the Titanic") > ## Higher survival rates in children? > apply(Titanic, c(3, 4), sum) Survived Age No Yes Child 52 57 Adult 1438 654 > ## Higher survival rates in females? > apply(Titanic, c(2, 4), sum) Survived Sex No Yes Male 1364 367 Female 126 344 > ## Use loglm() in package 'MASS' for further analysis ... > > > > cleanEx() > nameEx("ToothGrowth") > ### * ToothGrowth > > flush(stderr()); flush(stdout()) > > ### Name: ToothGrowth > ### Title: The Effect of Vitamin C on Tooth Growth in Guinea Pigs > ### Aliases: ToothGrowth > ### Keywords: datasets > > ### ** Examples > > require(graphics) > coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth, + xlab = "ToothGrowth data: length vs dose, given type of supplement") > > > > cleanEx() > nameEx("UCBAdmissions") > ### * UCBAdmissions > > flush(stderr()); flush(stdout()) > > ### Name: UCBAdmissions > ### Title: Student Admissions at UC Berkeley > ### Aliases: UCBAdmissions > ### Keywords: datasets > > ### ** Examples > > require(graphics) > ## Data aggregated over departments > apply(UCBAdmissions, c(1, 2), sum) Gender Admit Male Female Admitted 1198 557 Rejected 1493 1278 > mosaicplot(apply(UCBAdmissions, c(1, 2), sum), + main = "Student admissions at UC Berkeley") > ## Data for individual departments > opar <- par(mfrow = c(2, 3), oma = c(0, 0, 2, 0)) > for(i in 1:6) + mosaicplot(UCBAdmissions[,,i], + xlab = "Admit", ylab = "Sex", + main = paste("Department", LETTERS[i])) > mtext(expression(bold("Student admissions at UC Berkeley")), + outer = TRUE, cex = 1.5) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("UKDriverDeaths") > ### * UKDriverDeaths > > flush(stderr()); flush(stdout()) > > ### Name: UKDriverDeaths > ### Title: Road Casualties in Great Britain 1969-84 > ### Aliases: UKDriverDeaths Seatbelts > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## work with pre-seatbelt period to identify a model, use logs > work <- window(log10(UKDriverDeaths), end = 1982+11/12) > par(mfrow = c(3, 1)) > plot(work); acf(work); pacf(work) > par(mfrow = c(1, 1)) > (fit <- arima(work, c(1, 0, 0), seasonal = list(order = c(1, 0, 0)))) Call: arima(x = work, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0))) Coefficients: ar1 sar1 intercept 0.4378 0.6281 3.2274 s.e. 0.0764 0.0637 0.0131 sigma^2 estimated as 0.00157: log likelihood = 300.85, aic = -593.7 > z <- predict(fit, n.ahead = 24) > ts.plot(log10(UKDriverDeaths), z$pred, z$pred+2*z$se, z$pred-2*z$se, + lty = c(1, 3, 2, 2), col = c("black", "red", "blue", "blue")) > > ## now see the effect of the explanatory variables > X <- Seatbelts[, c("kms", "PetrolPrice", "law")] > X[, 1] <- log10(X[, 1]) - 4 > arima(log10(Seatbelts[, "drivers"]), c(1, 0, 0), + seasonal = list(order = c(1, 0, 0)), xreg = X) Call: arima(x = log10(Seatbelts[, "drivers"]), order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0)), xreg = X) Coefficients: ar1 sar1 intercept kms PetrolPrice law 0.3348 0.6672 3.3539 0.0082 -1.2224 -0.0963 s.e. 0.0775 0.0612 0.0441 0.0902 0.3839 0.0166 sigma^2 estimated as 0.001476: log likelihood = 349.73, aic = -685.46 > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("UKLungDeaths") > ### * UKLungDeaths > > flush(stderr()); flush(stdout()) > > ### Name: UKLungDeaths > ### Title: Monthly Deaths from Lung Diseases in the UK > ### Aliases: UKLungDeaths ldeaths fdeaths mdeaths > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) # for time > plot(ldeaths) > plot(mdeaths, fdeaths) > ## Better labels: > yr <- floor(tt <- time(mdeaths)) > plot(mdeaths, fdeaths, + xy.labels = paste(month.abb[12*(tt - yr)], yr-1900, sep = "'")) > > > > cleanEx() > nameEx("UKgas") > ### * UKgas > > flush(stderr()); flush(stdout()) > > ### Name: UKgas > ### Title: UK Quarterly Gas Consumption > ### Aliases: UKgas > ### Keywords: datasets > > ### ** Examples > > ## maybe str(UKgas) ; plot(UKgas) ... > > > > cleanEx() > nameEx("USArrests") > ### * USArrests > > flush(stderr()); flush(stdout()) > > ### Name: USArrests > ### Title: Violent Crime Rates by US State > ### Aliases: USArrests > ### Keywords: datasets > > ### ** Examples > > summary(USArrests) Murder Assault UrbanPop Rape Min. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.30 1st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07 Median : 7.250 Median :159.0 Median :66.00 Median :20.10 Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.23 3rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18 Max. :17.400 Max. :337.0 Max. :91.00 Max. :46.00 > > require(graphics) > pairs(USArrests, panel = panel.smooth, main = "USArrests data") > > ## Difference between 'USArrests' and its correction > USArrests["Maryland", "UrbanPop"] # 67 -- the transcription error [1] 67 > UA.C <- USArrests > UA.C["Maryland", "UrbanPop"] <- 76.6 > > ## also +/- 0.5 to restore the original .5 percentages > s5u <- c("Colorado", "Florida", "Mississippi", "Wyoming") > s5d <- c("Nebraska", "Pennsylvania") > UA.C[s5u, "UrbanPop"] <- UA.C[s5u, "UrbanPop"] + 0.5 > UA.C[s5d, "UrbanPop"] <- UA.C[s5d, "UrbanPop"] - 0.5 > > ## ==> UA.C is now a *C*orrected version of USArrests > > > > cleanEx() > nameEx("USJudgeRatings") > ### * USJudgeRatings > > flush(stderr()); flush(stdout()) > > ### Name: USJudgeRatings > ### Title: Lawyers' Ratings of State Judges in the US Superior Court > ### Aliases: USJudgeRatings > ### Keywords: datasets > > ### ** Examples > > require(graphics) > pairs(USJudgeRatings, main = "USJudgeRatings data") > > > > cleanEx() > nameEx("USPersonalExpenditure") > ### * USPersonalExpenditure > > flush(stderr()); flush(stdout()) > > ### Name: USPersonalExpenditure > ### Title: Personal Expenditure Data > ### Aliases: USPersonalExpenditure > ### Keywords: datasets > > ### ** Examples > > require(stats) # for medpolish > USPersonalExpenditure 1940 1945 1950 1955 1960 Food and Tobacco 22.200 44.500 59.60 73.2 86.80 Household Operation 10.500 15.500 29.00 36.5 46.20 Medical and Health 3.530 5.760 9.71 14.0 21.10 Personal Care 1.040 1.980 2.45 3.4 5.40 Private Education 0.341 0.974 1.80 2.6 3.64 > medpolish(log10(USPersonalExpenditure)) 1: 1.126317 2: 1.032421 Final: 1.032421 Median Polish Results (Dataset: "log10(USPersonalExpenditure)") Overall: 0.9872192 Row Effects: Food and Tobacco Household Operation Medical and Health Personal Care 0.7880270 0.4327608 0.0000000 -0.5606543 Private Education -0.7319467 Column Effects: 1940 1945 1950 1955 1960 -0.4288933 -0.2267967 0.0000000 0.1423128 0.3058289 Residuals: 1940 1945 1950 1955 1960 Food and Tobacco 0.000000 0.0999105 0.000000 -0.053048 -0.142555 Household Operation 0.030103 -0.0028516 0.042418 0.000000 -0.061167 Medical and Health -0.010551 0.0000000 0.000000 0.016596 0.031234 Personal Care 0.019362 0.0968971 -0.037399 -0.037399 0.000000 Private Education -0.293625 -0.0399168 0.000000 0.017388 0.000000 > > > > cleanEx() > nameEx("VADeaths") > ### * VADeaths > > flush(stderr()); flush(stdout()) > > ### Name: VADeaths > ### Title: Death Rates in Virginia (1940) > ### Aliases: VADeaths > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > n <- length(dr <- c(VADeaths)) > nam <- names(VADeaths) > d.VAD <- data.frame( + Drate = dr, + age = rep(ordered(rownames(VADeaths)), length.out = n), + gender = gl(2, 5, n, labels = c("M", "F")), + site = gl(2, 10, labels = c("rural", "urban"))) > coplot(Drate ~ as.numeric(age) | gender * site, data = d.VAD, + panel = panel.smooth, xlab = "VADeaths data - Given: gender") > summary(aov.VAD <- aov(Drate ~ .^2, data = d.VAD)) Df Sum Sq Mean Sq F value Pr(>F) age 4 6288 1572.1 590.858 8.55e-06 *** gender 1 648 647.5 243.361 9.86e-05 *** site 1 77 76.8 28.876 0.00579 ** age:gender 4 86 21.6 8.100 0.03358 * age:site 4 43 10.6 3.996 0.10414 gender:site 1 73 73.0 27.422 0.00636 ** Residuals 4 11 2.7 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0)) > plot(aov.VAD) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("WWWusage") > ### * WWWusage > > flush(stderr()); flush(stdout()) > > ### Name: WWWusage > ### Title: Internet Usage per Minute > ### Aliases: WWWusage > ### Keywords: datasets > > ### ** Examples > > require(graphics) > work <- diff(WWWusage) > par(mfrow = c(2, 1)); plot(WWWusage); plot(work) > ## Not run: > ##D require(stats) > ##D aics <- matrix(, 6, 6, dimnames = list(p = 0:5, q = 0:5)) > ##D for(q in 1:5) aics[1, 1+q] <- arima(WWWusage, c(0, 1, q), > ##D optim.control = list(maxit = 500))$aic > ##D for(p in 1:5) > ##D for(q in 0:5) aics[1+p, 1+q] <- arima(WWWusage, c(p, 1, q), > ##D optim.control = list(maxit = 500))$aic > ##D round(aics - min(aics, na.rm = TRUE), 2) > ## End(Not run) > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("WorldPhones") > ### * WorldPhones > > flush(stderr()); flush(stdout()) > > ### Name: WorldPhones > ### Title: The World's Telephones > ### Aliases: WorldPhones > ### Keywords: datasets > > ### ** Examples > > require(graphics) > matplot(rownames(WorldPhones), WorldPhones, type = "b", log = "y", + xlab = "Year", ylab = "Number of telephones (1000's)") > legend(1951.5, 80000, colnames(WorldPhones), col = 1:6, lty = 1:5, + pch = rep(21, 7)) > title(main = "World phones data: log scale for response") > > > > cleanEx() > nameEx("ability.cov") > ### * ability.cov > > flush(stderr()); flush(stdout()) > > ### Name: ability.cov > ### Title: Ability and Intelligence Tests > ### Aliases: ability.cov > ### Keywords: datasets > > ### ** Examples > > > cleanEx() > nameEx("airmiles") > ### * airmiles > > flush(stderr()); flush(stdout()) > > ### Name: airmiles > ### Title: Passenger Miles on Commercial US Airlines, 1937-1960 > ### Aliases: airmiles > ### Keywords: datasets > > ### ** Examples > > require(graphics) > plot(airmiles, main = "airmiles data", + xlab = "Passenger-miles flown by U.S. commercial airlines", col = 4) > > > > cleanEx() > nameEx("airquality") > ### * airquality > > flush(stderr()); flush(stdout()) > > ### Name: airquality > ### Title: New York Air Quality Measurements > ### Aliases: airquality > ### Keywords: datasets > > ### ** Examples > > require(graphics) > pairs(airquality, panel = panel.smooth, main = "airquality data") > > > > cleanEx() > nameEx("anscombe") > ### * anscombe > > flush(stderr()); flush(stdout()) > > ### Name: anscombe > ### Title: Anscombe's Quartet of 'Identical' Simple Linear Regressions > ### Aliases: anscombe > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > summary(anscombe) x1 x2 x3 x4 y1 Min. : 4.0 Min. : 4.0 Min. : 4.0 Min. : 8 Min. : 4.260 1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 8 1st Qu.: 6.315 Median : 9.0 Median : 9.0 Median : 9.0 Median : 8 Median : 7.580 Mean : 9.0 Mean : 9.0 Mean : 9.0 Mean : 9 Mean : 7.501 3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.: 8 3rd Qu.: 8.570 Max. :14.0 Max. :14.0 Max. :14.0 Max. :19 Max. :10.840 y2 y3 y4 Min. :3.100 Min. : 5.39 Min. : 5.250 1st Qu.:6.695 1st Qu.: 6.25 1st Qu.: 6.170 Median :8.140 Median : 7.11 Median : 7.040 Mean :7.501 Mean : 7.50 Mean : 7.501 3rd Qu.:8.950 3rd Qu.: 7.98 3rd Qu.: 8.190 Max. :9.260 Max. :12.74 Max. :12.500 > > ##-- now some "magic" to do the 4 regressions in a loop: > ff <- y ~ x > mods <- setNames(as.list(1:4), paste0("lm", 1:4)) > for(i in 1:4) { + ff[2:3] <- lapply(paste0(c("y","x"), i), as.name) + ## or ff[[2]] <- as.name(paste0("y", i)) + ## ff[[3]] <- as.name(paste0("x", i)) + mods[[i]] <- lmi <- lm(ff, data = anscombe) + print(anova(lmi)) + } Analysis of Variance Table Response: y1 Df Sum Sq Mean Sq F value Pr(>F) x1 1 27.510 27.5100 17.99 0.00217 ** Residuals 9 13.763 1.5292 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Analysis of Variance Table Response: y2 Df Sum Sq Mean Sq F value Pr(>F) x2 1 27.500 27.5000 17.966 0.002179 ** Residuals 9 13.776 1.5307 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Analysis of Variance Table Response: y3 Df Sum Sq Mean Sq F value Pr(>F) x3 1 27.470 27.4700 17.972 0.002176 ** Residuals 9 13.756 1.5285 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Analysis of Variance Table Response: y4 Df Sum Sq Mean Sq F value Pr(>F) x4 1 27.490 27.4900 18.003 0.002165 ** Residuals 9 13.742 1.5269 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > ## See how close they are (numerically!) > sapply(mods, coef) lm1 lm2 lm3 lm4 (Intercept) 3.0000909 3.000909 3.0024545 3.0017273 x1 0.5000909 0.500000 0.4997273 0.4999091 > lapply(mods, function(fm) coef(summary(fm))) $lm1 Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0000909 1.1247468 2.667348 0.025734051 x1 0.5000909 0.1179055 4.241455 0.002169629 $lm2 Estimate Std. Error t value Pr(>|t|) (Intercept) 3.000909 1.1253024 2.666758 0.025758941 x2 0.500000 0.1179637 4.238590 0.002178816 $lm3 Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0024545 1.1244812 2.670080 0.025619109 x3 0.4997273 0.1178777 4.239372 0.002176305 $lm4 Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0017273 1.1239211 2.670763 0.025590425 x4 0.4999091 0.1178189 4.243028 0.002164602 > > ## Now, do what you should have done in the first place: PLOTS > op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma = c(0, 0, 2, 0)) > for(i in 1:4) { + ff[2:3] <- lapply(paste0(c("y","x"), i), as.name) + plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2, + xlim = c(3, 19), ylim = c(3, 13)) + abline(mods[[i]], col = "blue") + } > mtext("Anscombe's 4 Regression data sets", outer = TRUE, cex = 1.5) > par(op) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("attenu") > ### * attenu > > flush(stderr()); flush(stdout()) > > ### Name: attenu > ### Title: The Joyner-Boore Attenuation Data > ### Aliases: attenu > ### Keywords: datasets > > ### ** Examples > > require(graphics) > ## check the data class of the variables > sapply(attenu, data.class) event mag station dist accel "numeric" "numeric" "factor" "numeric" "numeric" > summary(attenu) event mag station dist Min. : 1.00 Min. :5.000 117 : 5 Min. : 0.50 1st Qu.: 9.00 1st Qu.:5.300 1028 : 4 1st Qu.: 11.32 Median :18.00 Median :6.100 113 : 4 Median : 23.40 Mean :14.74 Mean :6.084 112 : 3 Mean : 45.60 3rd Qu.:20.00 3rd Qu.:6.600 135 : 3 3rd Qu.: 47.55 Max. :23.00 Max. :7.700 (Other):147 Max. :370.00 NA's : 16 accel Min. :0.00300 1st Qu.:0.04425 Median :0.11300 Mean :0.15422 3rd Qu.:0.21925 Max. :0.81000 > pairs(attenu, main = "attenu data") > coplot(accel ~ dist | as.factor(event), data = attenu, show.given = FALSE) > coplot(log(accel) ~ log(dist) | as.factor(event), + data = attenu, panel = panel.smooth, show.given = FALSE) > > > > cleanEx() > nameEx("attitude") > ### * attitude > > flush(stderr()); flush(stdout()) > > ### Name: attitude > ### Title: The Chatterjee-Price Attitude Data > ### Aliases: attitude > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > pairs(attitude, main = "attitude data") > summary(attitude) rating complaints privileges learning raises Min. :40.00 Min. :37.0 Min. :30.00 Min. :34.00 Min. :43.00 1st Qu.:58.75 1st Qu.:58.5 1st Qu.:45.00 1st Qu.:47.00 1st Qu.:58.25 Median :65.50 Median :65.0 Median :51.50 Median :56.50 Median :63.50 Mean :64.63 Mean :66.6 Mean :53.13 Mean :56.37 Mean :64.63 3rd Qu.:71.75 3rd Qu.:77.0 3rd Qu.:62.50 3rd Qu.:66.75 3rd Qu.:71.00 Max. :85.00 Max. :90.0 Max. :83.00 Max. :75.00 Max. :88.00 critical advance Min. :49.00 Min. :25.00 1st Qu.:69.25 1st Qu.:35.00 Median :77.50 Median :41.00 Mean :74.77 Mean :42.93 3rd Qu.:80.00 3rd Qu.:47.75 Max. :92.00 Max. :72.00 > summary(fm1 <- lm(rating ~ ., data = attitude)) Call: lm(formula = rating ~ ., data = attitude) Residuals: Min 1Q Median 3Q Max -10.9418 -4.3555 0.3158 5.5425 11.5990 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.78708 11.58926 0.931 0.361634 complaints 0.61319 0.16098 3.809 0.000903 *** privileges -0.07305 0.13572 -0.538 0.595594 learning 0.32033 0.16852 1.901 0.069925 . raises 0.08173 0.22148 0.369 0.715480 critical 0.03838 0.14700 0.261 0.796334 advance -0.21706 0.17821 -1.218 0.235577 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 7.068 on 23 degrees of freedom Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628 F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0), + mar = c(4.1, 4.1, 2.1, 1.1)) > plot(fm1) > summary(fm2 <- lm(rating ~ complaints, data = attitude)) Call: lm(formula = rating ~ complaints, data = attitude) Residuals: Min 1Q Median 3Q Max -12.8799 -5.9905 0.1783 6.2978 9.6294 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 14.37632 6.61999 2.172 0.0385 * complaints 0.75461 0.09753 7.737 1.99e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.993 on 28 degrees of freedom Multiple R-squared: 0.6813, Adjusted R-squared: 0.6699 F-statistic: 59.86 on 1 and 28 DF, p-value: 1.988e-08 > plot(fm2) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("beavers") > ### * beavers > > flush(stderr()); flush(stdout()) > > ### Name: beavers > ### Title: Body Temperature Series of Two Beavers > ### Aliases: beavers beaver1 beaver2 > ### Keywords: datasets > > ### ** Examples > > require(graphics) > (yl <- range(beaver1$temp, beaver2$temp)) [1] 36.33 38.35 > > beaver.plot <- function(bdat, ...) { + nam <- deparse(substitute(bdat)) + with(bdat, { + # Hours since start of day: + hours <- time %/% 100 + 24*(day - day[1]) + (time %% 100)/60 + plot (hours, temp, type = "l", ..., + main = paste(nam, "body temperature")) + abline(h = 37.5, col = "gray", lty = 2) + is.act <- activ == 1 + points(hours[is.act], temp[is.act], col = 2, cex = .8) + }) + } > op <- par(mfrow = c(2, 1), mar = c(3, 3, 4, 2), mgp = 0.9 * 2:0) > beaver.plot(beaver1, ylim = yl) > beaver.plot(beaver2, ylim = yl) > par(op) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("cars") > ### * cars > > flush(stderr()); flush(stdout()) > > ### Name: cars > ### Title: Speed and Stopping Distances of Cars > ### Aliases: cars > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)", + las = 1) > lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red") > title(main = "cars data") > plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)", + las = 1, log = "xy") > title(main = "cars data (logarithmic scales)") > lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red") > summary(fm1 <- lm(log(dist) ~ log(speed), data = cars)) Call: lm(formula = log(dist) ~ log(speed), data = cars) Residuals: Min 1Q Median 3Q Max -1.00215 -0.24578 -0.02898 0.20717 0.88289 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.7297 0.3758 -1.941 0.0581 . log(speed) 1.6024 0.1395 11.484 2.26e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4053 on 48 degrees of freedom Multiple R-squared: 0.7331, Adjusted R-squared: 0.7276 F-statistic: 131.9 on 1 and 48 DF, p-value: 2.259e-15 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0), + mar = c(4.1, 4.1, 2.1, 1.1)) > plot(fm1) > par(opar) > > ## An example of polynomial regression > plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)", + las = 1, xlim = c(0, 25)) > d <- seq(0, 25, length.out = 200) > for(degree in 1:4) { + fm <- lm(dist ~ poly(speed, degree), data = cars) + assign(paste("cars", degree, sep = "."), fm) + lines(d, predict(fm, data.frame(speed = d)), col = degree) + } > anova(cars.1, cars.2, cars.3, cars.4) Analysis of Variance Table Model 1: dist ~ poly(speed, degree) Model 2: dist ~ poly(speed, degree) Model 3: dist ~ poly(speed, degree) Model 4: dist ~ poly(speed, degree) Res.Df RSS Df Sum of Sq F Pr(>F) 1 48 11354 2 47 10825 1 528.81 2.3108 0.1355 3 46 10634 1 190.35 0.8318 0.3666 4 45 10298 1 336.55 1.4707 0.2316 > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("chickwts") > ### * chickwts > > flush(stderr()); flush(stdout()) > > ### Name: chickwts > ### Title: Chicken Weights by Feed Type > ### Aliases: chickwts > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > boxplot(weight ~ feed, data = chickwts, col = "lightgray", + varwidth = TRUE, notch = TRUE, main = "chickwt data", + ylab = "Weight at six weeks (gm)") Warning in bxp(list(stats = c(216, 271.5, 342, 373.5, 404, 108, 136, 151.5, : some notches went outside hinges ('box'): maybe set notch=FALSE > anova(fm1 <- lm(weight ~ feed, data = chickwts)) Analysis of Variance Table Response: weight Df Sum Sq Mean Sq F value Pr(>F) feed 5 231129 46226 15.365 5.936e-10 *** Residuals 65 195556 3009 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0), + mar = c(4.1, 4.1, 2.1, 1.1)) > plot(fm1) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("co2") > ### * co2 > > flush(stderr()); flush(stdout()) > > ### Name: co2 > ### Title: Mauna Loa Atmospheric CO2 Concentration > ### Aliases: co2 > ### Keywords: datasets > > ### ** Examples > > require(graphics) > plot(co2, ylab = expression("Atmospheric concentration of CO"[2]), + las = 1) > title(main = "co2 data set") > > > > cleanEx() > nameEx("crimtab") > ### * crimtab > > flush(stderr()); flush(stdout()) > > ### Name: crimtab > ### Title: Student's 3000 Criminals Data > ### Aliases: crimtab > ### Keywords: datasets > > ### ** Examples > > require(stats) > dim(crimtab) [1] 42 22 > utils::str(crimtab) 'table' int [1:42, 1:22] 0 0 0 0 0 0 1 0 0 0 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:42] "9.4" "9.5" "9.6" "9.7" ... ..$ : chr [1:22] "142.24" "144.78" "147.32" "149.86" ... > ## for nicer printing: > local({cT <- crimtab + colnames(cT) <- substring(colnames(cT), 2, 3) + print(cT, zero.print = " ") + }) 42 44 47 49 52 54 57 60 62 65 67 70 72 75 77 80 82 85 87 90 93 95 9.4 9.5 1 9.6 9.7 9.8 1 9.9 1 1 1 10 1 1 2 2 1 10.1 1 3 1 1 1 10.2 2 2 2 1 2 1 10.3 1 1 3 2 2 3 5 10.4 1 1 2 3 3 4 3 3 10.5 1 3 7 6 4 3 1 3 1 1 10.6 1 4 5 9 14 6 3 1 1 10.7 1 2 4 9 14 16 15 7 3 1 2 10.8 2 5 6 14 27 10 7 1 2 1 10.9 2 6 14 24 27 14 10 4 1 11 2 6 12 15 31 37 27 17 10 6 11.1 3 3 12 22 26 24 26 24 7 4 1 11.2 3 2 7 21 30 38 29 27 20 4 1 1 11.3 1 5 10 24 26 39 26 24 7 2 11.4 3 4 9 29 56 58 26 22 10 11 11.5 5 11 17 33 57 38 34 25 11 2 11.6 2 1 4 13 37 39 48 38 27 12 2 2 1 11.7 2 9 17 30 37 48 45 24 9 9 2 11.8 1 2 11 15 35 41 34 29 10 5 1 11.9 1 1 2 12 10 27 32 35 19 10 9 3 1 12 1 4 8 19 42 39 22 16 8 2 2 12.1 2 4 13 22 28 15 27 10 4 1 12.2 1 2 5 6 23 17 16 11 8 1 1 12.3 4 8 10 13 20 23 6 5 12.4 1 1 1 2 7 12 4 7 7 1 1 12.5 1 1 3 12 11 8 6 8 2 12.6 1 3 5 7 8 6 3 1 1 12.7 1 1 7 5 5 8 2 2 12.8 1 2 3 1 8 5 3 1 1 12.9 1 2 2 1 1 13 3 1 1 2 1 13.1 1 1 13.2 1 1 1 3 13.3 1 1 13.4 13.5 1 > > ## Repeat Student's experiment: > > # 1) Reconstitute 3000 raw data for heights in inches and rounded to > # nearest integer as in Student's paper: > > (heIn <- round(as.numeric(colnames(crimtab)) / 2.54)) [1] 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 > d.hei <- data.frame(height = rep(heIn, colSums(crimtab))) > > # 2) shuffle the data: > > set.seed(1) > d.hei <- d.hei[sample(1:3000), , drop = FALSE] > > # 3) Make 750 samples each of size 4: > > d.hei$sample <- as.factor(rep(1:750, each = 4)) > > # 4) Compute the means and standard deviations (n) for the 750 samples: > > h.mean <- with(d.hei, tapply(height, sample, FUN = mean)) > h.sd <- with(d.hei, tapply(height, sample, FUN = sd)) * sqrt(3/4) > > # 5) Compute the difference between the mean of each sample and > # the mean of the population and then divide by the > # standard deviation of the sample: > > zobs <- (h.mean - mean(d.hei[,"height"]))/h.sd > > # 6) Replace infinite values by +/- 6 as in Student's paper: > > zobs[infZ <- is.infinite(zobs)] # none of them named numeric(0) > zobs[infZ] <- 6 * sign(zobs[infZ]) > > # 7) Plot the distribution: > > require(grDevices); require(graphics) > hist(x = zobs, probability = TRUE, xlab = "Student's z", + col = grey(0.8), border = grey(0.5), + main = "Distribution of Student's z score for 'crimtab' data") > > > > cleanEx() > nameEx("discoveries") > ### * discoveries > > flush(stderr()); flush(stdout()) > > ### Name: discoveries > ### Title: Yearly Numbers of Important Discoveries > ### Aliases: discoveries > ### Keywords: datasets > > ### ** Examples > > require(graphics) > plot(discoveries, ylab = "Number of important discoveries", + las = 1) > title(main = "discoveries data set") > > > > cleanEx() > nameEx("esoph") > ### * esoph > > flush(stderr()); flush(stdout()) > > ### Name: esoph > ### Title: Smoking, Alcohol and (O)esophageal Cancer > ### Aliases: esoph > ### Keywords: datasets > > ### ** Examples > > require(stats) > require(graphics) # for mosaicplot > summary(esoph) agegp alcgp tobgp ncases ncontrols 25-34:15 0-39g/day:23 0-9g/day:24 Min. : 0.000 Min. : 1.00 35-44:15 40-79 :23 10-19 :24 1st Qu.: 0.000 1st Qu.: 3.00 45-54:16 80-119 :21 20-29 :20 Median : 1.000 Median : 6.00 55-64:16 120+ :21 30+ :20 Mean : 2.273 Mean :11.08 65-74:15 3rd Qu.: 4.000 3rd Qu.:14.00 75+ :11 Max. :17.000 Max. :60.00 > ## effects of alcohol, tobacco and interaction, age-adjusted > model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp, + data = esoph, family = binomial()) > anova(model1) Analysis of Deviance Table Model: binomial, link: logit Response: cbind(ncases, ncontrols) Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 87 227.241 agegp 5 88.128 82 139.112 tobgp 3 19.085 79 120.028 alcgp 3 66.054 76 53.973 tobgp:alcgp 9 6.489 67 47.484 > ## Try a linear effect of alcohol and tobacco > model2 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) + + unclass(alcgp), + data = esoph, family = binomial()) > summary(model2) Call: glm(formula = cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) + unclass(alcgp), family = binomial(), data = esoph) Deviance Residuals: Min 1Q Median 3Q Max -1.7628 -0.6426 -0.2709 0.3043 2.0421 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.01097 0.31224 -12.846 < 2e-16 *** agegp.L 2.96113 0.65092 4.549 5.39e-06 *** agegp.Q -1.33735 0.58918 -2.270 0.02322 * agegp.C 0.15292 0.44792 0.341 0.73281 agegp^4 0.06668 0.30776 0.217 0.82848 agegp^5 -0.20288 0.19523 -1.039 0.29872 unclass(tobgp) 0.26162 0.08198 3.191 0.00142 ** unclass(alcgp) 0.65308 0.08452 7.727 1.10e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 227.241 on 87 degrees of freedom Residual deviance: 59.277 on 80 degrees of freedom AIC: 222.76 Number of Fisher Scoring iterations: 6 > ## Re-arrange data for a mosaic plot > ttt <- table(esoph$agegp, esoph$alcgp, esoph$tobgp) > o <- with(esoph, order(tobgp, alcgp, agegp)) > ttt[ttt == 1] <- esoph$ncases[o] > tt1 <- table(esoph$agegp, esoph$alcgp, esoph$tobgp) > tt1[tt1 == 1] <- esoph$ncontrols[o] > tt <- array(c(ttt, tt1), c(dim(ttt),2), + c(dimnames(ttt), list(c("Cancer", "control")))) > mosaicplot(tt, main = "esoph data set", color = TRUE) > > > > cleanEx() > nameEx("euro") > ### * euro > > flush(stderr()); flush(stdout()) > > ### Name: euro > ### Title: Conversion Rates of Euro Currencies > ### Aliases: euro euro.cross > ### Keywords: datasets > > ### ** Examples > > cbind(euro) euro ATS 13.760300 BEF 40.339900 DEM 1.955830 ESP 166.386000 FIM 5.945730 FRF 6.559570 IEP 0.787564 ITL 1936.270000 LUF 40.339900 NLG 2.203710 PTE 200.482000 > > ## These relations hold: > euro == signif(euro, 6) # [6 digit precision in Euro's definition] ATS BEF DEM ESP FIM FRF IEP ITL LUF NLG PTE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE > all(euro.cross == outer(1/euro, euro)) [1] TRUE > > ## Convert 20 Euro to Belgian Franc > 20 * euro["BEF"] BEF 806.798 > ## Convert 20 Austrian Schilling to Euro > 20 / euro["ATS"] ATS 1.453457 > ## Convert 20 Spanish Pesetas to Italian Lira > 20 * euro.cross["ESP", "ITL"] [1] 232.7443 > > require(graphics) > dotchart(euro, + main = "euro data: 1 Euro in currency unit") > dotchart(1/euro, + main = "euro data: 1 currency unit in Euros") > dotchart(log(euro, 10), + main = "euro data: log10(1 Euro in currency unit)") > > > > cleanEx() > nameEx("faithful") > ### * faithful > > flush(stderr()); flush(stdout()) > > ### Name: faithful > ### Title: Old Faithful Geyser Data > ### Aliases: faithful > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > f.tit <- "faithful data: Eruptions of Old Faithful" > > ne60 <- round(e60 <- 60 * faithful$eruptions) > all.equal(e60, ne60) # relative diff. ~ 1/10000 [1] "Mean relative difference: 9.515332e-05" > table(zapsmall(abs(e60 - ne60))) # 0, 0.02 or 0.04 0 0.02 0.04 106 163 3 > faithful$better.eruptions <- ne60 / 60 > te <- table(ne60) > te[te >= 4] # (too) many multiples of 5 ! ne60 105 108 110 112 113 120 216 230 240 245 249 250 255 260 261 262 265 270 272 275 6 4 7 8 4 4 4 5 6 5 4 4 4 5 4 4 4 8 5 4 276 282 288 4 6 6 > plot(names(te), te, type = "h", main = f.tit, xlab = "Eruption time (sec)") > > plot(faithful[, -3], main = f.tit, + xlab = "Eruption time (min)", + ylab = "Waiting time to next eruption (min)") > lines(lowess(faithful$eruptions, faithful$waiting, f = 2/3, iter = 3), + col = "red") > > > > cleanEx() > nameEx("freeny") > ### * freeny > > flush(stderr()); flush(stdout()) > > ### Name: freeny > ### Title: Freeny's Revenue Data > ### Aliases: freeny freeny.x freeny.y > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > summary(freeny) y lag.quarterly.revenue price.index income.level Min. :8.791 Min. :8.791 Min. :4.278 Min. :5.821 1st Qu.:9.045 1st Qu.:9.020 1st Qu.:4.392 1st Qu.:5.948 Median :9.314 Median :9.284 Median :4.510 Median :6.061 Mean :9.306 Mean :9.281 Mean :4.496 Mean :6.039 3rd Qu.:9.591 3rd Qu.:9.561 3rd Qu.:4.605 3rd Qu.:6.139 Max. :9.794 Max. :9.775 Max. :4.710 Max. :6.200 market.potential Min. :12.97 1st Qu.:13.01 Median :13.07 Mean :13.07 3rd Qu.:13.12 Max. :13.17 > pairs(freeny, main = "freeny data") > # gives warning: freeny$y has class "ts" > > summary(fm1 <- lm(y ~ ., data = freeny)) Call: lm(formula = y ~ ., data = freeny) Residuals: Min 1Q Median 3Q Max -0.0259426 -0.0101033 0.0003824 0.0103236 0.0267124 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -10.4726 6.0217 -1.739 0.0911 . lag.quarterly.revenue 0.1239 0.1424 0.870 0.3904 price.index -0.7542 0.1607 -4.693 4.28e-05 *** income.level 0.7675 0.1339 5.730 1.93e-06 *** market.potential 1.3306 0.5093 2.613 0.0133 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.01473 on 34 degrees of freedom Multiple R-squared: 0.9981, Adjusted R-squared: 0.9978 F-statistic: 4354 on 4 and 34 DF, p-value: < 2.2e-16 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0), + mar = c(4.1, 4.1, 2.1, 1.1)) > plot(fm1) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("infert") > ### * infert > > flush(stderr()); flush(stdout()) > > ### Name: infert > ### Title: Infertility after Spontaneous and Induced Abortion > ### Aliases: infert > ### Keywords: datasets > > ### ** Examples > > require(stats) > model1 <- glm(case ~ spontaneous+induced, data = infert, family = binomial()) > summary(model1) Call: glm(formula = case ~ spontaneous + induced, family = binomial(), data = infert) Deviance Residuals: Min 1Q Median 3Q Max -1.6678 -0.8360 -0.5772 0.9030 1.9362 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.7079 0.2677 -6.380 1.78e-10 *** spontaneous 1.1972 0.2116 5.657 1.54e-08 *** induced 0.4181 0.2056 2.033 0.042 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 316.17 on 247 degrees of freedom Residual deviance: 279.61 on 245 degrees of freedom AIC: 285.61 Number of Fisher Scoring iterations: 4 > ## adjusted for other potential confounders: > summary(model2 <- glm(case ~ age+parity+education+spontaneous+induced, + data = infert, family = binomial())) Call: glm(formula = case ~ age + parity + education + spontaneous + induced, family = binomial(), data = infert) Deviance Residuals: Min 1Q Median 3Q Max -1.7603 -0.8162 -0.4956 0.8349 2.6536 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.14924 1.41220 -0.814 0.4158 age 0.03958 0.03120 1.269 0.2046 parity -0.82828 0.19649 -4.215 2.49e-05 *** education6-11yrs -1.04424 0.79255 -1.318 0.1876 education12+ yrs -1.40321 0.83416 -1.682 0.0925 . spontaneous 2.04591 0.31016 6.596 4.21e-11 *** induced 1.28876 0.30146 4.275 1.91e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 316.17 on 247 degrees of freedom Residual deviance: 257.80 on 241 degrees of freedom AIC: 271.8 Number of Fisher Scoring iterations: 4 > ## Really should be analysed by conditional logistic regression > ## which is in the survival package > > > > cleanEx() > nameEx("iris") > ### * iris > > flush(stderr()); flush(stdout()) > > ### Name: iris > ### Title: Edgar Anderson's Iris Data > ### Aliases: iris iris3 > ### Keywords: datasets > > ### ** Examples > > dni3 <- dimnames(iris3) > ii <- data.frame(matrix(aperm(iris3, c(1,3,2)), ncol = 4, + dimnames = list(NULL, sub(" L.",".Length", + sub(" W.",".Width", dni3[[2]])))), + Species = gl(3, 50, labels = sub("S", "s", sub("V", "v", dni3[[3]])))) > all.equal(ii, iris) # TRUE [1] TRUE > > > > cleanEx() > nameEx("islands") > ### * islands > > flush(stderr()); flush(stdout()) > > ### Name: islands > ### Title: Areas of the World's Major Landmasses > ### Aliases: islands > ### Keywords: datasets > > ### ** Examples > > require(graphics) > dotchart(log(islands, 10), + main = "islands data: log10(area) (log10(sq. miles))") > dotchart(log(islands[order(islands)], 10), + main = "islands data: log10(area) (log10(sq. miles))") > > > > cleanEx() > nameEx("longley") > ### * longley > > flush(stderr()); flush(stdout()) > > ### Name: longley > ### Title: Longley's Economic Regression Data > ### Aliases: longley > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## give the data set in the form it is used in S-PLUS: > longley.x <- data.matrix(longley[, 1:6]) > longley.y <- longley[, "Employed"] > pairs(longley, main = "longley data") > summary(fm1 <- lm(Employed ~ ., data = longley)) Call: lm(formula = Employed ~ ., data = longley) Residuals: Min 1Q Median 3Q Max -0.41011 -0.15767 -0.02816 0.10155 0.45539 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.482e+03 8.904e+02 -3.911 0.003560 ** GNP.deflator 1.506e-02 8.492e-02 0.177 0.863141 GNP -3.582e-02 3.349e-02 -1.070 0.312681 Unemployed -2.020e-02 4.884e-03 -4.136 0.002535 ** Armed.Forces -1.033e-02 2.143e-03 -4.822 0.000944 *** Population -5.110e-02 2.261e-01 -0.226 0.826212 Year 1.829e+00 4.555e-01 4.016 0.003037 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.3049 on 9 degrees of freedom Multiple R-squared: 0.9955, Adjusted R-squared: 0.9925 F-statistic: 330.3 on 6 and 9 DF, p-value: 4.984e-10 > opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0), + mar = c(4.1, 4.1, 2.1, 1.1)) > plot(fm1) > par(opar) > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("morley") > ### * morley > > flush(stderr()); flush(stdout()) > > ### Name: morley > ### Title: Michelson Speed of Light Data > ### Aliases: morley > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > michelson <- transform(morley, + Expt = factor(Expt), Run = factor(Run)) > xtabs(~ Expt + Run, data = michelson) # 5 x 20 balanced (two-way) Run Expt 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 > plot(Speed ~ Expt, data = michelson, + main = "Speed of Light Data", xlab = "Experiment No.") > fm <- aov(Speed ~ Run + Expt, data = michelson) > summary(fm) Df Sum Sq Mean Sq F value Pr(>F) Run 19 113344 5965 1.105 0.36321 Expt 4 94514 23629 4.378 0.00307 ** Residuals 76 410166 5397 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > fm0 <- update(fm, . ~ . - Run) > anova(fm0, fm) Analysis of Variance Table Model 1: Speed ~ Expt Model 2: Speed ~ Run + Expt Res.Df RSS Df Sum of Sq F Pr(>F) 1 95 523510 2 76 410166 19 113344 1.1053 0.3632 > > > > cleanEx() > nameEx("mtcars") > ### * mtcars > > flush(stderr()); flush(stdout()) > > ### Name: mtcars > ### Title: Motor Trend Car Road Tests > ### Aliases: mtcars > ### Keywords: datasets > > ### ** Examples > > require(graphics) > pairs(mtcars, main = "mtcars data", gap = 1/4) > coplot(mpg ~ disp | as.factor(cyl), data = mtcars, + panel = panel.smooth, rows = 1) > ## possibly more meaningful, e.g., for summary() or bivariate plots: > mtcars2 <- within(mtcars, { + vs <- factor(vs, labels = c("V", "S")) + am <- factor(am, labels = c("automatic", "manual")) + cyl <- ordered(cyl) + gear <- ordered(gear) + carb <- ordered(carb) + }) > summary(mtcars2) mpg cyl disp hp drat Min. :10.40 4:11 Min. : 71.1 Min. : 52.0 Min. :2.760 1st Qu.:15.43 6: 7 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080 Median :19.20 8:14 Median :196.3 Median :123.0 Median :3.695 Mean :20.09 Mean :230.7 Mean :146.7 Mean :3.597 3rd Qu.:22.80 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920 Max. :33.90 Max. :472.0 Max. :335.0 Max. :4.930 wt qsec vs am gear carb Min. :1.513 Min. :14.50 V:18 automatic:19 3:15 1: 7 1st Qu.:2.581 1st Qu.:16.89 S:14 manual :13 4:12 2:10 Median :3.325 Median :17.71 5: 5 3: 3 Mean :3.217 Mean :17.85 4:10 3rd Qu.:3.610 3rd Qu.:18.90 6: 1 Max. :5.424 Max. :22.90 8: 1 > > > > cleanEx() > nameEx("nhtemp") > ### * nhtemp > > flush(stderr()); flush(stdout()) > > ### Name: nhtemp > ### Title: Average Yearly Temperatures in New Haven > ### Aliases: nhtemp > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > plot(nhtemp, main = "nhtemp data", + ylab = "Mean annual temperature in New Haven, CT (deg. F)") > > > > cleanEx() > nameEx("nottem") > ### * nottem > > flush(stderr()); flush(stdout()) > > ### Name: nottem > ### Title: Average Monthly Temperatures at Nottingham, 1920-1939 > ### Aliases: nottem > ### Keywords: datasets > > ### ** Examples > > > > cleanEx() > nameEx("npk") > ### * npk > > flush(stderr()); flush(stdout()) > > ### Name: npk > ### Title: Classical N, P, K Factorial Experiment > ### Aliases: npk > ### Keywords: datasets > > ### ** Examples > > > cleanEx() > nameEx("occupationalStatus") > ### * occupationalStatus > > flush(stderr()); flush(stdout()) > > ### Name: occupationalStatus > ### Title: Occupational Status of Fathers and their Sons > ### Aliases: occupationalStatus > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > > plot(occupationalStatus) > > ## Fit a uniform association model separating diagonal effects > Diag <- as.factor(diag(1:8)) > Rscore <- scale(as.numeric(row(occupationalStatus)), scale = FALSE) > Cscore <- scale(as.numeric(col(occupationalStatus)), scale = FALSE) > modUnif <- glm(Freq ~ origin + destination + Diag + Rscore:Cscore, + family = poisson, data = occupationalStatus) > > summary(modUnif) Call: glm(formula = Freq ~ origin + destination + Diag + Rscore:Cscore, family = poisson, data = occupationalStatus) Deviance Residuals: Min 1Q Median 3Q Max -2.6521 -0.6267 0.0000 0.5913 2.0964 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.568592 0.183358 3.101 0.001929 ** origin2 0.431314 0.149415 2.887 0.003893 ** origin3 1.461862 0.131141 11.147 < 2e-16 *** origin4 1.788731 0.126588 14.130 < 2e-16 *** origin5 0.441011 0.144844 3.045 0.002329 ** origin6 2.491735 0.121219 20.556 < 2e-16 *** origin7 1.127536 0.129032 8.738 < 2e-16 *** origin8 0.796445 0.131863 6.040 1.54e-09 *** destination2 0.873202 0.166902 5.232 1.68e-07 *** destination3 1.813992 0.153861 11.790 < 2e-16 *** destination4 2.082515 0.150887 13.802 < 2e-16 *** destination5 1.366383 0.155590 8.782 < 2e-16 *** destination6 2.816369 0.146054 19.283 < 2e-16 *** destination7 1.903918 0.147810 12.881 < 2e-16 *** destination8 1.398585 0.151658 9.222 < 2e-16 *** Diag1 1.665495 0.237383 7.016 2.28e-12 *** Diag2 0.959681 0.212122 4.524 6.06e-06 *** Diag3 0.021750 0.156530 0.139 0.889490 Diag4 0.226399 0.124364 1.820 0.068689 . Diag5 0.808646 0.229754 3.520 0.000432 *** Diag6 0.132277 0.077191 1.714 0.086597 . Diag7 0.506709 0.115936 4.371 1.24e-05 *** Diag8 0.221880 0.134803 1.646 0.099771 . Rscore:Cscore 0.136974 0.007489 18.289 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 4679.004 on 63 degrees of freedom Residual deviance: 58.436 on 40 degrees of freedom AIC: 428.78 Number of Fisher Scoring iterations: 4 > plot(modUnif) # 4 plots, with warning about h_ii ~= 1 Warning: not plotting observations with leverage one: 1, 10, 19, 28, 37, 46, 55, 64 Warning: not plotting observations with leverage one: 1, 10, 19, 28, 37, 46, 55, 64 > > > > cleanEx() > nameEx("precip") > ### * precip > > flush(stderr()); flush(stdout()) > > ### Name: precip > ### Title: Annual Precipitation in US Cities > ### Aliases: precip > ### Keywords: datasets > > ### ** Examples > > require(graphics) > dotchart(precip[order(precip)], main = "precip data") > title(sub = "Average annual precipitation (in.)") > > ## Old ("wrong") version of dataset (just name change): > precip.O <- local({ + p <- precip; names(p)[names(p) == "Cincinnati"] <- "Cincinati" ; p }) > stopifnot(all(precip == precip.O), + match("Cincinnati", names(precip)) == 46, + identical(names(precip)[-46], names(precip.O)[-46])) > > > > cleanEx() > nameEx("presidents") > ### * presidents > > flush(stderr()); flush(stdout()) > > ### Name: presidents > ### Title: Quarterly Approval Ratings of US Presidents > ### Aliases: presidents > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > plot(presidents, las = 1, ylab = "Approval rating (%)", + main = "presidents data") > > > > cleanEx() > nameEx("pressure") > ### * pressure > > flush(stderr()); flush(stdout()) > > ### Name: pressure > ### Title: Vapor Pressure of Mercury as a Function of Temperature > ### Aliases: pressure > ### Keywords: datasets > > ### ** Examples > > require(graphics) > plot(pressure, xlab = "Temperature (deg C)", + ylab = "Pressure (mm of Hg)", + main = "pressure data: Vapor Pressure of Mercury") > plot(pressure, xlab = "Temperature (deg C)", log = "y", + ylab = "Pressure (mm of Hg)", + main = "pressure data: Vapor Pressure of Mercury") > > > > cleanEx() > nameEx("quakes") > ### * quakes > > flush(stderr()); flush(stdout()) > > ### Name: quakes > ### Title: Locations of Earthquakes off Fiji > ### Aliases: quakes > ### Keywords: datasets > > ### ** Examples > > require(graphics) > pairs(quakes, main = "Fiji Earthquakes, N = 1000", cex.main = 1.2, pch = ".") > > > > cleanEx() > nameEx("randu") > ### * randu > > flush(stderr()); flush(stdout()) > > ### Name: randu > ### Title: Random Numbers from Congruential Generator RANDU > ### Aliases: randu > ### Keywords: datasets > > ### ** Examples > > > > > cleanEx() > nameEx("sleep") > ### * sleep > > flush(stderr()); flush(stdout()) > > ### Name: sleep > ### Title: Student's Sleep Data > ### Aliases: sleep > ### Keywords: datasets > > ### ** Examples > > require(stats) > ## Student's paired t-test > with(sleep, + t.test(extra[group == 1], + extra[group == 2], paired = TRUE)) Paired t-test data: extra[group == 1] and extra[group == 2] t = -4.0621, df = 9, p-value = 0.002833 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -2.4598858 -0.7001142 sample estimates: mean of the differences -1.58 > > ## The sleep *prolongations* > sleep1 <- with(sleep, extra[group == 2] - extra[group == 1]) > summary(sleep1) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 1.05 1.30 1.58 1.70 4.60 > stripchart(sleep1, method = "stack", xlab = "hours", + main = "Sleep prolongation (n = 10)") > boxplot(sleep1, horizontal = TRUE, add = TRUE, + at = .6, pars = list(boxwex = 0.5, staplewex = 0.25)) > > > > cleanEx() > nameEx("stackloss") > ### * stackloss > > flush(stderr()); flush(stdout()) > > ### Name: stackloss > ### Title: Brownlee's Stack Loss Plant Data > ### Aliases: stackloss stack.loss stack.x > ### Keywords: datasets > > ### ** Examples > > require(stats) > summary(lm.stack <- lm(stack.loss ~ stack.x)) Call: lm(formula = stack.loss ~ stack.x) Residuals: Min 1Q Median 3Q Max -7.2377 -1.7117 -0.4551 2.3614 5.6978 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -39.9197 11.8960 -3.356 0.00375 ** stack.xAir.Flow 0.7156 0.1349 5.307 5.8e-05 *** stack.xWater.Temp 1.2953 0.3680 3.520 0.00263 ** stack.xAcid.Conc. -0.1521 0.1563 -0.973 0.34405 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.243 on 17 degrees of freedom Multiple R-squared: 0.9136, Adjusted R-squared: 0.8983 F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09 > > > > cleanEx() > nameEx("sunspot.month") > ### * sunspot.month > > flush(stderr()); flush(stdout()) > > ### Name: sunspot.month > ### Title: Monthly Sunspot Data, from 1749 to "Present" > ### Aliases: sunspot.month > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## Compare the monthly series > plot (sunspot.month, + main="sunspot.month & sunspots [package'datasets']", col=2) > lines(sunspots) # -> faint differences where they overlap > > ## Now look at the difference : > all(tsp(sunspots) [c(1,3)] == + tsp(sunspot.month)[c(1,3)]) ## Start & Periodicity are the same [1] TRUE > n1 <- length(sunspots) > table(eq <- sunspots == sunspot.month[1:n1]) #> 132 are different ! FALSE TRUE 143 2677 > i <- which(!eq) > rug(time(eq)[i]) > s1 <- sunspots[i] ; s2 <- sunspot.month[i] > cbind(i = i, time = time(sunspots)[i], sunspots = s1, ss.month = s2, + perc.diff = round(100*2*abs(s1-s2)/(s1+s2), 1)) i time sunspots ss.month perc.diff [1,] 55 1753.500 22.2 22.0 0.9 [2,] 838 1818.750 31.7 31.6 0.3 [3,] 841 1819.000 32.5 32.8 0.9 [4,] 862 1820.750 9.0 8.9 1.1 [5,] 864 1820.917 9.7 9.1 6.4 [6,] 866 1821.083 4.3 4.2 2.4 [7,] 876 1821.917 0.0 0.2 200.0 [8,] 901 1824.000 21.6 21.7 0.5 [9,] 917 1825.333 15.4 15.5 0.6 [10,] 920 1825.583 25.4 25.7 1.2 [11,] 943 1827.500 42.9 42.3 1.4 [12,] 946 1827.750 57.2 56.1 1.9 [13,] 955 1828.500 54.3 54.2 0.2 [14,] 960 1828.917 46.6 46.9 0.6 [15,] 965 1829.333 67.5 67.4 0.1 [16,] 968 1829.583 78.3 77.6 0.9 [17,] 976 1830.250 107.1 106.3 0.7 [18,] 988 1831.250 54.6 54.5 0.2 [19,] 992 1831.583 54.9 55.0 0.2 [20,] 994 1831.750 46.2 46.3 0.2 [21,] 998 1832.083 55.5 55.6 0.2 [22,] 1003 1832.500 13.9 14.0 0.7 [23,] 1047 1836.167 98.1 98.2 0.1 [24,] 1061 1837.333 111.3 111.7 0.4 [25,] 1081 1839.000 107.6 105.6 1.9 [26,] 1087 1839.500 84.7 84.8 0.1 [27,] 1090 1839.750 90.8 90.9 0.1 [28,] 1092 1839.917 63.6 63.7 0.2 [29,] 1095 1840.167 55.5 67.8 20.0 [30,] 1102 1840.750 49.8 55.0 9.9 [31,] 1105 1841.000 24.0 24.1 0.4 [32,] 1108 1841.250 42.6 40.2 5.8 [33,] 1109 1841.333 67.4 67.5 0.1 [34,] 1113 1841.667 35.1 36.5 3.9 [35,] 1124 1842.583 26.5 26.6 0.4 [36,] 1125 1842.667 18.5 18.4 0.5 [37,] 1132 1843.250 8.8 9.5 7.7 [38,] 1145 1844.333 12.0 11.6 3.4 [39,] 1149 1844.667 6.9 7.0 1.4 [40,] 1156 1845.250 56.9 57.0 0.2 [41,] 1168 1846.250 69.2 69.3 0.1 [42,] 1185 1847.667 161.2 160.9 0.2 [43,] 1191 1848.167 108.9 108.6 0.3 [44,] 1194 1848.417 123.8 129.0 4.1 [45,] 1196 1848.583 132.5 132.6 0.1 [46,] 1200 1848.917 159.9 159.5 0.3 [47,] 1201 1849.000 156.7 157.0 0.2 [48,] 1202 1849.083 131.7 131.8 0.1 [49,] 1203 1849.167 96.5 96.2 0.3 [50,] 1206 1849.417 81.2 81.1 0.1 [51,] 1208 1849.583 61.3 67.7 9.9 [52,] 1211 1849.833 99.7 99.0 0.7 [53,] 1224 1850.917 60.0 61.0 1.7 [54,] 1235 1851.833 50.9 51.0 0.2 [55,] 1238 1852.083 67.5 66.4 1.6 [56,] 1243 1852.500 42.0 42.1 0.2 [57,] 1256 1853.583 50.4 50.5 0.2 [58,] 1258 1853.750 42.3 42.4 0.2 [59,] 1264 1854.250 26.4 26.5 0.4 [60,] 1270 1854.750 12.7 12.6 0.8 [61,] 1272 1854.917 21.4 21.6 0.9 [62,] 1282 1855.750 9.7 9.6 1.0 [63,] 1283 1855.833 4.3 4.2 2.4 [64,] 1290 1856.417 5.0 5.2 3.9 [65,] 1301 1857.333 29.2 28.5 2.4 [66,] 1333 1860.000 81.5 82.4 1.1 [67,] 1334 1860.083 88.0 88.3 0.3 [68,] 1346 1861.083 77.8 77.7 0.1 [69,] 1350 1861.417 87.8 88.1 0.3 [70,] 1366 1862.750 42.0 41.9 0.2 [71,] 1407 1866.167 24.6 24.5 0.4 [72,] 1424 1867.583 4.9 4.8 2.1 [73,] 1427 1867.833 9.3 9.6 3.2 [74,] 1429 1868.000 15.6 15.5 0.6 [75,] 1430 1868.083 15.8 15.7 0.6 [76,] 1435 1868.500 28.6 29.0 1.4 [77,] 1437 1868.667 43.8 47.2 7.5 [78,] 1438 1868.750 61.7 61.6 0.2 [79,] 1442 1869.083 59.3 59.9 1.0 [80,] 1445 1869.333 104.0 103.9 0.1 [81,] 1450 1869.750 59.4 59.3 0.2 [82,] 1451 1869.833 77.4 78.1 0.9 [83,] 1452 1869.917 104.3 104.4 0.1 [84,] 1455 1870.167 159.4 157.5 1.2 [85,] 1472 1871.583 110.0 110.1 0.1 [86,] 1476 1871.917 90.3 90.4 0.1 [87,] 1486 1872.750 103.5 102.6 0.9 [88,] 1497 1873.667 47.5 47.1 0.8 [89,] 1498 1873.750 47.4 47.1 0.6 [90,] 1514 1875.083 22.2 21.5 3.2 [91,] 1527 1876.167 31.2 30.6 1.9 [92,] 1539 1877.167 11.7 11.9 1.7 [93,] 1541 1877.333 21.2 21.6 1.9 [94,] 1542 1877.417 13.4 14.2 5.8 [95,] 1543 1877.500 5.9 6.0 1.7 [96,] 1545 1877.667 16.4 16.9 3.0 [97,] 1547 1877.833 14.5 14.2 2.1 [98,] 1548 1877.917 2.3 2.2 4.4 [99,] 1550 1878.083 6.0 6.6 9.5 [100,] 1553 1878.333 5.8 5.9 1.7 [101,] 1561 1879.000 0.8 1.0 22.2 [102,] 1571 1879.833 12.9 13.1 1.5 [103,] 1572 1879.917 7.2 7.3 1.4 [104,] 1574 1880.083 27.5 27.2 1.1 [105,] 1575 1880.167 19.5 19.3 1.0 [106,] 1576 1880.250 19.3 19.5 1.0 [107,] 1588 1881.250 51.7 51.6 0.2 [108,] 1592 1881.583 58.0 58.4 0.7 [109,] 1594 1881.750 64.0 64.4 0.6 [110,] 1598 1882.083 69.3 69.5 0.3 [111,] 1599 1882.167 67.5 66.8 1.0 [112,] 1613 1883.333 32.1 31.5 1.9 [113,] 1614 1883.417 76.5 76.3 0.3 [114,] 1623 1884.167 86.8 87.5 0.8 [115,] 1643 1885.833 33.3 30.9 7.5 [116,] 1656 1886.917 12.4 13.0 4.7 [117,] 1663 1887.500 23.3 23.4 0.4 [118,] 1683 1889.167 7.0 6.7 4.4 [119,] 1687 1889.500 9.7 9.4 3.1 [120,] 1712 1891.583 33.2 33.0 0.6 [121,] 1716 1891.917 32.3 32.5 0.6 [122,] 1723 1892.500 76.8 76.5 0.4 [123,] 1734 1893.417 88.2 89.9 1.9 [124,] 1735 1893.500 88.8 88.6 0.2 [125,] 1738 1893.750 79.7 80.0 0.4 [126,] 1774 1896.750 28.4 28.7 1.1 [127,] 1837 1902.000 5.2 5.5 5.6 [128,] 2126 1926.083 70.0 69.9 0.1 [129,] 2151 1928.167 85.4 85.5 0.1 [130,] 2153 1928.333 76.9 77.0 0.1 [131,] 2162 1929.083 64.1 62.8 2.0 [132,] 2174 1930.083 49.2 49.9 1.4 [133,] 2233 1935.000 18.9 18.6 1.6 [134,] 2315 1941.833 38.3 38.4 0.3 [135,] 2329 1943.000 12.4 12.5 0.8 [136,] 2378 1947.083 113.4 133.4 16.2 [137,] 2427 1951.167 59.9 55.9 6.9 [138,] 2498 1957.083 130.2 130.3 0.1 [139,] 2552 1961.583 55.9 55.8 0.2 [140,] 2556 1961.917 40.0 39.9 0.3 [141,] 2594 1965.083 14.2 14.3 0.7 [142,] 2790 1981.417 90.0 90.9 1.0 [143,] 2819 1983.833 33.3 33.4 0.3 > > ## How to recreate the "old" sunspot.month (R <= 3.0.3): > .sunspot.diff <- cbind( + i = c(1202L, 1256L, 1258L, 1301L, 1407L, 1429L, 1452L, 1455L, + 1663L, 2151L, 2329L, 2498L, 2594L, 2694L, 2819L), + res10 = c(1L, 1L, 1L, -1L, -1L, -1L, 1L, -1L, + 1L, 1L, 1L, 1L, 1L, 20L, 1L)) > ssm0 <- sunspot.month[1:2988] > with(as.data.frame(.sunspot.diff), ssm0[i] <<- ssm0[i] - res10/10) > sunspot.month.0 <- ts(ssm0, start = 1749, frequency = 12) > > > > cleanEx() > nameEx("sunspot.year") > ### * sunspot.year > > flush(stderr()); flush(stdout()) > > ### Name: sunspot.year > ### Title: Yearly Sunspot Data, 1700-1988 > ### Aliases: sunspot.year > ### Keywords: datasets > > ### ** Examples > > utils::str(sm <- sunspots)# the monthly version we keep unchanged Time-Series [1:2820] from 1749 to 1984: 58 62.6 70 55.7 85 83.5 94.8 66.3 75.9 75.5 ... > utils::str(sy <- sunspot.year) Time-Series [1:289] from 1700 to 1988: 5 11 16 23 36 58 29 20 10 8 ... > ## The common time interval > (t1 <- c(max(start(sm), start(sy)), 1)) # Jan 1749 [1] 1749 1 > (t2 <- c(min( end(sm)[1],end(sy)[1]), 12)) # Dec 1983 [1] 1983 12 > s.m <- window(sm, start=t1, end=t2) > s.y <- window(sy, start=t1, end=t2[1]) # {irrelevant warning} > stopifnot(length(s.y) * 12 == length(s.m), + ## The yearly series *is* close to the averages of the monthly one: + all.equal(s.y, aggregate(s.m, FUN = mean), tol = 0.0020)) > ## NOTE: Strangely, correctly weighting the number of days per month > ## (using 28.25 for February) is *not* closer than the simple mean: > ndays <- c(31, 28.25, rep(c(31,30, 31,30, 31), 2)) > all.equal(s.y, aggregate(s.m, FUN = mean)) # 0.0013 [1] "Mean relative difference: 0.001312539" > all.equal(s.y, aggregate(s.m, FUN = weighted.mean, w = ndays)) # 0.0017 [1] "Mean relative difference: 0.001692215" > > > > cleanEx() > nameEx("sunspots") > ### * sunspots > > flush(stderr()); flush(stdout()) > > ### Name: sunspots > ### Title: Monthly Sunspot Numbers, 1749-1983 > ### Aliases: sunspots > ### Keywords: datasets > > ### ** Examples > > require(graphics) > plot(sunspots, main = "sunspots data", xlab = "Year", + ylab = "Monthly sunspot numbers") > > > > cleanEx() > nameEx("swiss") > ### * swiss > > flush(stderr()); flush(stdout()) > > ### Name: swiss > ### Title: Swiss Fertility and Socioeconomic Indicators (1888) Data > ### Aliases: swiss > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > pairs(swiss, panel = panel.smooth, main = "swiss data", + col = 3 + (swiss$Catholic > 50)) > summary(lm(Fertility ~ . , data = swiss)) Call: lm(formula = Fertility ~ ., data = swiss) Residuals: Min 1Q Median 3Q Max -15.2743 -5.2617 0.5032 4.1198 15.3213 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 66.91518 10.70604 6.250 1.91e-07 *** Agriculture -0.17211 0.07030 -2.448 0.01873 * Examination -0.25801 0.25388 -1.016 0.31546 Education -0.87094 0.18303 -4.758 2.43e-05 *** Catholic 0.10412 0.03526 2.953 0.00519 ** Infant.Mortality 1.07705 0.38172 2.822 0.00734 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 7.165 on 41 degrees of freedom Multiple R-squared: 0.7067, Adjusted R-squared: 0.671 F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10 > > > > cleanEx() > nameEx("trees") > ### * trees > > flush(stderr()); flush(stdout()) > > ### Name: trees > ### Title: Diameter, Height and Volume for Black Cherry Trees > ### Aliases: trees > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > pairs(trees, panel = panel.smooth, main = "trees data") > plot(Volume ~ Girth, data = trees, log = "xy") > coplot(log(Volume) ~ log(Girth) | Height, data = trees, + panel = panel.smooth) > summary(fm1 <- lm(log(Volume) ~ log(Girth), data = trees)) Call: lm(formula = log(Volume) ~ log(Girth), data = trees) Residuals: Min 1Q Median 3Q Max -0.205999 -0.068702 0.001011 0.072585 0.247963 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.35332 0.23066 -10.20 4.18e-11 *** log(Girth) 2.19997 0.08983 24.49 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.115 on 29 degrees of freedom Multiple R-squared: 0.9539, Adjusted R-squared: 0.9523 F-statistic: 599.7 on 1 and 29 DF, p-value: < 2.2e-16 > summary(fm2 <- update(fm1, ~ . + log(Height), data = trees)) Call: lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees) Residuals: Min 1Q Median 3Q Max -0.168561 -0.048488 0.002431 0.063637 0.129223 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -6.63162 0.79979 -8.292 5.06e-09 *** log(Girth) 1.98265 0.07501 26.432 < 2e-16 *** log(Height) 1.11712 0.20444 5.464 7.81e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.08139 on 28 degrees of freedom Multiple R-squared: 0.9777, Adjusted R-squared: 0.9761 F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16 > step(fm2) Start: AIC=-152.69 log(Volume) ~ log(Girth) + log(Height) Df Sum of Sq RSS AIC 0.1855 -152.685 - log(Height) 1 0.1978 0.3832 -132.185 - log(Girth) 1 4.6275 4.8130 -53.743 Call: lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees) Coefficients: (Intercept) log(Girth) log(Height) -6.632 1.983 1.117 > ## i.e., Volume ~= c * Height * Girth^2 seems reasonable > > > > cleanEx() > nameEx("uspop") > ### * uspop > > flush(stderr()); flush(stdout()) > > ### Name: uspop > ### Title: Populations Recorded by the US Census > ### Aliases: uspop > ### Keywords: datasets > > ### ** Examples > > require(graphics) > plot(uspop, log = "y", main = "uspop data", xlab = "Year", + ylab = "U.S. Population (millions)") > > > > cleanEx() > nameEx("volcano") > ### * volcano > > flush(stderr()); flush(stdout()) > > ### Name: volcano > ### Title: Topographic Information on Auckland's Maunga Whau Volcano > ### Aliases: volcano > ### Keywords: datasets > > ### ** Examples > > require(grDevices); require(graphics) > filled.contour(volcano, color.palette = terrain.colors, asp = 1) > title(main = "volcano data: filled contour map") > > > > cleanEx() > nameEx("warpbreaks") > ### * warpbreaks > > flush(stderr()); flush(stdout()) > > ### Name: warpbreaks > ### Title: The Number of Breaks in Yarn during Weaving > ### Aliases: warpbreaks > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > summary(warpbreaks) breaks wool tension Min. :10.00 A:27 L:18 1st Qu.:18.25 B:27 M:18 Median :26.00 H:18 Mean :28.15 3rd Qu.:34.00 Max. :70.00 > opar <- par(mfrow = c(1, 2), oma = c(0, 0, 1.1, 0)) > plot(breaks ~ tension, data = warpbreaks, col = "lightgray", + varwidth = TRUE, subset = wool == "A", main = "Wool A") > plot(breaks ~ tension, data = warpbreaks, col = "lightgray", + varwidth = TRUE, subset = wool == "B", main = "Wool B") > mtext("warpbreaks data", side = 3, outer = TRUE) > par(opar) > summary(fm1 <- lm(breaks ~ wool*tension, data = warpbreaks)) Call: lm(formula = breaks ~ wool * tension, data = warpbreaks) Residuals: Min 1Q Median 3Q Max -19.5556 -6.8889 -0.6667 7.1944 25.4444 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 44.556 3.647 12.218 2.43e-16 *** woolB -16.333 5.157 -3.167 0.002677 ** tensionM -20.556 5.157 -3.986 0.000228 *** tensionH -20.000 5.157 -3.878 0.000320 *** woolB:tensionM 21.111 7.294 2.895 0.005698 ** woolB:tensionH 10.556 7.294 1.447 0.154327 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.94 on 48 degrees of freedom Multiple R-squared: 0.3778, Adjusted R-squared: 0.3129 F-statistic: 5.828 on 5 and 48 DF, p-value: 0.0002772 > anova(fm1) Analysis of Variance Table Response: breaks Df Sum Sq Mean Sq F value Pr(>F) wool 1 450.7 450.67 3.7653 0.0582130 . tension 2 2034.3 1017.13 8.4980 0.0006926 *** wool:tension 2 1002.8 501.39 4.1891 0.0210442 * Residuals 48 5745.1 119.69 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > > > graphics::par(get("par.postscript", pos = 'CheckExEnv')) > cleanEx() > nameEx("women") > ### * women > > flush(stderr()); flush(stdout()) > > ### Name: women > ### Title: Average Heights and Weights for American Women > ### Aliases: women > ### Keywords: datasets > > ### ** Examples > > require(graphics) > plot(women, xlab = "Height (in)", ylab = "Weight (lb)", + main = "women data: American women aged 30-39") > > > > cleanEx() > nameEx("zCO2") > ### * zCO2 > > flush(stderr()); flush(stdout()) > > ### Name: CO2 > ### Title: Carbon Dioxide Uptake in Grass Plants > ### Aliases: CO2 > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## Don't show: > options(show.nls.convergence=FALSE) > ## End(Don't show) > coplot(uptake ~ conc | Plant, data = CO2, show.given = FALSE, type = "b") > ## fit the data for the first plant > fm1 <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0), + data = CO2, subset = Plant == "Qn1") > summary(fm1) Formula: uptake ~ SSasymp(conc, Asym, lrc, c0) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 38.1398 0.9164 41.620 1.99e-06 *** lrc -34.2766 18.9661 -1.807 0.145 c0 -4.3806 0.2042 -21.457 2.79e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.663 on 4 degrees of freedom > ## fit each plant separately > fmlist <- list() > for (pp in levels(CO2$Plant)) { + fmlist[[pp]] <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0), + data = CO2, subset = Plant == pp) + } > ## check the coefficients by plant > print(sapply(fmlist, coef), digits = 3) Qn1 Qn2 Qn3 Qc1 Qc3 Qc2 Mn3 Mn2 Mn1 Mc2 Mc3 Asym 38.14 42.87 44.23 36.43 40.68 39.82 28.48 32.13 34.08 13.56 18.54 lrc -34.28 -29.66 -37.63 -9.90 -11.54 -51.53 -17.37 -29.04 -8.81 -1.98 -136.11 c0 -4.38 -4.67 -4.49 -4.86 -4.95 -4.46 -4.59 -4.47 -5.06 -4.56 -3.47 Mc1 Asym 21.79 lrc 2.45 c0 -5.14 > > > > ### *