R version 2.15.2 Patched (2012-11-27 r61171) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-unknown-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') > > assign(".oldSearch", search(), 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 > > require(graphics) > pairs(USArrests, panel = panel.smooth, main = "USArrests data") > > > > 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)] # 3 of them 73 312 674 -Inf Inf -Inf > 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 > if(require(survival)){ + model3 <- clogit(case ~ spontaneous+induced+strata(stratum), data = infert) + print(summary(model3)) + detach() # survival (conflicts) + } Loading required package: survival Loading required package: splines Call: coxph(formula = Surv(rep(1, 248L), case) ~ spontaneous + induced + strata(stratum), data = infert, method = "exact") n= 248, number of events= 83 coef exp(coef) se(coef) z Pr(>|z|) spontaneous 1.9859 7.2854 0.3524 5.635 1.75e-08 *** induced 1.4090 4.0919 0.3607 3.906 9.38e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 exp(coef) exp(-coef) lower .95 upper .95 spontaneous 7.285 0.1373 3.651 14.536 induced 4.092 0.2444 2.018 8.298 Rsquare= 0.193 (max possible= 0.519 ) Likelihood ratio test= 53.15 on 2 df, p=2.869e-12 Wald test = 31.84 on 2 df, p=1.221e-07 Score (logrank) test = 48.44 on 2 df, p=3.032e-11 > > > > cleanEx() detaching ‘package:splines’ > 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) > morley$Expt <- factor(morley$Expt) > morley$Run <- factor(morley$Run) > > xtabs(~ Expt + Run, data = morley) # 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 = morley, + main = "Speed of Light Data", xlab = "Experiment No.") > fm <- aov(Speed ~ Run + Expt, data = morley) > 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") > coplot(mpg ~ disp | as.factor(cyl), data = mtcars, + panel = panel.smooth, rows = 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 > > ## Not run: > ##D require(stats); require(graphics) > ##D nott <- window(nottem, end = c(1936,12)) > ##D fit <- arima(nott, order = c(1,0,0), list(order = c(2,1,0), period = 12)) > ##D nott.fore <- predict(fit, n.ahead = 36) > ##D ts.plot(nott, nott.fore$pred, nott.fore$pred+2*nott.fore$se, > ##D nott.fore$pred-2*nott.fore$se, gpars = list(col = c(1,1,4,4))) > ## End(Not run) > > > 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.)") > > > > 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 > > ## Not run: > ##D ## We could re-generate the dataset by the following R code > ##D seed <- as.double(1) > ##D RANDU <- function() { > ##D seed <<- ((2^16 + 3) * seed) %% (2^31) > ##D seed/(2^31) > ##D } > ##D for(i in 1:400) { > ##D U <- c(RANDU(), RANDU(), RANDU(), RANDU(), RANDU()) > ##D print(round(U[1:3], 6)) > ##D } > ## End(Not run) > > > > 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 > > > > 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, 1749-1997 > ### Aliases: sunspot.month > ### Keywords: datasets > > ### ** Examples > > require(stats); require(graphics) > ## Compare the monthly series > plot (sunspot.month, main = "sunspot.month [stats]", col = 2) > lines(sunspots) # "very barely" see something > > ## 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 132 2688 > i <- which(!eq) > rug(time(eq)[i]) > s1 <- sunspots[i] ; s2 <- sunspot.month[i] > cbind(i = i, sunspots = s1, ss.month = s2, + perc.diff = round(100*2*abs(s1-s2)/(s1+s2), 1)) i sunspots ss.month perc.diff [1,] 55 22.2 22.0 0.9 [2,] 838 31.7 31.6 0.3 [3,] 841 32.5 32.8 0.9 [4,] 862 9.0 8.9 1.1 [5,] 864 9.7 9.1 6.4 [6,] 866 4.3 4.2 2.4 [7,] 876 0.0 0.2 200.0 [8,] 901 21.6 21.7 0.5 [9,] 917 15.4 15.5 0.6 [10,] 920 25.4 25.7 1.2 [11,] 943 42.9 42.3 1.4 [12,] 946 57.2 56.1 1.9 [13,] 955 54.3 54.2 0.2 [14,] 960 46.6 46.9 0.6 [15,] 965 67.5 67.4 0.1 [16,] 968 78.3 77.6 0.9 [17,] 976 107.1 106.3 0.7 [18,] 988 54.6 54.5 0.2 [19,] 992 54.9 55.0 0.2 [20,] 994 46.2 46.3 0.2 [21,] 998 55.5 55.6 0.2 [22,] 1003 13.9 14.0 0.7 [23,] 1047 98.1 98.2 0.1 [24,] 1061 111.3 111.7 0.4 [25,] 1081 107.6 105.6 1.9 [26,] 1087 84.7 84.8 0.1 [27,] 1090 90.8 90.9 0.1 [28,] 1092 63.6 63.7 0.2 [29,] 1095 55.5 67.8 20.0 [30,] 1102 49.8 55.0 9.9 [31,] 1105 24.0 24.1 0.4 [32,] 1108 42.6 40.2 5.8 [33,] 1109 67.4 67.5 0.1 [34,] 1113 35.1 36.5 3.9 [35,] 1124 26.5 26.6 0.4 [36,] 1125 18.5 18.4 0.5 [37,] 1132 8.8 9.5 7.7 [38,] 1145 12.0 11.6 3.4 [39,] 1149 6.9 7.0 1.4 [40,] 1156 56.9 57.0 0.2 [41,] 1168 69.2 69.3 0.1 [42,] 1185 161.2 160.9 0.2 [43,] 1191 108.9 108.6 0.3 [44,] 1194 123.8 129.0 4.1 [45,] 1196 132.5 132.6 0.1 [46,] 1200 159.9 159.5 0.3 [47,] 1201 156.7 157.0 0.2 [48,] 1203 96.5 96.2 0.3 [49,] 1206 81.2 81.1 0.1 [50,] 1208 61.3 67.7 9.9 [51,] 1211 99.7 99.0 0.7 [52,] 1224 60.0 61.0 1.7 [53,] 1235 50.9 51.0 0.2 [54,] 1238 67.5 66.4 1.6 [55,] 1243 42.0 42.1 0.2 [56,] 1264 26.4 26.5 0.4 [57,] 1270 12.7 12.6 0.8 [58,] 1272 21.4 21.6 0.9 [59,] 1282 9.7 9.6 1.0 [60,] 1283 4.3 4.2 2.4 [61,] 1290 5.0 5.2 3.9 [62,] 1301 29.2 28.6 2.1 [63,] 1333 81.5 82.4 1.1 [64,] 1334 88.0 88.3 0.3 [65,] 1346 77.8 77.7 0.1 [66,] 1350 87.8 88.1 0.3 [67,] 1366 42.0 41.9 0.2 [68,] 1424 4.9 4.8 2.1 [69,] 1427 9.3 9.6 3.2 [70,] 1430 15.8 15.7 0.6 [71,] 1435 28.6 29.0 1.4 [72,] 1437 43.8 47.2 7.5 [73,] 1438 61.7 61.6 0.2 [74,] 1442 59.3 59.9 1.0 [75,] 1445 104.0 103.9 0.1 [76,] 1450 59.4 59.3 0.2 [77,] 1451 77.4 78.1 0.9 [78,] 1455 159.4 157.6 1.1 [79,] 1472 110.0 110.1 0.1 [80,] 1476 90.3 90.4 0.1 [81,] 1486 103.5 102.6 0.9 [82,] 1497 47.5 47.1 0.8 [83,] 1498 47.4 47.1 0.6 [84,] 1514 22.2 21.5 3.2 [85,] 1527 31.2 30.6 1.9 [86,] 1539 11.7 11.9 1.7 [87,] 1541 21.2 21.6 1.9 [88,] 1542 13.4 14.2 5.8 [89,] 1543 5.9 6.0 1.7 [90,] 1545 16.4 16.9 3.0 [91,] 1547 14.5 14.2 2.1 [92,] 1548 2.3 2.2 4.4 [93,] 1550 6.0 6.6 9.5 [94,] 1553 5.8 5.9 1.7 [95,] 1561 0.8 1.0 22.2 [96,] 1571 12.9 13.1 1.5 [97,] 1572 7.2 7.3 1.4 [98,] 1574 27.5 27.2 1.1 [99,] 1575 19.5 19.3 1.0 [100,] 1576 19.3 19.5 1.0 [101,] 1588 51.7 51.6 0.2 [102,] 1592 58.0 58.4 0.7 [103,] 1594 64.0 64.4 0.6 [104,] 1598 69.3 69.5 0.3 [105,] 1599 67.5 66.8 1.0 [106,] 1613 32.1 31.5 1.9 [107,] 1614 76.5 76.3 0.3 [108,] 1623 86.8 87.5 0.8 [109,] 1643 33.3 30.9 7.5 [110,] 1656 12.4 13.0 4.7 [111,] 1683 7.0 6.7 4.4 [112,] 1687 9.7 9.4 3.1 [113,] 1712 33.2 33.0 0.6 [114,] 1716 32.3 32.5 0.6 [115,] 1723 76.8 76.5 0.4 [116,] 1734 88.2 89.9 1.9 [117,] 1735 88.8 88.6 0.2 [118,] 1738 79.7 80.0 0.4 [119,] 1774 28.4 28.7 1.1 [120,] 1837 5.2 5.5 5.6 [121,] 2126 70.0 69.9 0.1 [122,] 2153 76.9 77.0 0.1 [123,] 2162 64.1 62.8 2.0 [124,] 2174 49.2 49.9 1.4 [125,] 2233 18.9 18.6 1.6 [126,] 2315 38.3 38.4 0.3 [127,] 2378 113.4 133.4 16.2 [128,] 2427 59.9 55.9 6.9 [129,] 2552 55.9 55.8 0.2 [130,] 2556 40.0 39.9 0.3 [131,] 2694 39.5 37.5 5.2 [132,] 2790 90.0 90.9 1.0 > > > > 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: Girth, 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 > > > > ### *