The R Project SVN
R
(root)/trunk/tests/Examples/datasets-Ex.Rout.save – Rev 81048
Rev 80819 | Rev 81944 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R Under development (unstable) (2021-10-13 r81047) -- "Unsuffered Consequences"Copyright (C) 2021 The R Foundation for Statistical ComputingPlatform: 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 localeR 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 lrc19.14258 -0.63282> fm1Nonlinear regression modelmodel: demand ~ A * (1 - exp(-exp(lrc) * Time))data: BODA lrc19.143 -0.633residual sum-of-squares: 26> # using the plinear algorithm (trace o/p differs by platform)> ## IGNORE_RDIFF_BEGIN> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,+ start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)32.946 (7.09e-01): par = (-1.0498 22.126)25.992 (1.29e-02): par = (-0.62572 19.103)25.990 (2.13e-04): par = (-0.6327 19.142)25.990 (3.55e-06): par = (-0.63282 19.143)> ## IGNORE_RDIFF_END> # 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 ‘ ’ 1Residual 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.654B 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 ‘ ’ 1Residual standard error: 0.01981 on 12 degrees of freedom> anova(fm1, fm2)Analysis of Variance TableModel 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.00478962 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.011743Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 0.005086 0.007834 0.649 0.552carb 0.876286 0.013535 64.744 3.41e-07 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.008649 on 4 degrees of freedomMultiple R-squared: 0.999, Adjusted R-squared: 0.9988F-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)> xEyeHair Brown Blue Hazel GreenBlack 68 20 15 5Brown 119 84 54 29Red 26 17 14 14Blond 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 weight0.158 0.135 0.190 0.187 0.760bitro.diameter chest.girth chest.width0.829 0.877 0.801Loadings:Factor1height 0.918arm.span 0.930forearm 0.900lower.leg 0.902weight 0.490bitro.diameter 0.413chest.girth 0.351chest.width 0.446Factor1SS loadings 4.064Proportion Var 0.508Test 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 weight0.170 0.107 0.166 0.199 0.089bitro.diameter chest.girth chest.width0.364 0.416 0.537Loadings:Factor1 Factor2height 0.865 0.287arm.span 0.927 0.181forearm 0.895 0.179lower.leg 0.859 0.252weight 0.233 0.925bitro.diameter 0.194 0.774chest.girth 0.134 0.752chest.width 0.278 0.621Factor1 Factor2SS loadings 3.335 2.617Proportion Var 0.417 0.327Cumulative Var 0.417 0.744Test 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-11Call:factanal(factors = factors, covmat = Harman23.cor)Uniquenesses:height arm.span forearm lower.leg weight0.127 0.005 0.193 0.157 0.090bitro.diameter chest.girth chest.width0.359 0.411 0.490Loadings:Factor1 Factor2 Factor3height 0.886 0.267 -0.130arm.span 0.937 0.195 0.280forearm 0.874 0.188lower.leg 0.877 0.230 -0.145weight 0.242 0.916 -0.106bitro.diameter 0.193 0.777chest.girth 0.137 0.755chest.width 0.261 0.646 0.159Factor1 Factor2 Factor3SS loadings 3.379 2.628 0.162Proportion Var 0.422 0.329 0.020Cumulative Var 0.422 0.751 0.771Test 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.00184Call:factanal(factors = factors, covmat = Harman23.cor)Uniquenesses:height arm.span forearm lower.leg weight0.137 0.005 0.191 0.116 0.138bitro.diameter chest.girth chest.width0.283 0.178 0.488Loadings:Factor1 Factor2 Factor3 Factor4height 0.879 0.277 -0.115arm.span 0.937 0.194 0.277forearm 0.875 0.191lower.leg 0.887 0.209 0.135 -0.188weight 0.246 0.882 0.111 -0.109bitro.diameter 0.187 0.822chest.girth 0.117 0.729 0.526chest.width 0.263 0.644 0.141Factor1 Factor2 Factor3 Factor4SS loadings 3.382 2.595 0.323 0.165Proportion Var 0.423 0.324 0.040 0.021Cumulative Var 0.423 0.747 0.787 0.808Test 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 PaperFormBoard0.677 0.866 0.830Flags GeneralInformation PargraphComprehension0.768 0.487 0.491SentenceCompletion WordClassification WordMeaning0.500 0.514 0.474Addition Code CountingDots0.818 0.731 0.824StraightCurvedCapitals WordRecognition NumberRecognition0.681 0.833 0.863FigureRecognition ObjectNumber NumberFigure0.775 0.812 0.778FigureWord Deduction NumericalPuzzles0.816 0.612 0.676ProblemReasoning SeriesCompletion ArithmeticProblems0.619 0.524 0.593Loadings:Factor1VisualPerception 0.569Cubes 0.366PaperFormBoard 0.412Flags 0.482GeneralInformation 0.716PargraphComprehension 0.713SentenceCompletion 0.707WordClassification 0.697WordMeaning 0.725Addition 0.426Code 0.519CountingDots 0.419StraightCurvedCapitals 0.565WordRecognition 0.408NumberRecognition 0.370FigureRecognition 0.474ObjectNumber 0.434NumberFigure 0.471FigureWord 0.429Deduction 0.623NumericalPuzzles 0.569ProblemReasoning 0.617SeriesCompletion 0.690ArithmeticProblems 0.638Factor1SS loadings 7.438Proportion Var 0.310Test 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 PaperFormBoard0.650 0.864 0.844Flags GeneralInformation PargraphComprehension0.778 0.375 0.316SentenceCompletion WordClassification WordMeaning0.319 0.503 0.258Addition Code CountingDots0.670 0.608 0.581StraightCurvedCapitals WordRecognition NumberRecognition0.567 0.832 0.850FigureRecognition ObjectNumber NumberFigure0.743 0.770 0.625FigureWord Deduction NumericalPuzzles0.792 0.629 0.579ProblemReasoning SeriesCompletion ArithmeticProblems0.634 0.539 0.553Loadings:Factor1 Factor2VisualPerception 0.506 0.306Cubes 0.304 0.209PaperFormBoard 0.297 0.260Flags 0.327 0.339GeneralInformation 0.240 0.753PargraphComprehension 0.171 0.809SentenceCompletion 0.163 0.809WordClassification 0.344 0.615WordMeaning 0.148 0.849Addition 0.563 0.115Code 0.591 0.207CountingDots 0.647StraightCurvedCapitals 0.612 0.241WordRecognition 0.315 0.263NumberRecognition 0.328 0.205FigureRecognition 0.457 0.218ObjectNumber 0.431 0.209NumberFigure 0.601 0.116FigureWord 0.399 0.222Deduction 0.379 0.477NumericalPuzzles 0.604 0.237ProblemReasoning 0.390 0.462SeriesCompletion 0.486 0.474ArithmeticProblems 0.544 0.389Factor1 Factor2SS loadings 4.573 4.548Proportion Var 0.191 0.190Cumulative Var 0.191 0.380Test 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-13Call:factanal(factors = factors, covmat = Harman74.cor)Uniquenesses:VisualPerception Cubes PaperFormBoard0.500 0.793 0.662Flags GeneralInformation PargraphComprehension0.694 0.352 0.316SentenceCompletion WordClassification WordMeaning0.300 0.502 0.256Addition Code CountingDots0.200 0.586 0.494StraightCurvedCapitals WordRecognition NumberRecognition0.569 0.838 0.848FigureRecognition ObjectNumber NumberFigure0.643 0.780 0.635FigureWord Deduction NumericalPuzzles0.788 0.590 0.580ProblemReasoning SeriesCompletion ArithmeticProblems0.597 0.498 0.500Loadings:Factor1 Factor2 Factor3VisualPerception 0.176 0.656 0.198Cubes 0.122 0.428PaperFormBoard 0.145 0.563Flags 0.239 0.487 0.107GeneralInformation 0.745 0.191 0.237PargraphComprehension 0.780 0.249 0.118SentenceCompletion 0.802 0.175 0.160WordClassification 0.571 0.327 0.256WordMeaning 0.821 0.248Addition 0.162 -0.118 0.871Code 0.198 0.219 0.572CountingDots 0.179 0.688StraightCurvedCapitals 0.190 0.381 0.499WordRecognition 0.231 0.253 0.210NumberRecognition 0.158 0.299 0.195FigureRecognition 0.108 0.557 0.186ObjectNumber 0.178 0.267 0.342NumberFigure 0.427 0.424FigureWord 0.167 0.355 0.240Deduction 0.392 0.472 0.181NumericalPuzzles 0.178 0.406 0.473ProblemReasoning 0.382 0.473 0.182SeriesCompletion 0.379 0.528 0.283ArithmeticProblems 0.377 0.226 0.554Factor1 Factor2 Factor3SS loadings 3.802 3.488 3.186Proportion Var 0.158 0.145 0.133Cumulative Var 0.158 0.304 0.436Test 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-05Call:factanal(factors = factors, covmat = Harman74.cor)Uniquenesses:VisualPerception Cubes PaperFormBoard0.438 0.780 0.644Flags GeneralInformation PargraphComprehension0.651 0.352 0.312SentenceCompletion WordClassification WordMeaning0.283 0.485 0.257Addition Code CountingDots0.240 0.551 0.435StraightCurvedCapitals WordRecognition NumberRecognition0.491 0.646 0.696FigureRecognition ObjectNumber NumberFigure0.549 0.598 0.593FigureWord Deduction NumericalPuzzles0.762 0.592 0.583ProblemReasoning SeriesCompletion ArithmeticProblems0.601 0.497 0.500Loadings:Factor1 Factor2 Factor3 Factor4VisualPerception 0.160 0.689 0.187 0.160Cubes 0.117 0.436PaperFormBoard 0.137 0.570 0.110Flags 0.233 0.527GeneralInformation 0.739 0.185 0.213 0.150PargraphComprehension 0.767 0.205 0.233SentenceCompletion 0.806 0.197 0.153WordClassification 0.569 0.339 0.242 0.132WordMeaning 0.806 0.201 0.227Addition 0.167 -0.118 0.831 0.166Code 0.180 0.120 0.512 0.374CountingDots 0.210 0.716StraightCurvedCapitals 0.188 0.438 0.525WordRecognition 0.197 0.553NumberRecognition 0.122 0.116 0.520FigureRecognition 0.408 0.525ObjectNumber 0.142 0.219 0.574NumberFigure 0.293 0.336 0.456FigureWord 0.148 0.239 0.161 0.365Deduction 0.378 0.402 0.118 0.301NumericalPuzzles 0.175 0.381 0.438 0.223ProblemReasoning 0.366 0.399 0.123 0.301SeriesCompletion 0.369 0.500 0.244 0.239ArithmeticProblems 0.370 0.158 0.496 0.304Factor1 Factor2 Factor3 Factor4SS loadings 3.647 2.872 2.657 2.290Proportion Var 0.152 0.120 0.111 0.095Cumulative Var 0.152 0.272 0.382 0.478Test 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.0224Call:factanal(factors = factors, covmat = Harman74.cor)Uniquenesses:VisualPerception Cubes PaperFormBoard0.450 0.781 0.639Flags GeneralInformation PargraphComprehension0.649 0.357 0.288SentenceCompletion WordClassification WordMeaning0.277 0.485 0.262Addition Code CountingDots0.215 0.386 0.444StraightCurvedCapitals WordRecognition NumberRecognition0.256 0.639 0.706FigureRecognition ObjectNumber NumberFigure0.550 0.614 0.596FigureWord Deduction NumericalPuzzles0.764 0.521 0.564ProblemReasoning SeriesCompletion ArithmeticProblems0.580 0.442 0.478Loadings:Factor1 Factor2 Factor3 Factor4 Factor5VisualPerception 0.161 0.658 0.136 0.182 0.199Cubes 0.113 0.435 0.107PaperFormBoard 0.135 0.562 0.107 0.116Flags 0.231 0.533GeneralInformation 0.736 0.188 0.192 0.162PargraphComprehension 0.775 0.187 0.251 0.113SentenceCompletion 0.809 0.208 0.136WordClassification 0.568 0.348 0.223 0.131WordMeaning 0.800 0.215 0.224Addition 0.175 -0.100 0.844 0.176Code 0.185 0.438 0.451 0.426CountingDots 0.222 0.690 0.101 0.140StraightCurvedCapitals 0.186 0.425 0.458 0.559WordRecognition 0.197 0.557NumberRecognition 0.121 0.130 0.508FigureRecognition 0.400 0.529ObjectNumber 0.145 0.208 0.562NumberFigure 0.306 0.325 0.452FigureWord 0.147 0.242 0.145 0.364Deduction 0.370 0.452 0.139 0.287 -0.190NumericalPuzzles 0.170 0.402 0.439 0.230ProblemReasoning 0.358 0.423 0.126 0.302SeriesCompletion 0.360 0.549 0.256 0.223 -0.107ArithmeticProblems 0.371 0.185 0.502 0.307Factor1 Factor2 Factor3 Factor4 Factor5SS loadings 3.632 2.964 2.456 2.345 0.663Proportion Var 0.151 0.124 0.102 0.098 0.028Cumulative Var 0.151 0.275 0.377 0.475 0.503Test 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 Factor5VisualPerception 0.831 -0.127 0.230Cubes 0.534PaperFormBoard 0.736 -0.290 0.136Flags 0.647 -0.104SeriesCompletion 0.555 0.126 0.127GeneralInformation 0.764PargraphComprehension 0.845 -0.140 0.140SentenceCompletion 0.872 -0.140WordClassification 0.277 0.505 0.104WordMeaning 0.846 -0.108Addition -0.334 1.012CountingDots 0.206 -0.200 0.722 0.185ArithmeticProblems 0.197 0.500 0.139WordRecognition -0.126 0.127 -0.103 0.657NumberRecognition 0.568FigureRecognition 0.399 -0.142 -0.207 0.562ObjectNumber -0.108 0.107 0.613StraightCurvedCapitals 0.542 0.247 0.618Code 0.112 0.288 0.486 0.424NumberFigure 0.255 -0.230 0.211 0.413FigureWord 0.187 0.347Deduction 0.404 0.169 0.117 -0.203NumericalPuzzles 0.393 0.368ProblemReasoning 0.381 0.188 0.169Factor1 Factor2 Factor3 Factor4 Factor5SS loadings 3.529 3.311 2.367 2.109 0.762Proportion Var 0.147 0.138 0.099 0.088 0.032Cumulative 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.7509Coefficients: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.125530dpi -0.0003369 0.0009311 -0.362 0.719173ddpi 0.4096949 0.1961971 2.088 0.042471 *---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 3.803 on 45 degrees of freedomMultiple R-squared: 0.3385, Adjusted R-squared: 0.2797F-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 2Call:ar(x = Nile)Coefficients:1 20.4081 0.1812Order 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 intercept0.4096 0.1987 919.8397s.e. 0.0974 0.0990 35.6410sigma^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 intercept0.3622 0.1678 918.3103s.e. 0.1273 0.1323 39.5037sigma^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 epsilon1469 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 fittedCall:StructTS(x = Nile, type = "trend")Variances:level slope epsilon1427 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 epsilon685.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 (function (z, notch = FALSE, width = NULL, varwidth = FALSE, :some notches went outside hinges ('box'): maybe set notch=FALSE> anova(lm(weight ~ group, data = PlantGrowth))Analysis of Variance TableResponse: weightDf 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 ‘ ’ 1Residual 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 ‘ ’ 1Residual 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.51lCl -3.2861 0.1448 -22.695 1.51e-08 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual 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)SurvivedAge No YesChild 52 57Adult 1438 654> ## Higher survival rates in females?> apply(Titanic, c(2, 4), sum)SurvivedSex No YesMale 1364 367Female 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)GenderAdmit Male FemaleAdmitted 1198 557Rejected 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 intercept0.4378 0.6281 3.2274s.e. 0.0764 0.0637 0.0131sigma^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 law0.3348 0.6672 3.3539 0.0082 -1.2224 -0.0963s.e. 0.0775 0.0612 0.0441 0.0902 0.3839 0.0166sigma^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 RapeMin. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.301st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07Median : 7.250 Median :159.0 Median :66.00 Median :20.10Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.233rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18Max. :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 <n>.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> USPersonalExpenditure1940 1945 1950 1955 1960Food and Tobacco 22.200 44.500 59.60 73.2 86.80Household Operation 10.500 15.500 29.00 36.5 46.20Medical and Health 3.530 5.760 9.71 14.0 21.10Personal Care 1.040 1.980 2.45 3.4 5.40Private Education 0.341 0.974 1.80 2.6 3.64> medpolish(log10(USPersonalExpenditure))1: 1.1263172: 1.032421Final: 1.032421Median Polish Results (Dataset: "log10(USPersonalExpenditure)")Overall: 0.9872192Row Effects:Food and Tobacco Household Operation Medical and Health Personal Care0.7880270 0.4327608 0.0000000 -0.5606543Private Education-0.7319467Column Effects:1940 1945 1950 1955 1960-0.4288933 -0.2267967 0.0000000 0.1423128 0.3058289Residuals:1940 1945 1950 1955 1960Food and Tobacco 0.000000 0.0999105 0.000000 -0.053048 -0.142555Household Operation 0.030103 -0.0028516 0.042418 0.000000 -0.061167Medical and Health -0.010551 0.0000000 0.000000 0.016596 0.031234Personal Care 0.019362 0.0968971 -0.037399 -0.037399 0.000000Private 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.10414gender: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 y1Min. : 4.0 Min. : 4.0 Min. : 4.0 Min. : 8 Min. : 4.2601st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 8 1st Qu.: 6.315Median : 9.0 Median : 9.0 Median : 9.0 Median : 8 Median : 7.580Mean : 9.0 Mean : 9.0 Mean : 9.0 Mean : 9 Mean : 7.5013rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.: 8 3rd Qu.: 8.570Max. :14.0 Max. :14.0 Max. :14.0 Max. :19 Max. :10.840y2 y3 y4Min. :3.100 Min. : 5.39 Min. : 5.2501st Qu.:6.695 1st Qu.: 6.25 1st Qu.: 6.170Median :8.140 Median : 7.11 Median : 7.040Mean :7.501 Mean : 7.50 Mean : 7.5013rd Qu.:8.950 3rd Qu.: 7.98 3rd Qu.: 8.190Max. :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 TableResponse: y1Df 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 ‘ ’ 1Analysis of Variance TableResponse: y2Df 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 ‘ ’ 1Analysis of Variance TableResponse: y3Df 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 ‘ ’ 1Analysis of Variance TableResponse: y4Df 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.0017273x1 0.5000909 0.500000 0.4997273 0.4999091> lapply(mods, function(fm) coef(summary(fm)))$lm1Estimate Std. Error t value Pr(>|t|)(Intercept) 3.0000909 1.1247468 2.667348 0.025734051x1 0.5000909 0.1179055 4.241455 0.002169629$lm2Estimate Std. Error t value Pr(>|t|)(Intercept) 3.000909 1.1253024 2.666758 0.025758941x2 0.500000 0.1179637 4.238590 0.002178816$lm3Estimate Std. Error t value Pr(>|t|)(Intercept) 3.0024545 1.1244812 2.670080 0.025619109x3 0.4997273 0.1178777 4.239372 0.002176305$lm4Estimate Std. Error t value Pr(>|t|)(Intercept) 3.0017273 1.1239211 2.670763 0.025590425x4 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 distMin. : 1.00 Min. :5.000 117 : 5 Min. : 0.501st Qu.: 9.00 1st Qu.:5.300 1028 : 4 1st Qu.: 11.32Median :18.00 Median :6.100 113 : 4 Median : 23.40Mean :14.74 Mean :6.084 112 : 3 Mean : 45.603rd Qu.:20.00 3rd Qu.:6.600 135 : 3 3rd Qu.: 47.55Max. :23.00 Max. :7.700 (Other):147 Max. :370.00NA's : 16accelMin. :0.003001st Qu.:0.04425Median :0.11300Mean :0.154223rd Qu.:0.21925Max. :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 raisesMin. :40.00 Min. :37.0 Min. :30.00 Min. :34.00 Min. :43.001st Qu.:58.75 1st Qu.:58.5 1st Qu.:45.00 1st Qu.:47.00 1st Qu.:58.25Median :65.50 Median :65.0 Median :51.50 Median :56.50 Median :63.50Mean :64.63 Mean :66.6 Mean :53.13 Mean :56.37 Mean :64.633rd Qu.:71.75 3rd Qu.:77.0 3rd Qu.:62.50 3rd Qu.:66.75 3rd Qu.:71.00Max. :85.00 Max. :90.0 Max. :83.00 Max. :75.00 Max. :88.00critical advanceMin. :49.00 Min. :25.001st Qu.:69.25 1st Qu.:35.00Median :77.50 Median :41.00Mean :74.77 Mean :42.933rd Qu.:80.00 3rd Qu.:47.75Max. :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.5990Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 10.78708 11.58926 0.931 0.361634complaints 0.61319 0.16098 3.809 0.000903 ***privileges -0.07305 0.13572 -0.538 0.595594learning 0.32033 0.16852 1.901 0.069925 .raises 0.08173 0.22148 0.369 0.715480critical 0.03838 0.14700 0.261 0.796334advance -0.21706 0.17821 -1.218 0.235577---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 7.068 on 23 degrees of freedomMultiple R-squared: 0.7326, Adjusted R-squared: 0.6628F-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.6294Coefficients: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 ‘ ’ 1Residual standard error: 6.993 on 28 degrees of freedomMultiple R-squared: 0.6813, Adjusted R-squared: 0.6699F-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.88289Coefficients: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 ‘ ’ 1Residual standard error: 0.4053 on 48 degrees of freedomMultiple R-squared: 0.7331, Adjusted R-squared: 0.7276F-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 TableModel 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 113542 47 10825 1 528.81 2.3108 0.13553 46 10634 1 190.35 0.8318 0.36664 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 (function (z, notch = FALSE, width = NULL, varwidth = FALSE, :some notches went outside hinges ('box'): maybe set notch=FALSE> anova(fm1 <- lm(weight ~ feed, data = chickwts))Analysis of Variance TableResponse: weightDf 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 959.49.5 19.69.79.8 19.9 1 1 110 1 1 2 2 110.1 1 3 1 1 110.2 2 2 2 1 2 110.3 1 1 3 2 2 3 510.4 1 1 2 3 3 4 3 310.5 1 3 7 6 4 3 1 3 1 110.6 1 4 5 9 14 6 3 1 110.7 1 2 4 9 14 16 15 7 3 1 210.8 2 5 6 14 27 10 7 1 2 110.9 2 6 14 24 27 14 10 4 111 2 6 12 15 31 37 27 17 10 611.1 3 3 12 22 26 24 26 24 7 4 111.2 3 2 7 21 30 38 29 27 20 4 1 111.3 1 5 10 24 26 39 26 24 7 211.4 3 4 9 29 56 58 26 22 10 1111.5 5 11 17 33 57 38 34 25 11 211.6 2 1 4 13 37 39 48 38 27 12 2 2 111.7 2 9 17 30 37 48 45 24 9 9 211.8 1 2 11 15 35 41 34 29 10 5 111.9 1 1 2 12 10 27 32 35 19 10 9 3 112 1 4 8 19 42 39 22 16 8 2 212.1 2 4 13 22 28 15 27 10 4 112.2 1 2 5 6 23 17 16 11 8 1 112.3 4 8 10 13 20 23 6 512.4 1 1 1 2 7 12 4 7 7 1 112.5 1 1 3 12 11 8 6 8 212.6 1 3 5 7 8 6 3 1 112.7 1 1 7 5 5 8 2 212.8 1 2 3 1 8 5 3 1 112.9 1 2 2 1 113 3 1 1 2 113.1 1 113.2 1 1 1 313.3 1 113.413.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 themnamed 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 ncontrols25-34:15 0-39g/day:23 0-9g/day:24 Min. : 0.000 Min. : 0.00035-44:15 40-79 :23 10-19 :24 1st Qu.: 0.000 1st Qu.: 1.00045-54:16 80-119 :21 20-29 :20 Median : 1.000 Median : 4.00055-64:16 120+ :21 30+ :20 Mean : 2.273 Mean : 8.80765-74:15 3rd Qu.: 4.000 3rd Qu.:10.00075+ :11 Max. :17.000 Max. :60.000> ## 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 TableModel: binomial, link: logitResponse: cbind(ncases, ncontrols)Terms added sequentially (first to last)Df Deviance Resid. Df Resid. DevNULL 87 367.95agegp 5 121.045 82 246.91tobgp 3 36.639 79 210.27alcgp 3 127.933 76 82.34tobgp:alcgp 9 5.451 67 76.89> ## 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-2.3018 -0.7234 -0.2306 0.5737 2.4290Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -4.92067 0.37344 -13.177 < 2e-16 ***agegp.L 3.81892 0.67862 5.627 1.83e-08 ***agegp.Q -1.49473 0.60671 -2.464 0.0138 *agegp.C 0.07923 0.46318 0.171 0.8642agegp^4 0.12136 0.32203 0.377 0.7063agegp^5 -0.24856 0.21153 -1.175 0.2400unclass(tobgp) 0.43955 0.09623 4.568 4.93e-06 ***unclass(alcgp) 1.06766 0.10493 10.175 < 2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 367.953 on 87 degrees of freedomResidual deviance: 91.121 on 80 degrees of freedomAIC: 222.18Number 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)euroATS 13.760300BEF 40.339900DEM 1.955830ESP 166.386000FIM 5.945730FRF 6.559570IEP 0.787564ITL 1936.270000LUF 40.339900NLG 2.203710PTE 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 PTETRUE 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"]BEF806.798> ## Convert 20 Austrian Schilling to Euro> 20 / euro["ATS"]ATS1.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.040 0.02 0.04106 163 3> faithful$better.eruptions <- ne60 / 60> te <- table(ne60)> te[te >= 4] # (too) many multiples of 5 !ne60105 108 110 112 113 120 216 230 240 245 249 250 255 260 261 262 265 270 272 2756 4 7 8 4 4 4 5 6 5 4 4 4 5 4 4 4 8 5 4276 282 2884 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.levelMin. :8.791 Min. :8.791 Min. :4.278 Min. :5.8211st Qu.:9.045 1st Qu.:9.020 1st Qu.:4.392 1st Qu.:5.948Median :9.314 Median :9.284 Median :4.510 Median :6.061Mean :9.306 Mean :9.281 Mean :4.496 Mean :6.0393rd Qu.:9.591 3rd Qu.:9.561 3rd Qu.:4.605 3rd Qu.:6.139Max. :9.794 Max. :9.775 Max. :4.710 Max. :6.200market.potentialMin. :12.971st Qu.:13.01Median :13.07Mean :13.073rd Qu.:13.12Max. :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.0267124Coefficients: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.3904price.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 ‘ ’ 1Residual standard error: 0.01473 on 34 degrees of freedomMultiple R-squared: 0.9981, Adjusted R-squared: 0.9978F-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.9362Coefficients: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 freedomResidual deviance: 279.61 on 245 degrees of freedomAIC: 285.61Number 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.6536Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -1.14924 1.41220 -0.814 0.4158age 0.03958 0.03120 1.269 0.2046parity -0.82828 0.19649 -4.215 2.49e-05 ***education6-11yrs -1.04424 0.79255 -1.318 0.1876education12+ 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 freedomResidual deviance: 257.80 on 241 degrees of freedomAIC: 271.8Number 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.45539Coefficients: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.863141GNP -3.582e-02 3.349e-02 -1.070 0.312681Unemployed -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.826212Year 1.829e+00 4.555e-01 4.016 0.003037 **---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.3049 on 9 degrees of freedomMultiple R-squared: 0.9955, Adjusted R-squared: 0.9925F-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)RunExpt 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 201 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 12 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 13 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 14 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 15 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.36321Expt 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 TableModel 1: Speed ~ ExptModel 2: Speed ~ Run + ExptRes.Df RSS Df Sum of Sq F Pr(>F)1 95 5235102 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 dratMin. :10.40 4:11 Min. : 71.1 Min. : 52.0 Min. :2.7601st Qu.:15.43 6: 7 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080Median :19.20 8:14 Median :196.3 Median :123.0 Median :3.695Mean :20.09 Mean :230.7 Mean :146.7 Mean :3.5973rd Qu.:22.80 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920Max. :33.90 Max. :472.0 Max. :335.0 Max. :4.930wt qsec vs am gear carbMin. :1.513 Min. :14.50 V:18 automatic:19 3:15 1: 71st Qu.:2.581 1st Qu.:16.89 S:14 manual :13 4:12 2:10Median :3.325 Median :17.71 5: 5 3: 3Mean :3.217 Mean :17.85 4:103rd Qu.:3.610 3rd Qu.:18.90 6: 1Max. :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.0964Coefficients: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.889490Diag4 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 freedomResidual deviance: 58.436 on 40 degrees of freedomAIC: 428.78Number of Fisher Scoring iterations: 4> plot(modUnif) # 4 plots, with warning about h_ii ~= 1Warning: 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-testdata: extra[group == 1] and extra[group == 2]t = -4.0621, df = 9, p-value = 0.002833alternative hypothesis: true mean difference is not equal to 095 percent confidence interval:-2.4598858 -0.7001142sample estimates:mean difference-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.6978Coefficients: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 ‘ ’ 1Residual standard error: 3.243 on 17 degrees of freedomMultiple R-squared: 0.9136, Adjusted R-squared: 0.8983F-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 TRUE143 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 unchangedTime-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), tolerance = 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.3213Coefficients: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.31546Education -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 ‘ ’ 1Residual standard error: 7.165 on 41 degrees of freedomMultiple R-squared: 0.7067, Adjusted R-squared: 0.671F-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.247963Coefficients: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 ‘ ’ 1Residual standard error: 0.115 on 29 degrees of freedomMultiple R-squared: 0.9539, Adjusted R-squared: 0.9523F-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.129223Coefficients: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 ‘ ’ 1Residual standard error: 0.08139 on 28 degrees of freedomMultiple R-squared: 0.9777, Adjusted R-squared: 0.9761F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16> step(fm2)Start: AIC=-152.69log(Volume) ~ log(Girth) + log(Height)Df Sum of Sq RSS AIC<none> 0.1855 -152.685- log(Height) 1 0.1978 0.3832 -132.185- log(Girth) 1 4.6275 4.8130 -53.743Call: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 tensionMin. :10.00 A:27 L:181st Qu.:18.25 B:27 M:18Median :26.00 H:18Mean :28.153rd Qu.:34.00Max. :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.4444Coefficients: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 ‘ ’ 1Residual standard error: 10.94 on 48 degrees of freedomMultiple R-squared: 0.3778, Adjusted R-squared: 0.3129F-statistic: 5.828 on 5 and 48 DF, p-value: 0.0002772> anova(fm1)Analysis of Variance TableResponse: breaksDf 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.145c0 -4.3806 0.2042 -21.457 2.79e-05 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual 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 Mc3Asym 38.14 42.87 44.23 36.43 40.68 39.82 28.48 32.13 34.08 13.56 18.54lrc -34.28 -29.66 -37.63 -9.90 -11.54 -51.53 -17.37 -29.04 -8.81 -1.98 -136.11c0 -4.38 -4.67 -4.49 -4.86 -4.95 -4.46 -4.59 -4.47 -5.06 -4.56 -3.47Mc1Asym 21.79lrc 2.45c0 -5.14>>>> ### * <FOOTER>> ###> cleanEx()> options(digits = 7L)> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")Time elapsed: 3.094 0.064 3.175 0 0> grDevices::dev.off()null device1> ###> ### Local variables: ***> ### mode: outline-minor ***> ### outline-regexp: "\\(> \\)?### [*]+" ***> ### End: ***> quit('no')