The R Project SVN
R
(root)/branches/R-bool/tests/Examples/datasets-Ex.Rout.save – Rev 87776
Rev 87725 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R Under development (unstable) (2025-02-20 r87772) -- "Unsuffered Consequences"Copyright (C) 2025 The R Foundation for Statistical ComputingPlatform: x86_64-pc-linux-gnuR is free software and comes with ABSOLUTELY NO WARRANTY.You are welcome to redistribute it under certain conditions.Type 'license()' or 'licence()' for distribution details.R is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R or R packages in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for 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>> ## The classic 'airline model', by full ML> (fit <- arima(log10(AirPassengers), c(0, 1, 1),+ seasonal = list(order = c(0, 1, 1), period = 12)))Call:arima(x = log10(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0,1, 1), period = 12))Coefficients:ma1 sma1-0.4018 -0.5569s.e. 0.0896 0.0731sigma^2 estimated as 0.0002543: log likelihood = 353.96, aic = -701.92> update(fit, method = "CSS")Call:arima(x = log10(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0,1, 1), period = 12), method = "CSS")Coefficients:ma1 sma1-0.3772 -0.5724s.e. 0.0883 0.0704sigma^2 estimated as 0.0002619: part log likelihood = 354.32> update(fit, x = window(log10(AirPassengers), start = 1954))Call:arima(x = window(log10(AirPassengers), start = 1954), order = c(0, 1, 1), seasonal = list(order = c(0,1, 1), period = 12))Coefficients:ma1 sma1-0.4797 -0.4460s.e. 0.1000 0.1514sigma^2 estimated as 0.0001603: log likelihood = 208.02, aic = -410.04> pred <- predict(fit, n.ahead = 24)> tl <- pred$pred - 1.96 * pred$se> tu <- pred$pred + 1.96 * pred$se> ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty = c(1, 2, 2))>> ## full ML fit is the same if the series is reversed, CSS fit is not> ap0 <- rev(log10(AirPassengers))> attributes(ap0) <- attributes(AirPassengers)> arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))Call:arima(x = ap0, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))Coefficients:ma1 sma1-0.4018 -0.5569s.e. 0.0896 0.0731sigma^2 estimated as 0.0002543: log likelihood = 353.95, aic = -701.91> arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),+ method = "CSS")Call:arima(x = ap0, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),method = "CSS")Coefficients:ma1 sma1-0.4188 -0.5224s.e. 0.0875 0.0695sigma^2 estimated as 0.0002673: part log likelihood = 353>> ## Structural Time Series> ap <- log10(AirPassengers) - 2> (fit <- StructTS(ap, type = "BSM"))Call:StructTS(x = ap, type = "BSM")Variances:level slope seas epsilon0.0001456 0.0000000 0.0002635 0.0000000> par(mfrow = c(1, 2))> plot(cbind(ap, fitted(fit)), plot.type = "single")> plot(cbind(ap, tsSmooth(fit)), plot.type = "single")>>>> graphics::par(get("par.postscript", pos = 'CheckExEnv'))> 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> ## IGNORE_RDIFF_BEGIN> # using the plinear algorithm (trace o/p differs by platform)> 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>> require(graphics)> coplot(weight ~ Time | Chick, data = ChickWeight,+ type = "b", show.given = FALSE)>>>> 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. Dev Pr(>Chi)NULL 87 367.95agegp 5 121.045 82 246.91 < 2.2e-16 ***tobgp 3 36.639 79 210.27 5.485e-08 ***alcgp 3 127.933 76 82.34 < 2.2e-16 ***tobgp:alcgp 9 5.451 67 76.89 0.7934---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1> ## 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)Coefficients: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("gait")> ### * gait>> flush(stderr()); flush(stdout())>> ### Name: gait> ### Title: Hip and Knee Angle while Walking> ### Aliases: gait> ### Keywords: datasets>> ### ** Examples>> plot(gait[, 1, ], type = "b",+ xlim = range(gait[,,1]), ylim = range(gait[,,2]),+ xlab = "Hip Angle", ylab = "Knee Angle", main = "'gait' data : Boy 1")> mtext("all other boys", col = "thistle"); grid()> matlines(gait[, -1, 1], gait[, -1, 2], type = "l", lty = 1, col = adjustcolor("thistle", 1/3))>> ## The data array, two matrices :> op <- options(width = 128) # on a wide console> aperm(gait, c(2:1, 3)), , Variable = Hip AngleTimeSubject 0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975boy1 37 36 33 29 23 18 15 12 9 6 4 6 13 22 31 38 43 44 40 35boy2 47 46 42 34 27 21 16 12 4 -1 -2 2 12 24 34 42 47 48 45 43boy3 46 44 39 34 33 27 21 15 7 1 -3 -1 10 22 34 43 51 53 50 49boy4 37 36 27 20 15 15 10 6 3 0 3 4 10 14 23 31 33 35 34 32boy5 20 18 11 8 7 5 1 0 -2 -4 -4 1 10 18 27 36 34 33 32 28boy6 57 48 44 35 31 27 21 18 14 10 8 8 11 21 32 40 49 55 56 55boy7 46 38 33 25 18 15 14 10 8 4 0 3 9 19 31 37 46 45 46 41boy8 46 46 43 40 36 30 24 17 12 7 1 3 12 24 35 44 49 49 47 43boy9 46 42 37 34 31 25 19 17 13 10 5 7 13 21 29 39 47 51 52 47boy10 35 34 29 28 19 15 13 6 0 -1 2 11 20 30 39 45 48 48 43 39boy11 38 37 33 29 26 20 14 9 3 -2 -1 5 15 24 36 44 47 45 42 40boy12 35 31 29 26 22 19 15 11 6 3 0 4 16 27 35 41 44 45 40 39boy13 34 31 27 23 19 15 13 10 6 1 0 3 8 16 27 34 40 41 39 36boy14 43 41 36 31 26 20 16 13 10 8 4 12 19 27 37 41 44 44 41 37boy15 43 37 35 28 26 21 20 15 13 7 3 8 16 25 35 43 49 50 45 46boy16 40 41 36 32 27 20 17 10 5 -2 -4 -2 7 24 38 49 57 59 54 46boy17 51 49 45 39 31 23 16 9 4 0 -5 1 10 21 33 44 51 55 56 51boy18 52 46 41 35 31 24 18 12 6 3 2 5 12 22 32 42 46 48 49 46boy19 36 33 28 22 18 13 11 7 1 -2 -3 -3 2 13 22 32 39 41 38 34boy20 35 37 33 27 22 14 9 6 4 -1 1 2 11 17 34 35 40 43 43 42boy21 46 38 30 23 17 13 7 3 0 -3 -3 1 11 20 36 48 55 57 56 50boy22 43 41 37 30 24 16 12 8 5 2 1 8 16 24 33 42 48 48 48 46boy23 55 51 47 41 35 30 26 23 20 11 8 10 19 35 41 52 57 61 63 58boy24 39 38 31 27 21 14 9 9 6 4 2 3 13 30 37 47 53 53 49 44boy25 37 34 30 27 26 19 15 10 4 0 0 4 12 22 31 40 43 43 38 38boy26 36 33 28 22 18 13 11 7 1 -2 -3 -3 2 13 22 32 39 41 38 34boy27 36 33 30 28 21 15 8 1 -5 -11 -12 -7 4 16 26 37 44 47 44 37boy28 42 40 40 34 23 15 11 7 5 6 8 12 22 33 43 51 57 58 54 46boy29 38 34 30 23 17 12 8 4 1 0 -4 -4 1 10 22 32 38 41 41 40boy30 46 47 44 37 29 23 19 14 8 3 -2 1 12 26 39 48 52 48 43 42boy31 54 48 44 37 30 27 21 18 15 11 12 20 30 41 50 56 61 59 57 58boy32 52 44 44 33 28 27 23 24 19 15 15 16 24 32 43 52 58 59 57 52boy33 32 28 26 22 19 13 8 5 -1 -6 -5 0 12 22 30 36 39 36 30 29boy34 46 41 38 31 25 20 13 7 1 -4 -3 -2 8 22 34 45 53 57 55 43boy35 46 44 40 35 31 25 19 15 10 5 3 6 14 27 40 48 53 53 50 47boy36 48 42 42 35 30 23 19 14 9 4 -1 3 9 25 37 45 52 53 52 46boy37 44 41 38 32 24 18 10 6 3 0 -2 -1 6 20 36 44 49 46 38 35boy38 55 56 51 46 41 36 30 25 21 15 8 5 9 19 31 43 52 56 59 59boy39 48 50 47 42 37 29 22 14 8 5 8 15 24 36 51 59 63 64 61 55, , Variable = Knee AngleTimeSubject 0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975boy1 10 15 18 18 15 14 13 13 15 21 26 39 52 64 70 66 57 40 22 11boy2 16 25 28 25 18 12 10 9 7 9 17 32 51 67 70 66 55 39 23 16boy3 18 27 32 32 28 23 16 11 8 11 18 33 56 72 79 77 67 46 22 14boy4 5 14 16 17 10 8 6 5 5 9 20 40 57 66 67 61 43 18 3 0boy5 2 6 6 6 5 6 6 7 8 15 24 41 59 69 71 66 44 23 5 4boy6 15 17 23 23 20 19 15 13 9 10 14 25 40 58 65 66 61 45 21 9boy7 13 16 22 17 12 9 8 10 14 16 22 36 51 64 71 68 65 40 17 3boy8 14 17 18 19 19 15 11 8 7 8 12 26 47 68 79 76 61 35 13 5boy9 15 20 23 26 25 21 16 13 8 8 10 22 40 56 65 67 61 47 26 10boy10 9 22 25 21 10 9 9 9 11 19 33 55 69 75 77 70 55 36 17 16boy11 13 24 27 23 18 13 12 15 16 19 27 43 60 70 75 69 50 25 13 6boy12 7 8 11 12 8 6 4 5 5 10 19 36 58 70 72 66 51 30 8 5boy13 9 14 16 15 15 12 13 14 14 15 20 31 45 60 71 70 62 45 26 12boy14 15 20 22 22 21 19 19 19 19 20 25 41 58 70 73 68 53 32 16 10boy15 6 11 20 18 13 9 9 8 11 12 20 38 56 67 71 65 48 25 5 8boy16 11 19 30 28 25 17 15 10 7 7 13 25 46 68 79 77 67 48 23 13boy17 24 32 35 33 29 24 18 14 12 14 20 36 55 69 75 72 60 43 26 16boy18 16 20 21 20 18 14 13 11 9 11 19 34 53 66 73 72 61 43 22 11boy19 16 20 22 21 20 20 20 19 16 16 21 31 49 69 74 71 60 41 23 15boy20 7 13 14 17 14 8 7 7 7 8 17 30 49 60 71 65 58 39 20 10boy21 21 24 25 21 16 9 4 4 5 9 19 36 56 70 80 80 68 48 27 16boy22 11 14 14 11 8 5 6 7 8 10 18 36 53 65 72 73 64 46 28 13boy23 12 17 20 20 18 12 11 9 11 10 16 30 49 68 69 69 63 44 28 8boy24 8 12 14 13 12 9 7 10 9 13 21 39 59 76 76 71 55 26 8 12boy25 11 20 22 21 21 17 16 17 15 19 28 42 59 70 73 66 49 24 6 9boy26 16 20 22 21 20 20 20 19 16 16 21 31 49 69 74 71 60 41 23 15boy27 19 26 28 28 24 18 12 8 7 8 16 33 52 68 78 80 74 57 35 18boy28 13 23 30 28 19 10 6 5 10 19 33 50 68 79 81 75 61 42 25 16boy29 11 15 19 20 18 17 15 11 9 9 12 23 40 58 71 70 56 36 19 11boy30 17 25 30 30 27 22 17 11 7 7 13 31 52 69 75 68 46 21 6 8boy31 20 20 22 16 10 10 10 12 13 17 28 47 64 72 72 63 47 23 9 19boy32 18 18 25 23 18 19 21 25 24 25 27 39 59 72 77 69 60 38 17 14boy33 9 12 16 15 14 11 8 6 5 9 18 39 56 68 71 67 54 30 9 6boy34 8 10 17 16 12 10 8 4 2 4 12 27 48 66 76 77 71 54 26 3boy35 9 18 19 19 19 15 12 11 10 12 20 36 56 72 80 75 59 34 10 9boy36 13 18 27 26 25 18 12 10 9 13 16 33 53 73 77 75 64 44 23 16boy37 19 23 26 25 21 18 14 13 14 17 24 37 56 73 82 76 59 31 7 7boy38 16 23 28 28 25 21 19 17 16 13 13 19 37 57 69 71 62 45 28 20boy39 14 25 32 34 30 20 12 7 6 11 24 45 67 79 82 76 65 46 25 15> options(op)>>>> 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)Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -1.7079 0.2677 -6.380 1.78e-10 ***spontaneous 1.1972 0.2116 5.657 1.54e-08 ***induced 0.4181 0.2056 2.033 0.042 *---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 316.17 on 247 degrees of 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)Coefficients: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>> summary(iris)Sepal.Length Sepal.Width Petal.Length Petal.WidthMin. :4.300 Min. :2.000 Min. :1.000 Min. :0.1001st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300Median :5.800 Median :3.000 Median :4.350 Median :1.300Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.1993rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500Speciessetosa :50versicolor:50virginica :50>> ## Fisher's (1936) research question: whether (compound measurements of)> ## Iris versicolor "differs twice as much from I. setosa as from I. virginica"> pairs(iris[1:4], col = iris$Species)> legend(0.5, 1, levels(iris$Species), fill = 1:3, bty = "n",+ horiz = TRUE, xjust = 0.5, yjust = 0, xpd = TRUE)>> ## equivalence of legacy array (iris3) and data.frame (iris) representation> 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]]))))> stopifnot(all.equal(ii, iris))>>>> 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 was 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>> require(stats); require(graphics)> nott <- window(nottem, end = c(1936,12))> fit <- arima(nott, order = c(1,0,0), list(order = c(2,1,0), period = 12))> nott.fore <- predict(fit, n.ahead = 36)> ts.plot(nott, nott.fore$pred, nott.fore$pred+2*nott.fore$se,+ nott.fore$pred-2*nott.fore$se, gpars = list(col = c(1,1,4,4)))>>>> 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)Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) 0.568592 0.183358 3.101 0.001929 **origin2 0.431314 0.149415 2.887 0.003893 **origin3 1.461862 0.131141 11.147 < 2e-16 ***origin4 1.788731 0.126588 14.130 < 2e-16 ***origin5 0.441011 0.144844 3.045 0.002329 **origin6 2.491735 0.121219 20.556 < 2e-16 ***origin7 1.127536 0.129032 8.738 < 2e-16 ***origin8 0.796445 0.131863 6.040 1.54e-09 ***destination2 0.873202 0.166902 5.232 1.68e-07 ***destination3 1.813992 0.153861 11.790 < 2e-16 ***destination4 2.082515 0.150887 13.802 < 2e-16 ***destination5 1.366383 0.155590 8.782 < 2e-16 ***destination6 2.816369 0.146054 19.283 < 2e-16 ***destination7 1.903918 0.147810 12.881 < 2e-16 ***destination8 1.398585 0.151658 9.222 < 2e-16 ***Diag1 1.665495 0.237383 7.016 2.28e-12 ***Diag2 0.959681 0.212122 4.524 6.06e-06 ***Diag3 0.021750 0.156530 0.139 0.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("penguins")> ### * penguins>> flush(stderr()); flush(stdout())>> ### Name: penguins> ### Title: Measurements of Penguins near Palmer Station, Antarctica> ### Aliases: penguins penguins_raw> ### Keywords: datasets>> ### ** Examples>> ## view summaries> summary(penguins)species island bill_len bill_depAdelie :152 Biscoe :168 Min. :32.10 Min. :13.10Chinstrap: 68 Dream :124 1st Qu.:39.23 1st Qu.:15.60Gentoo :124 Torgersen: 52 Median :44.45 Median :17.30Mean :43.92 Mean :17.153rd Qu.:48.50 3rd Qu.:18.70Max. :59.60 Max. :21.50NA's :2 NA's :2flipper_len body_mass sex yearMin. :172.0 Min. :2700 female:165 Min. :20071st Qu.:190.0 1st Qu.:3550 male :168 1st Qu.:2007Median :197.0 Median :4050 NA's : 11 Median :2008Mean :200.9 Mean :4202 Mean :20083rd Qu.:213.0 3rd Qu.:4750 3rd Qu.:2009Max. :231.0 Max. :6300 Max. :2009NA's :2 NA's :2> summary(penguins_raw) # not useful for character vectorsstudyName Sample Number Species RegionLength:344 Min. : 1.00 Length:344 Length:344Class :character 1st Qu.: 29.00 Class :character Class :characterMode :character Median : 58.00 Mode :character Mode :characterMean : 63.153rd Qu.: 95.25Max. :152.00Island Stage Individual ID Clutch CompletionLength:344 Length:344 Length:344 Length:344Class :character Class :character Class :character Class :characterMode :character Mode :character Mode :character Mode :characterDate Egg Culmen Length (mm) Culmen Depth (mm) Flipper Length (mm)Min. :2007-11-09 Min. :32.10 Min. :13.10 Min. :172.01st Qu.:2007-11-28 1st Qu.:39.23 1st Qu.:15.60 1st Qu.:190.0Median :2008-11-09 Median :44.45 Median :17.30 Median :197.0Mean :2008-11-27 Mean :43.92 Mean :17.15 Mean :200.93rd Qu.:2009-11-16 3rd Qu.:48.50 3rd Qu.:18.70 3rd Qu.:213.0Max. :2009-12-01 Max. :59.60 Max. :21.50 Max. :231.0NA's :2 NA's :2 NA's :2Body Mass (g) Sex Delta 15 N (o/oo) Delta 13 C (o/oo)Min. :2700 Length:344 Min. : 7.632 Min. :-27.021st Qu.:3550 Class :character 1st Qu.: 8.300 1st Qu.:-26.32Median :4050 Mode :character Median : 8.652 Median :-25.83Mean :4202 Mean : 8.733 Mean :-25.693rd Qu.:4750 3rd Qu.: 9.172 3rd Qu.:-25.06Max. :6300 Max. :10.025 Max. :-23.79NA's :2 NA's :14 NA's :13CommentsLength:344Class :characterMode :character> ## convert character vectors to factors first> dFactor <- function(dat) {+ dat[] <- lapply(dat, \(.) if (is.character(.)) as.factor(.) else .)+ dat+ }> summary(dFactor(penguins_raw))studyName Sample Number SpeciesPAL0708:110 Min. : 1.00 Adelie Penguin (Pygoscelis adeliae) :152PAL0809:114 1st Qu.: 29.00 Chinstrap penguin (Pygoscelis antarctica): 68PAL0910:120 Median : 58.00 Gentoo penguin (Pygoscelis papua) :124Mean : 63.153rd Qu.: 95.25Max. :152.00Region Island Stage Individual IDAnvers:344 Biscoe :168 Adult, 1 Egg Stage:344 N13A1 : 3Dream :124 N13A2 : 3Torgersen: 52 N18A1 : 3N18A2 : 3N21A1 : 3N21A2 : 3(Other):326Clutch Completion Date Egg Culmen Length (mm) Culmen Depth (mm)No : 36 Min. :2007-11-09 Min. :32.10 Min. :13.10Yes:308 1st Qu.:2007-11-28 1st Qu.:39.23 1st Qu.:15.60Median :2008-11-09 Median :44.45 Median :17.30Mean :2008-11-27 Mean :43.92 Mean :17.153rd Qu.:2009-11-16 3rd Qu.:48.50 3rd Qu.:18.70Max. :2009-12-01 Max. :59.60 Max. :21.50NA's :2 NA's :2Flipper Length (mm) Body Mass (g) Sex Delta 15 N (o/oo)Min. :172.0 Min. :2700 FEMALE:165 Min. : 7.6321st Qu.:190.0 1st Qu.:3550 MALE :168 1st Qu.: 8.300Median :197.0 Median :4050 NA's : 11 Median : 8.652Mean :200.9 Mean :4202 Mean : 8.7333rd Qu.:213.0 3rd Qu.:4750 3rd Qu.: 9.172Max. :231.0 Max. :6300 Max. :10.025NA's :2 NA's :2 NA's :14Delta 13 C (o/oo) CommentsMin. :-27.02 Nest never observed with full clutch.: 341st Qu.:-26.32 Not enough blood for isotopes. : 7Median :-25.83 Sexing primers did not amplify. : 4Mean :-25.69 No blood sample obtained for sexing. : 23rd Qu.:-25.06 No blood sample obtained. : 2Max. :-23.79 (Other) : 5NA's :13 NA's :290>> ## visualise distribution across factors> plot(island ~ species, data = penguins)> plot(sex ~ interaction(island, species, sep = "\n"), data = penguins)>> ## bill depth vs. length by species (color) and sex (symbol):> ## positive correlations for all species, males tend to have bigger bills> sym <- c(1, 16)> pal <- c("darkorange","purple","cyan4")> plot(bill_dep ~ bill_len, data = penguins, pch = sym[sex], col = pal[species])>> ## simplified sex dimorphism analysis for Adelie species:> ## proportion of males increases with several size measurements> adelie <- subset(penguins, species == "Adelie")> plot(sex ~ bill_len, data = adelie)> plot(sex ~ bill_dep, data = adelie)> plot(sex ~ body_mass, data = adelie)> m <- glm(sex ~ bill_len + bill_dep + body_mass, data = adelie, family = binomial)> summary(m)Call:glm(formula = sex ~ bill_len + bill_dep + body_mass, family = binomial,data = adelie)Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -85.088438 17.811711 -4.777 1.78e-06 ***bill_len 0.840401 0.237656 3.536 0.000406 ***bill_dep 1.305989 0.423911 3.081 0.002064 **body_mass 0.007790 0.001969 3.957 7.59e-05 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 202.399 on 145 degrees of freedomResidual deviance: 54.024 on 142 degrees of freedom(6 observations deleted due to missingness)AIC: 62.024Number of Fisher Scoring iterations: 8>> ## Produce the long variable names as from {palmerpenguins} pkg:> long_nms <- sub("len", "length_mm",+ sub("dep","depth_mm",+ sub("mass", "mass_g", colnames(penguins))))> ## compare long and short names:> noquote(rbind(long_nms, nms = colnames(penguins)))[,1] [,2] [,3] [,4] [,5]long_nms species island bill_length_mm bill_depth_mm flipper_length_mmnms species island bill_len bill_dep flipper_len[,6] [,7] [,8]long_nms body_mass_g sex yearnms body_mass sex year>> ## Not run:> ##D # << keeping shorter 'penguins' names in this example:> ##D colnames(penguins) <- long_nms> ## End(Not run)>>>> 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>> ## We could re-generate the dataset by the following R code> seed <- as.double(1)> RANDU <- function() {+ seed <<- ((2^16 + 3) * seed) %% (2^31)+ seed/(2^31)+ }> myrandu <- matrix(NA_real_, 400, 3, dimnames = list(NULL, c("x","y","z")))> for(i in 1:400) {+ U <- c(RANDU(), RANDU(), RANDU(), RANDU(), RANDU())+ myrandu[i,] <- round(U[1:3], 6)+ }> stopifnot(all.equal(randu, as.data.frame(myrandu), tolerance = 1e-5))>>>> 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("state")> ### * state>> flush(stderr()); flush(stdout())>> ### Encoding: UTF-8> ### Name: state> ### Title: US State Facts and Figures> ### Aliases: state state.abb state.area state.center state.division> ### state.name state.region state.x77> ### Keywords: datasets>> ### ** Examples>> (dst <- dxy <- data.frame(state.center, row.names=state.abb))x yAL -86.7509 32.5901AK -127.2500 49.2500AZ -111.6250 34.2192AR -92.2992 34.7336CA -119.7730 36.5341CO -105.5130 38.6777CT -72.3573 41.5928DE -74.9841 38.6777FL -81.6850 27.8744GA -83.3736 32.3329HI -126.2500 31.7500ID -113.9300 43.5648IL -89.3776 40.0495IN -86.0808 40.0495IA -93.3714 41.9358KS -98.1156 38.4204KY -84.7674 37.3915LA -92.2724 30.6181ME -68.9801 45.6226MD -76.6459 39.2778MA -71.5800 42.3645MI -84.6870 43.1361MN -94.6043 46.3943MS -89.8065 32.6758MO -92.5137 38.3347MT -109.3200 46.8230NE -99.5898 41.3356NV -116.8510 39.1063NH -71.3924 43.3934NJ -74.2336 39.9637NM -105.9420 34.4764NY -75.1449 43.1361NC -78.4686 35.4195ND -100.0990 47.2517OH -82.5963 40.2210OK -97.1239 35.5053OR -120.0680 43.9078PA -77.4500 40.9069RI -71.1244 41.5928SC -80.5056 33.6190SD -99.7238 44.3365TN -86.4560 35.6767TX -98.7857 31.3897UT -111.3300 39.1063VT -72.5450 44.2508VA -78.2005 37.5630WA -119.7460 47.4231WV -80.6665 38.4204WI -89.9941 44.5937WY -107.2560 43.0504> ## Alaska and Hawaii are placed just off the West Coast (for compact map drawing):> dst[c("AK", "HI"),]x yAK -127.25 49.25HI -126.25 31.75> ## state.center2 := version of state.center with "correct" coordinates for AK & HI:> ## From https://pubs.usgs.gov/gip/Elevations-Distances/elvadist.html#Geographic%20Centers> ## Alaska 63°50' N., 152°00' W., 60 miles northwest of Mount McKinley> ## Hawaii 20°15' N., 156°20' W., off Maui Island> dxy["AK",] <- c(-152. , 63.83) # or c(-152.11, 65.17)> dxy["HI",] <- c(-156.33, 20.25) # or c(-156.69, 20.89)> state.center2 <- as.list(dxy)>> plot(dxy, asp=1.2, pch=3, col=2)> text(state.center2, state.abb, cex=1/2, pos=4, offset=1/4)> i <- c("AK","HI")> do.call(arrows, c(setNames(c(dst[i,], dxy[i,]), c("x0","y0", "x1","y1")),+ col=adjustcolor(4, .7), length=1/8))> points(dst[i,], col=2)> if(FALSE) { # if(require("maps")) {+ map("state", interior = FALSE, add = TRUE)+ map("state", boundary = FALSE, lty = 2, add = TRUE)+ }>>>> cleanEx()> nameEx("sunspot.month")> ### * sunspot.month>> flush(stderr()); flush(stdout())>> ### Name: sunspot.month> ### Title: Monthly Sunspot Data, from 1749 to "Present"> ### Aliases: sunspot.m2014 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) # -> clear recalibration (to *larger* values)> at. <- seq(1750, 2030, by=10)> atyr <- at.[at. %% 50 != 0] #%> Axis(time(sunspot.month), at = atyr, side = 1,+ tck= -1/100, padj = -1.5, cex.axis = 3/4)>> ## 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.m2014[1:n1]) #> 143 are different !FALSE TRUE143 2677> i <- which(!eq)> rug(time(eq)[i])> s1 <- sunspots[i] ; s2 <- sunspot.m2014[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.month.0> .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.m2014[1:2988]> with(as.data.frame(.sunspot.diff), ssm0[i] <<- ssm0[i] - res10/10)> sunspot.month.0 <- ts(ssm0, start = 1749, frequency = 12)> stopifnot(length(sunspot.month.0) == 2988)>>>> 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 (will not be updated!)[1] 1983 12> s.m <- window(sm, start=t1, end=t2)> s.y <- window(sy, start=t1, end=t2[1])> 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: 2.05 0.135 2.209 0 0> grDevices::dev.off()null device1> ###> ### Local variables: ***> ### mode: outline-minor ***> ### outline-regexp: "\\(> \\)?### [*]+" ***> ### End: ***> quit('no')