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