The R Project SVN
R
(root)/branches/R-2-15-branch/tests/Examples/datasets-Ex.Rout.save – Rev 61172
Rev 61155 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R version 2.15.2 Patched (2012-11-27 r61171) -- "Trick or Treat"Copyright (C) 2012 The R Foundation for Statistical ComputingISBN 3-900051-07-0Platform: 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 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')>> 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 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> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,+ start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)32.946 : -1.0498 22.126025.992 : -0.62572 19.1031925.99 : -0.6327 19.141925.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 ‘ ’ 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 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 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>> 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> 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 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 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)] # 3 of them73 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 ncontrols25-34:15 0-39g/day:23 0-9g/day:24 Min. : 0.000 Min. : 1.0035-44:15 40-79 :23 10-19 :24 1st Qu.: 0.000 1st Qu.: 3.0045-54:16 80-119 :21 20-29 :20 Median : 1.000 Median : 6.0055-64:16 120+ :21 30+ :20 Mean : 2.273 Mean :11.0865-74:15 3rd Qu.: 4.000 3rd Qu.:14.0075+ :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 TableModel: binomial, link: logitResponse: cbind(ncases, ncontrols)Terms added sequentially (first to last)Df Deviance Resid. Df Resid. DevNULL 87 227.241agegp 5 88.128 82 139.112tobgp 3 19.085 79 120.028alcgp 3 66.054 76 53.973tobgp: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.0421Coefficients: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.73281agegp^4 0.06668 0.30776 0.217 0.82848agegp^5 -0.20288 0.19523 -1.039 0.29872unclass(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 freedomResidual deviance: 59.277 on 80 degrees of freedomAIC: 222.76Number 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> if(require(survival)){+ model3 <- clogit(case ~ spontaneous+induced+strata(stratum), data = infert)+ print(summary(model3))+ detach() # survival (conflicts)+ }Loading required package: survivalLoading required package: splinesCall:coxph(formula = Surv(rep(1, 248L), case) ~ spontaneous + induced +strata(stratum), data = infert, method = "exact")n= 248, number of events= 83coef 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 ‘ ’ 1exp(coef) exp(-coef) lower .95 upper .95spontaneous 7.285 0.1373 3.651 14.536induced 4.092 0.2444 2.018 8.298Rsquare= 0.193 (max possible= 0.519 )Likelihood ratio test= 53.15 on 2 df, p=2.869e-12Wald test = 31.84 on 2 df, p=1.221e-07Score (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.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)> morley$Expt <- factor(morley$Expt)> morley$Run <- factor(morley$Run)>> xtabs(~ Expt + Run, data = morley) # 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 = 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.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")> 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.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, 64Warning: 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-testdata: extra[group == 1] and extra[group == 2]t = -4.0621, df = 9, p-value = 0.002833alternative hypothesis: true difference in means is not equal to 095 percent confidence interval:-2.4598858 -0.7001142sample 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.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, 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 TRUE132 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.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: 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.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>> ###> cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")Time elapsed: 2.722 0.057 2.802 0 0> grDevices::dev.off()null device1> ###> ### Local variables: ***> ### mode: outline-minor ***> ### outline-regexp: "\\(> \\)?### [*]+" ***> ### End: ***> quit('no')