R version 2.14.0 alpha (2011-10-03 r57134)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> pkgname <- "datasets"
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('datasets')
>
> assign(".oldSearch", search(), pos = 'CheckExEnv')
> cleanEx()
> nameEx("AirPassengers")
> ### * AirPassengers
>
> flush(stderr()); flush(stdout())
>
> ### Name: AirPassengers
> ### Title: Monthly Airline Passenger Numbers 1949-1960
> ### Aliases: AirPassengers
> ### Keywords: datasets
>
> ### ** Examples
>
> ## Not run:
> ##D ## These are quite slow and so not run by example(AirPassengers)
> ##D
> ##D ## The classic 'airline model', by full ML
> ##D (fit <- arima(log10(AirPassengers), c(0, 1, 1),
> ##D seasonal = list(order=c(0, 1 ,1), period=12)))
> ##D update(fit, method = "CSS")
> ##D update(fit, x=window(log10(AirPassengers), start = 1954))
> ##D pred <- predict(fit, n.ahead = 24)
> ##D tl <- pred$pred - 1.96 * pred$se
> ##D tu <- pred$pred + 1.96 * pred$se
> ##D ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty = c(1,2,2))
> ##D
> ##D ## full ML fit is the same if the series is reversed, CSS fit is not
> ##D ap0 <- rev(log10(AirPassengers))
> ##D attributes(ap0) <- attributes(AirPassengers)
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order=c(0, 1 ,1), period=12))
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order=c(0, 1 ,1), period=12),
> ##D method = "CSS")
> ##D
> ##D ## Structural Time Series
> ##D ap <- log10(AirPassengers) - 2
> ##D (fit <- StructTS(ap, type= "BSM"))
> ##D par(mfrow=c(1,2))
> ##D plot(cbind(ap, fitted(fit)), plot.type = "single")
> ##D plot(cbind(ap, tsSmooth(fit)), plot.type = "single")
> ## End(Not run)
>
>
> cleanEx()
> nameEx("BOD")
> ### * BOD
>
> flush(stderr()); flush(stdout())
>
> ### Name: BOD
> ### Title: Biochemical Oxygen Demand
> ### Aliases: BOD
> ### Keywords: datasets
>
> ### ** Examples
>
> ## Don't show:
> options(show.nls.convergence=FALSE)
> old <- options(digits = 5)
> ## End Don't show
> require(stats)
> # simplest form of fitting a first-order model to these data
> fm1 <- nls(demand ~ A*(1-exp(-exp(lrc)*Time)), data = BOD,
+ start = c(A = 20, lrc = log(.35)))
> coef(fm1)
A lrc
19.14258 -0.63282
> fm1
Nonlinear regression model
model: demand ~ A * (1 - exp(-exp(lrc) * Time))
data: BOD
A lrc
19.143 -0.633
residual sum-of-squares: 26
> # using the plinear algorithm
> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,
+ start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)
32.946 : -1.0498 22.1260
25.992 : -0.62572 19.10319
25.99 : -0.6327 19.1419
25.99 : -0.63282 19.14256
> # using a self-starting model
> fm3 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD)
> summary(fm3)
Formula: demand ~ SSasympOrig(Time, A, lrc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
A 19.143 2.496 7.67 0.0016 **
lrc -0.633 0.382 -1.65 0.1733
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.55 on 4 degrees of freedom
> ## Don't show:
> options(old)
> ## End Don't show
>
>
>
> cleanEx()
> nameEx("ChickWeight")
> ### * ChickWeight
>
> flush(stderr()); flush(stdout())
>
> ### Name: ChickWeight
> ### Title: Weight versus age of chicks on different diets
> ### Aliases: ChickWeight
> ### Keywords: datasets
>
> ### ** Examples
>
>
> cleanEx()
> nameEx("DNase")
> ### * DNase
>
> flush(stderr()); flush(stdout())
>
> ### Name: DNase
> ### Title: Elisa assay of DNase
> ### Aliases: DNase
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> ## Don't show:
> options(show.nls.convergence=FALSE)
> ## End Don't show
> coplot(density ~ conc | Run, data = DNase,
+ show.given = FALSE, type = "b")
> coplot(density ~ log(conc) | Run, data = DNase,
+ show.given = FALSE, type = "b")
> ## fit a representative run
> fm1 <- nls(density ~ SSlogis( log(conc), Asym, xmid, scal ),
+ data = DNase, subset = Run == 1)
> ## compare with a four-parameter logistic
> fm2 <- nls(density ~ SSfpl( log(conc), A, B, xmid, scal ),
+ data = DNase, subset = Run == 1)
> summary(fm2)
Formula: density ~ SSfpl(log(conc), A, B, xmid, scal)
Parameters:
Estimate Std. Error t value Pr(>|t|)
A -0.007897 0.017200 -0.459 0.654
B 2.377239 0.109516 21.707 5.35e-11 ***
xmid 1.507403 0.102080 14.767 4.65e-09 ***
scal 1.062579 0.056996 18.643 3.16e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01981 on 12 degrees of freedom
> anova(fm1, fm2)
Analysis of Variance Table
Model 1: density ~ SSlogis(log(conc), Asym, xmid, scal)
Model 2: density ~ SSfpl(log(conc), A, B, xmid, scal)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 13 0.0047896
2 12 0.0047073 1 8.2314e-05 0.2098 0.6551
>
>
>
> cleanEx()
> nameEx("Formaldehyde")
> ### * Formaldehyde
>
> flush(stderr()); flush(stdout())
>
> ### Name: Formaldehyde
> ### Title: Determination of Formaldehyde
> ### Aliases: Formaldehyde
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> plot(optden ~ carb, data = Formaldehyde,
+ xlab = "Carbohydrate (ml)", ylab = "Optical Density",
+ main = "Formaldehyde data", col = 4, las = 1)
> abline(fm1 <- lm(optden ~ carb, data = Formaldehyde))
> summary(fm1)
Call:
lm(formula = optden ~ carb, data = Formaldehyde)
Residuals:
1 2 3 4 5 6
-0.006714 0.001029 0.002771 0.007143 0.007514 -0.011743
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.005086 0.007834 0.649 0.552
carb 0.876286 0.013535 64.744 3.41e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.008649 on 4 degrees of freedom
Multiple R-squared: 0.999, Adjusted R-squared: 0.9988
F-statistic: 4192 on 1 and 4 DF, p-value: 3.409e-07
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(fm1)
> par(opar)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("HairEyeColor")
> ### * HairEyeColor
>
> flush(stderr()); flush(stdout())
>
> ### Name: HairEyeColor
> ### Title: Hair and Eye Color of Statistics Students
> ### Aliases: HairEyeColor
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> ## Full mosaic
> mosaicplot(HairEyeColor)
> ## Aggregate over sex (as in Snee's original data)
> x <- apply(HairEyeColor, c(1, 2), sum)
> x
Eye
Hair Brown Blue Hazel Green
Black 68 20 15 5
Brown 119 84 54 29
Red 26 17 14 14
Blond 7 94 10 16
> mosaicplot(x, main = "Relation between hair and eye color")
>
>
>
> cleanEx()
> nameEx("Harman23.cor")
> ### * Harman23.cor
>
> flush(stderr()); flush(stdout())
>
> ### Name: Harman23.cor
> ### Title: Harman Example 2.3
> ### Aliases: Harman23.cor
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats)
> (Harman23.FA <- factanal(factors = 1, covmat = Harman23.cor))
Call:
factanal(factors = 1, covmat = Harman23.cor)
Uniquenesses:
height arm.span forearm lower.leg weight
0.158 0.135 0.190 0.187 0.760
bitro.diameter chest.girth chest.width
0.829 0.877 0.801
Loadings:
Factor1
height 0.918
arm.span 0.930
forearm 0.900
lower.leg 0.902
weight 0.490
bitro.diameter 0.413
chest.girth 0.351
chest.width 0.446
Factor1
SS loadings 4.064
Proportion Var 0.508
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 611.44 on 20 degrees of freedom.
The p-value is 1.12e-116
> for(factors in 2:4) print(update(Harman23.FA, factors = factors))
Call:
factanal(factors = factors, covmat = Harman23.cor)
Uniquenesses:
height arm.span forearm lower.leg weight
0.170 0.107 0.166 0.199 0.089
bitro.diameter chest.girth chest.width
0.364 0.416 0.537
Loadings:
Factor1 Factor2
height 0.865 0.287
arm.span 0.927 0.181
forearm 0.895 0.179
lower.leg 0.859 0.252
weight 0.233 0.925
bitro.diameter 0.194 0.774
chest.girth 0.134 0.752
chest.width 0.278 0.621
Factor1 Factor2
SS loadings 3.335 2.617
Proportion Var 0.417 0.327
Cumulative Var 0.417 0.744
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 75.74 on 13 degrees of freedom.
The p-value is 6.94e-11
Call:
factanal(factors = factors, covmat = Harman23.cor)
Uniquenesses:
height arm.span forearm lower.leg weight
0.127 0.005 0.193 0.157 0.090
bitro.diameter chest.girth chest.width
0.359 0.411 0.490
Loadings:
Factor1 Factor2 Factor3
height 0.886 0.267 -0.130
arm.span 0.937 0.195 0.280
forearm 0.874 0.188
lower.leg 0.877 0.230 -0.145
weight 0.242 0.916 -0.106
bitro.diameter 0.193 0.777
chest.girth 0.137 0.755
chest.width 0.261 0.646 0.159
Factor1 Factor2 Factor3
SS loadings 3.379 2.628 0.162
Proportion Var 0.422 0.329 0.020
Cumulative Var 0.422 0.751 0.771
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 22.81 on 7 degrees of freedom.
The p-value is 0.00184
Call:
factanal(factors = factors, covmat = Harman23.cor)
Uniquenesses:
height arm.span forearm lower.leg weight
0.137 0.005 0.191 0.116 0.138
bitro.diameter chest.girth chest.width
0.283 0.178 0.488
Loadings:
Factor1 Factor2 Factor3 Factor4
height 0.879 0.277 -0.115
arm.span 0.937 0.194 0.277
forearm 0.875 0.191
lower.leg 0.887 0.209 0.135 -0.188
weight 0.246 0.882 0.111 -0.109
bitro.diameter 0.187 0.822
chest.girth 0.117 0.729 0.526
chest.width 0.263 0.644 0.141
Factor1 Factor2 Factor3 Factor4
SS loadings 3.382 2.595 0.323 0.165
Proportion Var 0.423 0.324 0.040 0.021
Cumulative Var 0.423 0.747 0.787 0.808
Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 4.63 on 2 degrees of freedom.
The p-value is 0.0988
>
>
>
> cleanEx()
> nameEx("Harman74.cor")
> ### * Harman74.cor
>
> flush(stderr()); flush(stdout())
>
> ### Name: Harman74.cor
> ### Title: Harman Example 7.4
> ### Aliases: Harman74.cor
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats)
> (Harman74.FA <- factanal(factors = 1, covmat = Harman74.cor))
Call:
factanal(factors = 1, covmat = Harman74.cor)
Uniquenesses:
VisualPerception Cubes PaperFormBoard
0.677 0.866 0.830
Flags GeneralInformation PargraphComprehension
0.768 0.487 0.491
SentenceCompletion WordClassification WordMeaning
0.500 0.514 0.474
Addition Code CountingDots
0.818 0.731 0.824
StraightCurvedCapitals WordRecognition NumberRecognition
0.681 0.833 0.863
FigureRecognition ObjectNumber NumberFigure
0.775 0.812 0.778
FigureWord Deduction NumericalPuzzles
0.816 0.612 0.676
ProblemReasoning SeriesCompletion ArithmeticProblems
0.619 0.524 0.593
Loadings:
Factor1
VisualPerception 0.569
Cubes 0.366
PaperFormBoard 0.412
Flags 0.482
GeneralInformation 0.716
PargraphComprehension 0.713
SentenceCompletion 0.707
WordClassification 0.697
WordMeaning 0.725
Addition 0.426
Code 0.519
CountingDots 0.419
StraightCurvedCapitals 0.565
WordRecognition 0.408
NumberRecognition 0.370
FigureRecognition 0.474
ObjectNumber 0.434
NumberFigure 0.471
FigureWord 0.429
Deduction 0.623
NumericalPuzzles 0.569
ProblemReasoning 0.617
SeriesCompletion 0.690
ArithmeticProblems 0.638
Factor1
SS loadings 7.438
Proportion Var 0.310
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 622.91 on 252 degrees of freedom.
The p-value is 2.28e-33
> for(factors in 2:5) print(update(Harman74.FA, factors = factors))
Call:
factanal(factors = factors, covmat = Harman74.cor)
Uniquenesses:
VisualPerception Cubes PaperFormBoard
0.650 0.864 0.844
Flags GeneralInformation PargraphComprehension
0.778 0.375 0.316
SentenceCompletion WordClassification WordMeaning
0.319 0.503 0.258
Addition Code CountingDots
0.670 0.608 0.581
StraightCurvedCapitals WordRecognition NumberRecognition
0.567 0.832 0.850
FigureRecognition ObjectNumber NumberFigure
0.743 0.770 0.625
FigureWord Deduction NumericalPuzzles
0.792 0.629 0.579
ProblemReasoning SeriesCompletion ArithmeticProblems
0.634 0.539 0.553
Loadings:
Factor1 Factor2
VisualPerception 0.506 0.306
Cubes 0.304 0.209
PaperFormBoard 0.297 0.260
Flags 0.327 0.339
GeneralInformation 0.240 0.753
PargraphComprehension 0.171 0.809
SentenceCompletion 0.163 0.809
WordClassification 0.344 0.615
WordMeaning 0.148 0.849
Addition 0.563 0.115
Code 0.591 0.207
CountingDots 0.647
StraightCurvedCapitals 0.612 0.241
WordRecognition 0.315 0.263
NumberRecognition 0.328 0.205
FigureRecognition 0.457 0.218
ObjectNumber 0.431 0.209
NumberFigure 0.601 0.116
FigureWord 0.399 0.222
Deduction 0.379 0.477
NumericalPuzzles 0.604 0.237
ProblemReasoning 0.390 0.462
SeriesCompletion 0.486 0.474
ArithmeticProblems 0.544 0.389
Factor1 Factor2
SS loadings 4.573 4.548
Proportion Var 0.191 0.190
Cumulative Var 0.191 0.380
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 420.24 on 229 degrees of freedom.
The p-value is 2.01e-13
Call:
factanal(factors = factors, covmat = Harman74.cor)
Uniquenesses:
VisualPerception Cubes PaperFormBoard
0.500 0.793 0.662
Flags GeneralInformation PargraphComprehension
0.694 0.352 0.316
SentenceCompletion WordClassification WordMeaning
0.300 0.502 0.256
Addition Code CountingDots
0.200 0.586 0.494
StraightCurvedCapitals WordRecognition NumberRecognition
0.569 0.838 0.848
FigureRecognition ObjectNumber NumberFigure
0.643 0.780 0.635
FigureWord Deduction NumericalPuzzles
0.788 0.590 0.580
ProblemReasoning SeriesCompletion ArithmeticProblems
0.597 0.498 0.500
Loadings:
Factor1 Factor2 Factor3
VisualPerception 0.176 0.656 0.198
Cubes 0.122 0.428
PaperFormBoard 0.145 0.563
Flags 0.239 0.487 0.107
GeneralInformation 0.745 0.191 0.237
PargraphComprehension 0.780 0.249 0.118
SentenceCompletion 0.802 0.175 0.160
WordClassification 0.571 0.327 0.256
WordMeaning 0.821 0.248
Addition 0.162 -0.118 0.871
Code 0.198 0.219 0.572
CountingDots 0.179 0.688
StraightCurvedCapitals 0.190 0.381 0.499
WordRecognition 0.231 0.253 0.210
NumberRecognition 0.158 0.299 0.195
FigureRecognition 0.108 0.557 0.186
ObjectNumber 0.178 0.267 0.342
NumberFigure 0.427 0.424
FigureWord 0.167 0.355 0.240
Deduction 0.392 0.472 0.181
NumericalPuzzles 0.178 0.406 0.473
ProblemReasoning 0.382 0.473 0.182
SeriesCompletion 0.379 0.528 0.283
ArithmeticProblems 0.377 0.226 0.554
Factor1 Factor2 Factor3
SS loadings 3.802 3.488 3.186
Proportion Var 0.158 0.145 0.133
Cumulative Var 0.158 0.304 0.436
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 295.59 on 207 degrees of freedom.
The p-value is 5.12e-05
Call:
factanal(factors = factors, covmat = Harman74.cor)
Uniquenesses:
VisualPerception Cubes PaperFormBoard
0.438 0.780 0.644
Flags GeneralInformation PargraphComprehension
0.651 0.352 0.312
SentenceCompletion WordClassification WordMeaning
0.283 0.485 0.257
Addition Code CountingDots
0.240 0.551 0.435
StraightCurvedCapitals WordRecognition NumberRecognition
0.491 0.646 0.696
FigureRecognition ObjectNumber NumberFigure
0.549 0.598 0.593
FigureWord Deduction NumericalPuzzles
0.762 0.592 0.583
ProblemReasoning SeriesCompletion ArithmeticProblems
0.601 0.497 0.500
Loadings:
Factor1 Factor2 Factor3 Factor4
VisualPerception 0.160 0.689 0.187 0.160
Cubes 0.117 0.436
PaperFormBoard 0.137 0.570 0.110
Flags 0.233 0.527
GeneralInformation 0.739 0.185 0.213 0.150
PargraphComprehension 0.767 0.205 0.233
SentenceCompletion 0.806 0.197 0.153
WordClassification 0.569 0.339 0.242 0.132
WordMeaning 0.806 0.201 0.227
Addition 0.167 -0.118 0.831 0.166
Code 0.180 0.120 0.512 0.374
CountingDots 0.210 0.716
StraightCurvedCapitals 0.188 0.438 0.525
WordRecognition 0.197 0.553
NumberRecognition 0.122 0.116 0.520
FigureRecognition 0.408 0.525
ObjectNumber 0.142 0.219 0.574
NumberFigure 0.293 0.336 0.456
FigureWord 0.148 0.239 0.161 0.365
Deduction 0.378 0.402 0.118 0.301
NumericalPuzzles 0.175 0.381 0.438 0.223
ProblemReasoning 0.366 0.399 0.123 0.301
SeriesCompletion 0.369 0.500 0.244 0.239
ArithmeticProblems 0.370 0.158 0.496 0.304
Factor1 Factor2 Factor3 Factor4
SS loadings 3.647 2.872 2.657 2.290
Proportion Var 0.152 0.120 0.111 0.095
Cumulative Var 0.152 0.272 0.382 0.478
Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 226.68 on 186 degrees of freedom.
The p-value is 0.0224
Call:
factanal(factors = factors, covmat = Harman74.cor)
Uniquenesses:
VisualPerception Cubes PaperFormBoard
0.450 0.781 0.639
Flags GeneralInformation PargraphComprehension
0.649 0.357 0.288
SentenceCompletion WordClassification WordMeaning
0.277 0.485 0.262
Addition Code CountingDots
0.215 0.386 0.444
StraightCurvedCapitals WordRecognition NumberRecognition
0.256 0.639 0.706
FigureRecognition ObjectNumber NumberFigure
0.550 0.614 0.596
FigureWord Deduction NumericalPuzzles
0.764 0.521 0.564
ProblemReasoning SeriesCompletion ArithmeticProblems
0.580 0.442 0.478
Loadings:
Factor1 Factor2 Factor3 Factor4 Factor5
VisualPerception 0.161 0.658 0.136 0.182 0.199
Cubes 0.113 0.435 0.107
PaperFormBoard 0.135 0.562 0.107 0.116
Flags 0.231 0.533
GeneralInformation 0.736 0.188 0.192 0.162
PargraphComprehension 0.775 0.187 0.251 0.113
SentenceCompletion 0.809 0.208 0.136
WordClassification 0.568 0.348 0.223 0.131
WordMeaning 0.800 0.215 0.224
Addition 0.175 -0.100 0.844 0.176
Code 0.185 0.438 0.451 0.426
CountingDots 0.222 0.690 0.101 0.140
StraightCurvedCapitals 0.186 0.425 0.458 0.559
WordRecognition 0.197 0.557
NumberRecognition 0.121 0.130 0.508
FigureRecognition 0.400 0.529
ObjectNumber 0.145 0.208 0.562
NumberFigure 0.306 0.325 0.452
FigureWord 0.147 0.242 0.145 0.364
Deduction 0.370 0.452 0.139 0.287 -0.190
NumericalPuzzles 0.170 0.402 0.439 0.230
ProblemReasoning 0.358 0.423 0.126 0.302
SeriesCompletion 0.360 0.549 0.256 0.223 -0.107
ArithmeticProblems 0.371 0.185 0.502 0.307
Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings 3.632 2.964 2.456 2.345 0.663
Proportion Var 0.151 0.124 0.102 0.098 0.028
Cumulative Var 0.151 0.275 0.377 0.475 0.503
Test of the hypothesis that 5 factors are sufficient.
The chi square statistic is 186.82 on 166 degrees of freedom.
The p-value is 0.128
> Harman74.FA <- factanal(factors = 5, covmat = Harman74.cor,
+ rotation="promax")
> print(Harman74.FA$loadings, sort = TRUE)
Loadings:
Factor1 Factor2 Factor3 Factor4 Factor5
VisualPerception 0.831 -0.127 0.230
Cubes 0.534
PaperFormBoard 0.736 -0.290 0.136
Flags 0.647 -0.104
SeriesCompletion 0.555 0.126 0.127
GeneralInformation 0.764
PargraphComprehension 0.845 -0.140 0.140
SentenceCompletion 0.872 -0.140
WordClassification 0.277 0.505 0.104
WordMeaning 0.846 -0.108
Addition -0.334 1.012
CountingDots 0.206 -0.200 0.722 0.185
ArithmeticProblems 0.197 0.500 0.139
WordRecognition -0.126 0.127 -0.103 0.657
NumberRecognition 0.568
FigureRecognition 0.399 -0.142 -0.207 0.562
ObjectNumber -0.108 0.107 0.613
StraightCurvedCapitals 0.542 0.247 0.618
Code 0.112 0.288 0.486 0.424
NumberFigure 0.255 -0.230 0.211 0.413
FigureWord 0.187 0.347
Deduction 0.404 0.169 0.117 -0.203
NumericalPuzzles 0.393 0.368
ProblemReasoning 0.381 0.188 0.169
Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings 3.529 3.311 2.367 2.109 0.762
Proportion Var 0.147 0.138 0.099 0.088 0.032
Cumulative Var 0.147 0.285 0.384 0.471 0.503
>
>
>
> cleanEx()
> nameEx("InsectSprays")
> ### * InsectSprays
>
> flush(stderr()); flush(stdout())
>
> ### Name: InsectSprays
> ### Title: Effectiveness of Insect Sprays
> ### Aliases: InsectSprays
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> boxplot(count ~ spray, data = InsectSprays,
+ xlab = "Type of spray", ylab = "Insect count",
+ main = "InsectSprays data", varwidth = TRUE, col = "lightgray")
> fm1 <- aov(count ~ spray, data = InsectSprays)
> summary(fm1)
Df Sum Sq Mean Sq F value Pr(>F)
spray 5 2669 533.8 34.7 <2e-16 ***
Residuals 66 1015 15.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(fm1)
> fm2 <- aov(sqrt(count) ~ spray, data = InsectSprays)
> summary(fm2)
Df Sum Sq Mean Sq F value Pr(>F)
spray 5 88.44 17.688 44.8 <2e-16 ***
Residuals 66 26.06 0.395
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> plot(fm2)
> par(opar)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("JohnsonJohnson")
> ### * JohnsonJohnson
>
> flush(stderr()); flush(stdout())
>
> ### Name: JohnsonJohnson
> ### Title: Quarterly Earnings per Johnson & Johnson Share
> ### Aliases: JohnsonJohnson
> ### Keywords: datasets
>
> ### ** Examples
>
>
> cleanEx()
> nameEx("LifeCycleSavings")
> ### * LifeCycleSavings
>
> flush(stderr()); flush(stdout())
>
> ### Name: LifeCycleSavings
> ### Title: Intercountry Life-Cycle Savings Data
> ### Aliases: LifeCycleSavings
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> pairs(LifeCycleSavings, panel = panel.smooth,
+ main = "LifeCycleSavings data")
> fm1 <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
> summary(fm1)
Call:
lm(formula = sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
Residuals:
Min 1Q Median 3Q Max
-8.2422 -2.6857 -0.2488 2.4280 9.7509
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.5660865 7.3545161 3.884 0.000334 ***
pop15 -0.4611931 0.1446422 -3.189 0.002603 **
pop75 -1.6914977 1.0835989 -1.561 0.125530
dpi -0.0003369 0.0009311 -0.362 0.719173
ddpi 0.4096949 0.1961971 2.088 0.042471 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.803 on 45 degrees of freedom
Multiple R-squared: 0.3385, Adjusted R-squared: 0.2797
F-statistic: 5.756 on 4 and 45 DF, p-value: 0.0007904
>
>
>
> cleanEx()
> nameEx("Loblolly")
> ### * Loblolly
>
> flush(stderr()); flush(stdout())
>
> ### Name: Loblolly
> ### Title: Growth of Loblolly pine trees
> ### Aliases: Loblolly
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> plot(height ~ age, data = Loblolly, subset = Seed == 329,
+ xlab = "Tree age (yr)", las = 1,
+ ylab = "Tree height (ft)",
+ main = "Loblolly data and fitted curve (Seed 329 only)")
> fm1 <- nls(height ~ SSasymp(age, Asym, R0, lrc),
+ data = Loblolly, subset = Seed == 329)
> age <- seq(0, 30, length.out = 101)
> lines(age, predict(fm1, list(age = age)))
>
>
>
> cleanEx()
> nameEx("Nile")
> ### * Nile
>
> flush(stderr()); flush(stdout())
>
> ### Name: Nile
> ### Title: Flow of the River Nile
> ### Aliases: Nile
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> par(mfrow = c(2,2))
> plot(Nile)
> acf(Nile)
> pacf(Nile)
> ar(Nile) # selects order 2
Call:
ar(x = Nile)
Coefficients:
1 2
0.4081 0.1812
Order selected 2 sigma^2 estimated as 21247
> cpgram(ar(Nile)$resid)
> par(mfrow = c(1,1))
> arima(Nile, c(2, 0, 0))
Call:
arima(x = Nile, order = c(2, 0, 0))
Coefficients:
ar1 ar2 intercept
0.4096 0.1987 919.8397
s.e. 0.0974 0.0990 35.6410
sigma^2 estimated as 20291: log likelihood = -637.98, aic = 1283.96
>
> ## Now consider missing values, following Durbin & Koopman
> NileNA <- Nile
> NileNA[c(21:40, 61:80)] <- NA
> arima(NileNA, c(2, 0, 0))
Call:
arima(x = NileNA, order = c(2, 0, 0))
Coefficients:
ar1 ar2 intercept
0.3622 0.1678 918.3103
s.e. 0.1273 0.1323 39.5037
sigma^2 estimated as 23676: log likelihood = -387.7, aic = 783.41
> plot(NileNA)
> pred <-
+ predict(arima(window(NileNA, 1871, 1890), c(2,0,0)), n.ahead = 20)
> lines(pred$pred, lty = 3, col = "red")
> lines(pred$pred + 2*pred$se, lty=2, col="blue")
> lines(pred$pred - 2*pred$se, lty=2, col="blue")
> pred <-
+ predict(arima(window(NileNA, 1871, 1930), c(2,0,0)), n.ahead = 20)
> lines(pred$pred, lty = 3, col = "red")
> lines(pred$pred + 2*pred$se, lty=2, col="blue")
> lines(pred$pred - 2*pred$se, lty=2, col="blue")
>
> ## Structural time series models
> par(mfrow = c(3, 1))
> plot(Nile)
> ## local level model
> (fit <- StructTS(Nile, type = "level"))
Call:
StructTS(x = Nile, type = "level")
Variances:
level epsilon
1469 15099
> lines(fitted(fit), lty = 2) # contemporaneous smoothing
> lines(tsSmooth(fit), lty = 2, col = 4) # fixed-interval smoothing
> plot(residuals(fit)); abline(h = 0, lty = 3)
> ## local trend model
> (fit2 <- StructTS(Nile, type = "trend")) ## constant trend fitted
Call:
StructTS(x = Nile, type = "trend")
Variances:
level slope epsilon
1427 0 15047
> pred <- predict(fit, n.ahead = 30)
> ## with 50% confidence interval
> ts.plot(Nile, pred$pred,
+ pred$pred + 0.67*pred$se, pred$pred -0.67*pred$se)
>
> ## Now consider missing values
> plot(NileNA)
> (fit3 <- StructTS(NileNA, type = "level"))
Call:
StructTS(x = NileNA, type = "level")
Variances:
level epsilon
685.8 17899.8
> lines(fitted(fit3), lty = 2)
> lines(tsSmooth(fit3), lty = 3)
> plot(residuals(fit3)); abline(h = 0, lty = 3)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("Orange")
> ### * Orange
>
> flush(stderr()); flush(stdout())
>
> ### Name: Orange
> ### Title: Growth of Orange Trees
> ### Aliases: Orange
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> coplot(circumference ~ age | Tree, data = Orange, show.given = FALSE)
> fm1 <- nls(circumference ~ SSlogis(age, Asym, xmid, scal),
+ data = Orange, subset = Tree == 3)
> plot(circumference ~ age, data = Orange, subset = Tree == 3,
+ xlab = "Tree age (days since 1968/12/31)",
+ ylab = "Tree circumference (mm)", las = 1,
+ main = "Orange tree data and fitted model (Tree 3 only)")
> age <- seq(0, 1600, length.out = 101)
> lines(age, predict(fm1, list(age = age)))
>
>
>
> cleanEx()
> nameEx("OrchardSprays")
> ### * OrchardSprays
>
> flush(stderr()); flush(stdout())
>
> ### Name: OrchardSprays
> ### Title: Potency of Orchard Sprays
> ### Aliases: OrchardSprays
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> pairs(OrchardSprays, main = "OrchardSprays data")
>
>
>
> cleanEx()
> nameEx("PlantGrowth")
> ### * PlantGrowth
>
> flush(stderr()); flush(stdout())
>
> ### Name: PlantGrowth
> ### Title: Results from an Experiment on Plant Growth
> ### Aliases: PlantGrowth
> ### Keywords: datasets
>
> ### ** Examples
>
> ## One factor ANOVA example from Dobson's book, cf. Table 7.4:
> require(stats); require(graphics)
> boxplot(weight ~ group, data = PlantGrowth, main = "PlantGrowth data",
+ ylab = "Dried weight of plants", col = "lightgray",
+ notch = TRUE, varwidth = TRUE)
Warning in bxp(list(stats = c(4.17, 4.53, 5.155, 5.33, 6.11, 3.59, 4.17, :
some notches went outside hinges ('box'): maybe set notch=FALSE
> anova(lm(weight ~ group, data = PlantGrowth))
Analysis of Variance Table
Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.7663 1.8832 4.8461 0.01591 *
Residuals 27 10.4921 0.3886
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
> cleanEx()
> nameEx("Puromycin")
> ### * Puromycin
>
> flush(stderr()); flush(stdout())
>
> ### Name: Puromycin
> ### Title: Reaction Velocity of an Enzymatic Reaction
> ### Aliases: Puromycin
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> ## Don't show:
> options(show.nls.convergence=FALSE)
> ## End Don't show
> plot(rate ~ conc, data = Puromycin, las = 1,
+ xlab = "Substrate concentration (ppm)",
+ ylab = "Reaction velocity (counts/min/min)",
+ pch = as.integer(Puromycin$state),
+ col = as.integer(Puromycin$state),
+ main = "Puromycin data and fitted Michaelis-Menten curves")
> ## simplest form of fitting the Michaelis-Menten model to these data
> fm1 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
+ subset = state == "treated",
+ start = c(Vm = 200, K = 0.05))
> fm2 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
+ subset = state == "untreated",
+ start = c(Vm = 160, K = 0.05))
> summary(fm1)
Formula: rate ~ Vm * conc/(K + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Vm 2.127e+02 6.947e+00 30.615 3.24e-11 ***
K 6.412e-02 8.281e-03 7.743 1.57e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.93 on 10 degrees of freedom
> summary(fm2)
Formula: rate ~ Vm * conc/(K + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Vm 1.603e+02 6.480e+00 24.734 1.38e-09 ***
K 4.771e-02 7.782e-03 6.131 0.000173 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 9.773 on 9 degrees of freedom
> ## add fitted lines to the plot
> conc <- seq(0, 1.2, length.out = 101)
> lines(conc, predict(fm1, list(conc = conc)), lty = 1, col = 1)
> lines(conc, predict(fm2, list(conc = conc)), lty = 2, col = 2)
> legend(0.8, 120, levels(Puromycin$state),
+ col = 1:2, lty = 1:2, pch = 1:2)
>
> ## using partial linearity
> fm3 <- nls(rate ~ conc/(K + conc), data = Puromycin,
+ subset = state == "treated", start = c(K = 0.05),
+ algorithm = "plinear")
>
>
>
> cleanEx()
> nameEx("Theoph")
> ### * Theoph
>
> flush(stderr()); flush(stdout())
>
> ### Name: Theoph
> ### Title: Pharmacokinetics of Theophylline
> ### Aliases: Theoph
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> ## Don't show:
> options(show.nls.convergence=FALSE)
> ## End Don't show
> coplot(conc ~ Time | Subject, data = Theoph, show.given = FALSE)
> Theoph.4 <- subset(Theoph, Subject == 4)
> fm1 <- nls(conc ~ SSfol(Dose, Time, lKe, lKa, lCl),
+ data = Theoph.4)
> summary(fm1)
Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl)
Parameters:
Estimate Std. Error t value Pr(>|t|)
lKe -2.4365 0.2257 -10.797 4.77e-06 ***
lKa 0.1583 0.2297 0.689 0.51
lCl -3.2861 0.1448 -22.695 1.51e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8465 on 8 degrees of freedom
> plot(conc ~ Time, data = Theoph.4,
+ xlab = "Time since drug administration (hr)",
+ ylab = "Theophylline concentration (mg/L)",
+ main = "Observed concentrations and fitted model",
+ sub = "Theophylline data - Subject 4 only",
+ las = 1, col = 4)
> xvals <- seq(0, par("usr")[2], length.out = 55)
> lines(xvals, predict(fm1, newdata = list(Time = xvals)),
+ col = 4)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("Titanic")
> ### * Titanic
>
> flush(stderr()); flush(stdout())
>
> ### Name: Titanic
> ### Title: Survival of passengers on the Titanic
> ### Aliases: Titanic
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> mosaicplot(Titanic, main = "Survival on the Titanic")
> ## Higher survival rates in children?
> apply(Titanic, c(3, 4), sum)
Survived
Age No Yes
Child 52 57
Adult 1438 654
> ## Higher survival rates in females?
> apply(Titanic, c(2, 4), sum)
Survived
Sex No Yes
Male 1364 367
Female 126 344
> ## Use loglm() in package 'MASS' for further analysis ...
>
>
>
> cleanEx()
> nameEx("ToothGrowth")
> ### * ToothGrowth
>
> flush(stderr()); flush(stdout())
>
> ### Name: ToothGrowth
> ### Title: The Effect of Vitamin C on Tooth Growth in Guinea Pigs
> ### Aliases: ToothGrowth
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
+ xlab = "ToothGrowth data: length vs dose, given type of supplement")
>
>
>
> cleanEx()
> nameEx("UCBAdmissions")
> ### * UCBAdmissions
>
> flush(stderr()); flush(stdout())
>
> ### Name: UCBAdmissions
> ### Title: Student Admissions at UC Berkeley
> ### Aliases: UCBAdmissions
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> ## Data aggregated over departments
> apply(UCBAdmissions, c(1, 2), sum)
Gender
Admit Male Female
Admitted 1198 557
Rejected 1493 1278
> mosaicplot(apply(UCBAdmissions, c(1, 2), sum),
+ main = "Student admissions at UC Berkeley")
> ## Data for individual departments
> opar <- par(mfrow = c(2, 3), oma = c(0, 0, 2, 0))
> for(i in 1:6)
+ mosaicplot(UCBAdmissions[,,i],
+ xlab = "Admit", ylab = "Sex",
+ main = paste("Department", LETTERS[i]))
> mtext(expression(bold("Student admissions at UC Berkeley")),
+ outer = TRUE, cex = 1.5)
> par(opar)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("UKDriverDeaths")
> ### * UKDriverDeaths
>
> flush(stderr()); flush(stdout())
>
> ### Name: UKDriverDeaths
> ### Title: Road Casualties in Great Britain 1969-84
> ### Aliases: UKDriverDeaths Seatbelts
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> ## work with pre-seatbelt period to identify a model, use logs
> work <- window(log10(UKDriverDeaths), end = 1982+11/12)
> par(mfrow = c(3,1))
> plot(work); acf(work); pacf(work)
> par(mfrow = c(1,1))
> (fit <- arima(work, c(1,0,0), seasonal = list(order= c(1,0,0))))
Call:
arima(x = work, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0)))
Coefficients:
ar1 sar1 intercept
0.4378 0.6281 3.2274
s.e. 0.0764 0.0637 0.0131
sigma^2 estimated as 0.00157: log likelihood = 300.85, aic = -593.7
> z <- predict(fit, n.ahead = 24)
> ts.plot(log10(UKDriverDeaths), z$pred, z$pred+2*z$se, z$pred-2*z$se,
+ lty = c(1,3,2,2), col = c("black", "red", "blue", "blue"))
>
> ## now see the effect of the explanatory variables
> X <- Seatbelts[, c("kms", "PetrolPrice", "law")]
> X[, 1] <- log10(X[, 1]) - 4
> arima(log10(Seatbelts[, "drivers"]), c(1,0,0),
+ seasonal = list(order= c(1,0,0)), xreg = X)
Call:
arima(x = log10(Seatbelts[, "drivers"]), order = c(1, 0, 0), seasonal = list(order = c(1,
0, 0)), xreg = X)
Coefficients:
ar1 sar1 intercept kms PetrolPrice law
0.3348 0.6672 3.3539 0.0082 -1.2224 -0.0963
s.e. 0.0775 0.0612 0.0441 0.0902 0.3839 0.0166
sigma^2 estimated as 0.001476: log likelihood = 349.73, aic = -685.46
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("UKLungDeaths")
> ### * UKLungDeaths
>
> flush(stderr()); flush(stdout())
>
> ### Name: UKLungDeaths
> ### Title: Monthly Deaths from Lung Diseases in the UK
> ### Aliases: UKLungDeaths ldeaths fdeaths mdeaths
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics) # for time
> plot(ldeaths)
> plot(mdeaths, fdeaths)
> ## Better labels:
> yr <- floor(tt <- time(mdeaths))
> plot(mdeaths, fdeaths,
+ xy.labels = paste(month.abb[12*(tt - yr)], yr-1900, sep="'"))
>
>
>
> cleanEx()
> nameEx("UKgas")
> ### * UKgas
>
> flush(stderr()); flush(stdout())
>
> ### Name: UKgas
> ### Title: UK Quarterly Gas Consumption
> ### Aliases: UKgas
> ### Keywords: datasets
>
> ### ** Examples
>
> ## maybe str(UKgas) ; plot(UKgas) ...
>
>
>
> cleanEx()
> nameEx("USArrests")
> ### * USArrests
>
> flush(stderr()); flush(stdout())
>
> ### Name: USArrests
> ### Title: Violent Crime Rates by US State
> ### Aliases: USArrests
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> pairs(USArrests, panel = panel.smooth, main = "USArrests data")
>
>
>
> cleanEx()
> nameEx("USJudgeRatings")
> ### * USJudgeRatings
>
> flush(stderr()); flush(stdout())
>
> ### Name: USJudgeRatings
> ### Title: Lawyers' Ratings of State Judges in the US Superior Court
> ### Aliases: USJudgeRatings
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> pairs(USJudgeRatings, main = "USJudgeRatings data")
>
>
>
> cleanEx()
> nameEx("USPersonalExpenditure")
> ### * USPersonalExpenditure
>
> flush(stderr()); flush(stdout())
>
> ### Name: USPersonalExpenditure
> ### Title: Personal Expenditure Data
> ### Aliases: USPersonalExpenditure
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats) # for medpolish
> USPersonalExpenditure
1940 1945 1950 1955 1960
Food and Tobacco 22.200 44.500 59.60 73.2 86.80
Household Operation 10.500 15.500 29.00 36.5 46.20
Medical and Health 3.530 5.760 9.71 14.0 21.10
Personal Care 1.040 1.980 2.45 3.4 5.40
Private Education 0.341 0.974 1.80 2.6 3.64
> medpolish(log10(USPersonalExpenditure))
1 : 1.126317
2 : 1.032421
Final: 1.032421
Median Polish Results (Dataset: "log10(USPersonalExpenditure)")
Overall: 0.9872192
Row Effects:
Food and Tobacco Household Operation Medical and Health Personal Care
0.7880270 0.4327608 0.0000000 -0.5606543
Private Education
-0.7319467
Column Effects:
1940 1945 1950 1955 1960
-0.4288933 -0.2267967 0.0000000 0.1423128 0.3058289
Residuals:
1940 1945 1950 1955 1960
Food and Tobacco 0.000000 0.0999105 0.000000 -0.053048 -0.142555
Household Operation 0.030103 -0.0028516 0.042418 0.000000 -0.061167
Medical and Health -0.010551 0.0000000 0.000000 0.016596 0.031234
Personal Care 0.019362 0.0968971 -0.037399 -0.037399 0.000000
Private Education -0.293625 -0.0399168 0.000000 0.017388 0.000000
>
>
>
> cleanEx()
> nameEx("VADeaths")
> ### * VADeaths
>
> flush(stderr()); flush(stdout())
>
> ### Name: VADeaths
> ### Title: Death Rates in Virginia (1940)
> ### Aliases: VADeaths
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> n <- length(dr <- c(VADeaths))
> nam <- names(VADeaths)
> d.VAD <- data.frame(
+ Drate = dr,
+ age = rep(ordered(rownames(VADeaths)),length.out=n),
+ gender= gl(2,5,n, labels= c("M", "F")),
+ site = gl(2,10, labels= c("rural", "urban")))
> coplot(Drate ~ as.numeric(age) | gender * site, data = d.VAD,
+ panel = panel.smooth, xlab = "VADeaths data - Given: gender")
> summary(aov.VAD <- aov(Drate ~ .^2, data = d.VAD))
Df Sum Sq Mean Sq F value Pr(>F)
age 4 6288 1572.1 590.858 8.55e-06 ***
gender 1 648 647.5 243.361 9.86e-05 ***
site 1 77 76.8 28.876 0.00579 **
age:gender 4 86 21.6 8.100 0.03358 *
age:site 4 43 10.6 3.996 0.10414
gender:site 1 73 73.0 27.422 0.00636 **
Residuals 4 11 2.7
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(aov.VAD)
> par(opar)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("WWWusage")
> ### * WWWusage
>
> flush(stderr()); flush(stdout())
>
> ### Name: WWWusage
> ### Title: Internet Usage per Minute
> ### Aliases: WWWusage
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> work <- diff(WWWusage)
> par(mfrow = c(2,1)); plot(WWWusage); plot(work)
> ## Not run:
> ##D require(stats)
> ##D aics <- matrix(, 6, 6, dimnames=list(p=0:5, q=0:5))
> ##D for(q in 1:5) aics[1, 1+q] <- arima(WWWusage, c(0,1,q),
> ##D optim.control = list(maxit = 500))$aic
> ##D for(p in 1:5)
> ##D for(q in 0:5) aics[1+p, 1+q] <- arima(WWWusage, c(p,1,q),
> ##D optim.control = list(maxit = 500))$aic
> ##D round(aics - min(aics, na.rm=TRUE), 2)
> ## End(Not run)
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("WorldPhones")
> ### * WorldPhones
>
> flush(stderr()); flush(stdout())
>
> ### Name: WorldPhones
> ### Title: The World's Telephones
> ### Aliases: WorldPhones
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> matplot(rownames(WorldPhones), WorldPhones, type = "b", log = "y",
+ xlab = "Year", ylab = "Number of telephones (1000's)")
> legend(1951.5, 80000, colnames(WorldPhones), col = 1:6, lty = 1:5,
+ pch = rep(21, 7))
> title(main = "World phones data: log scale for response")
>
>
>
> cleanEx()
> nameEx("ability.cov")
> ### * ability.cov
>
> flush(stderr()); flush(stdout())
>
> ### Name: ability.cov
> ### Title: Ability and Intelligence Tests
> ### Aliases: ability.cov
> ### Keywords: datasets
>
> ### ** Examples
>
>
> cleanEx()
> nameEx("airmiles")
> ### * airmiles
>
> flush(stderr()); flush(stdout())
>
> ### Name: airmiles
> ### Title: Passenger Miles on Commercial US Airlines, 1937-1960
> ### Aliases: airmiles
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> plot(airmiles, main = "airmiles data",
+ xlab = "Passenger-miles flown by U.S. commercial airlines", col = 4)
>
>
>
> cleanEx()
> nameEx("airquality")
> ### * airquality
>
> flush(stderr()); flush(stdout())
>
> ### Name: airquality
> ### Title: New York Air Quality Measurements
> ### Aliases: airquality
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> pairs(airquality, panel = panel.smooth, main = "airquality data")
>
>
>
> cleanEx()
> nameEx("anscombe")
> ### * anscombe
>
> flush(stderr()); flush(stdout())
>
> ### Name: anscombe
> ### Title: Anscombe's Quartet of 'Identical' Simple Linear Regressions
> ### Aliases: anscombe
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> summary(anscombe)
x1 x2 x3 x4 y1
Min. : 4.0 Min. : 4.0 Min. : 4.0 Min. : 8 Min. : 4.260
1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 8 1st Qu.: 6.315
Median : 9.0 Median : 9.0 Median : 9.0 Median : 8 Median : 7.580
Mean : 9.0 Mean : 9.0 Mean : 9.0 Mean : 9 Mean : 7.501
3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.: 8 3rd Qu.: 8.570
Max. :14.0 Max. :14.0 Max. :14.0 Max. :19 Max. :10.840
y2 y3 y4
Min. :3.100 Min. : 5.39 Min. : 5.250
1st Qu.:6.695 1st Qu.: 6.25 1st Qu.: 6.170
Median :8.140 Median : 7.11 Median : 7.040
Mean :7.501 Mean : 7.50 Mean : 7.501
3rd Qu.:8.950 3rd Qu.: 7.98 3rd Qu.: 8.190
Max. :9.260 Max. :12.74 Max. :12.500
>
> ##-- now some "magic" to do the 4 regressions in a loop:
> ff <- y ~ x
> for(i in 1:4) {
+ ff[2:3] <- lapply(paste(c("y","x"), i, sep=""), as.name)
+ ## or ff[[2]] <- as.name(paste("y", i, sep=""))
+ ## ff[[3]] <- as.name(paste("x", i, sep=""))
+ assign(paste("lm.",i,sep=""), 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(objects(pattern="lm\\.[1-4]$"), function(n) coef(get(n)))
lm.1 lm.2 lm.3 lm.4
(Intercept) 3.0000909 3.000909 3.0024545 3.0017273
x1 0.5000909 0.500000 0.4997273 0.4999091
> lapply(objects(pattern="lm\\.[1-4]$"),
+ function(n) coef(summary(get(n))))
[[1]]
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
[[2]]
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
[[3]]
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
[[4]]
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=.1+c(4,4,1,1), oma= c(0,0,2,0))
> for(i in 1:4) {
+ ff[2:3] <- lapply(paste(c("y","x"), i, sep=""), as.name)
+ plot(ff, data =anscombe, col="red", pch=21, bg = "orange", cex = 1.2,
+ xlim=c(3,19), ylim=c(3,13))
+ abline(get(paste("lm.",i,sep="")), 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 = .9* 2:0)
> beaver.plot(beaver1, ylim = yl)
> beaver.plot(beaver2, ylim = yl)
> par(op)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("cars")
> ### * cars
>
> flush(stderr()); flush(stdout())
>
> ### Name: cars
> ### Title: Speed and Stopping Distances of Cars
> ### Aliases: cars
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
+ las = 1)
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
> title(main = "cars data")
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
+ las = 1, log = "xy")
> title(main = "cars data (logarithmic scales)")
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
> summary(fm1 <- lm(log(dist) ~ log(speed), data = cars))
Call:
lm(formula = log(dist) ~ log(speed), data = cars)
Residuals:
Min 1Q Median 3Q Max
-1.00215 -0.24578 -0.02898 0.20717 0.88289
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.7297 0.3758 -1.941 0.0581 .
log(speed) 1.6024 0.1395 11.484 2.26e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4053 on 48 degrees of freedom
Multiple R-squared: 0.7331, Adjusted R-squared: 0.7276
F-statistic: 131.9 on 1 and 48 DF, p-value: 2.259e-15
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+ mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
>
> ## An example of polynomial regression
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
+ las = 1, xlim = c(0, 25))
> d <- seq(0, 25, length.out = 200)
> for(degree in 1:4) {
+ fm <- lm(dist ~ poly(speed, degree), data = cars)
+ assign(paste("cars", degree, sep="."), fm)
+ lines(d, predict(fm, data.frame(speed=d)), col = degree)
+ }
> anova(cars.1, cars.2, cars.3, cars.4)
Analysis of Variance Table
Model 1: dist ~ poly(speed, degree)
Model 2: dist ~ poly(speed, degree)
Model 3: dist ~ poly(speed, degree)
Model 4: dist ~ poly(speed, degree)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 48 11354
2 47 10825 1 528.81 2.3108 0.1355
3 46 10634 1 190.35 0.8318 0.3666
4 45 10298 1 336.55 1.4707 0.2316
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("chickwts")
> ### * chickwts
>
> flush(stderr()); flush(stdout())
>
> ### Name: chickwts
> ### Title: Chicken Weights by Feed Type
> ### Aliases: chickwts
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> boxplot(weight ~ feed, data = chickwts, col = "lightgray",
+ varwidth = TRUE, notch = TRUE, main = "chickwt data",
+ ylab = "Weight at six weeks (gm)")
Warning in bxp(list(stats = c(216, 271.5, 342, 373.5, 404, 108, 136, 151.5, :
some notches went outside hinges ('box'): maybe set notch=FALSE
> anova(fm1 <- lm(weight ~ feed, data = chickwts))
Analysis of Variance Table
Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
feed 5 231129 46226 15.365 5.936e-10 ***
Residuals 65 195556 3009
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+ mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("co2")
> ### * co2
>
> flush(stderr()); flush(stdout())
>
> ### Name: co2
> ### Title: Mauna Loa Atmospheric CO2 Concentration
> ### Aliases: co2
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> plot(co2, ylab = expression("Atmospheric concentration of CO"[2]),
+ las = 1)
> title(main = "co2 data set")
>
>
>
> cleanEx()
> nameEx("crimtab")
> ### * crimtab
>
> flush(stderr()); flush(stdout())
>
> ### Name: crimtab
> ### Title: Student's 3000 Criminals Data
> ### Aliases: crimtab
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats)
> dim(crimtab)
[1] 42 22
> utils::str(crimtab)
'table' int [1:42, 1:22] 0 0 0 0 0 0 1 0 0 0 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:42] "9.4" "9.5" "9.6" "9.7" ...
..$ : chr [1:22] "142.24" "144.78" "147.32" "149.86" ...
> ## for nicer printing:
> local({cT <- crimtab
+ colnames(cT) <- substring(colnames(cT), 2,3)
+ print(cT, zero.print = " ")
+ })
42 44 47 49 52 54 57 60 62 65 67 70 72 75 77 80 82 85 87 90 93 95
9.4
9.5 1
9.6
9.7
9.8 1
9.9 1 1 1
10 1 1 2 2 1
10.1 1 3 1 1 1
10.2 2 2 2 1 2 1
10.3 1 1 3 2 2 3 5
10.4 1 1 2 3 3 4 3 3
10.5 1 3 7 6 4 3 1 3 1 1
10.6 1 4 5 9 14 6 3 1 1
10.7 1 2 4 9 14 16 15 7 3 1 2
10.8 2 5 6 14 27 10 7 1 2 1
10.9 2 6 14 24 27 14 10 4 1
11 2 6 12 15 31 37 27 17 10 6
11.1 3 3 12 22 26 24 26 24 7 4 1
11.2 3 2 7 21 30 38 29 27 20 4 1 1
11.3 1 5 10 24 26 39 26 24 7 2
11.4 3 4 9 29 56 58 26 22 10 11
11.5 5 11 17 33 57 38 34 25 11 2
11.6 2 1 4 13 37 39 48 38 27 12 2 2 1
11.7 2 9 17 30 37 48 45 24 9 9 2
11.8 1 2 11 15 35 41 34 29 10 5 1
11.9 1 1 2 12 10 27 32 35 19 10 9 3 1
12 1 4 8 19 42 39 22 16 8 2 2
12.1 2 4 13 22 28 15 27 10 4 1
12.2 1 2 5 6 23 17 16 11 8 1 1
12.3 4 8 10 13 20 23 6 5
12.4 1 1 1 2 7 12 4 7 7 1 1
12.5 1 1 3 12 11 8 6 8 2
12.6 1 3 5 7 8 6 3 1 1
12.7 1 1 7 5 5 8 2 2
12.8 1 2 3 1 8 5 3 1 1
12.9 1 2 2 1 1
13 3 1 1 2 1
13.1 1 1
13.2 1 1 1 3
13.3 1 1
13.4
13.5 1
>
> ## Repeat Student's experiment:
>
> # 1) Reconstitute 3000 raw data for heights in inches and rounded to
> # nearest integer as in Student's paper:
>
> (heIn <- round(as.numeric(colnames(crimtab)) / 2.54))
[1] 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
> d.hei <- data.frame(height = rep(heIn, colSums(crimtab)))
>
> # 2) shuffle the data:
>
> set.seed(1)
> d.hei <- d.hei[sample(1:3000), , drop = FALSE]
>
> # 3) Make 750 samples each of size 4:
>
> d.hei$sample <- as.factor(rep(1:750, each = 4))
>
> # 4) Compute the means and standard deviations (n) for the 750 samples:
>
> h.mean <- with(d.hei, tapply(height, sample, FUN = mean))
> h.sd <- with(d.hei, tapply(height, sample, FUN = sd)) * sqrt(3/4)
>
> # 5) Compute the difference between the mean of each sample and
> # the mean of the population and then divide by the
> # standard deviation of the sample:
>
> zobs <- (h.mean - mean(d.hei[,"height"]))/h.sd
>
> # 6) Replace infinite values by +/- 6 as in Student's paper:
>
> zobs[infZ <- is.infinite(zobs)] # 3 of them
73 312 674
-Inf Inf -Inf
> zobs[infZ] <- 6 * sign(zobs[infZ])
>
> # 7) Plot the distribution:
>
> require(grDevices); require(graphics)
> hist(x = zobs, probability = TRUE, xlab = "Student's z",
+ col = grey(0.8), border = grey(0.5),
+ main = "Distribution of Student's z score for 'crimtab' data")
>
>
>
> cleanEx()
> nameEx("discoveries")
> ### * discoveries
>
> flush(stderr()); flush(stdout())
>
> ### Name: discoveries
> ### Title: Yearly Numbers of Important Discoveries
> ### Aliases: discoveries
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> plot(discoveries, ylab = "Number of important discoveries",
+ las = 1)
> title(main = "discoveries data set")
>
>
>
> cleanEx()
> nameEx("esoph")
> ### * esoph
>
> flush(stderr()); flush(stdout())
>
> ### Name: esoph
> ### Title: Smoking, Alcohol and (O)esophageal Cancer
> ### Aliases: esoph
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats)
> require(graphics) # for mosaicplot
> summary(esoph)
agegp alcgp tobgp ncases ncontrols
25-34:15 0-39g/day:23 0-9g/day:24 Min. : 0.000 Min. : 1.00
35-44:15 40-79 :23 10-19 :24 1st Qu.: 0.000 1st Qu.: 3.00
45-54:16 80-119 :21 20-29 :20 Median : 1.000 Median : 6.00
55-64:16 120+ :21 30+ :20 Mean : 2.273 Mean :11.08
65-74:15 3rd Qu.: 4.000 3rd Qu.:14.00
75+ :11 Max. :17.000 Max. :60.00
> ## effects of alcohol, tobacco and interaction, age-adjusted
> model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp,
+ data = esoph, family = binomial())
> anova(model1)
Analysis of Deviance Table
Model: binomial, link: logit
Response: cbind(ncases, ncontrols)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev
NULL 87 227.241
agegp 5 88.128 82 139.112
tobgp 3 19.085 79 120.028
alcgp 3 66.054 76 53.973
tobgp:alcgp 9 6.489 67 47.484
> ## Try a linear effect of alcohol and tobacco
> model2 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp)
+ + unclass(alcgp),
+ data = esoph, family = binomial())
> summary(model2)
Call:
glm(formula = cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) +
unclass(alcgp), family = binomial(), data = esoph)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7628 -0.6426 -0.2709 0.3043 2.0421
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.01097 0.31224 -12.846 < 2e-16 ***
agegp.L 2.96113 0.65092 4.549 5.39e-06 ***
agegp.Q -1.33735 0.58918 -2.270 0.02322 *
agegp.C 0.15292 0.44792 0.341 0.73281
agegp^4 0.06668 0.30776 0.217 0.82848
agegp^5 -0.20288 0.19523 -1.039 0.29872
unclass(tobgp) 0.26162 0.08198 3.191 0.00142 **
unclass(alcgp) 0.65308 0.08452 7.727 1.10e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 227.241 on 87 degrees of freedom
Residual deviance: 59.277 on 80 degrees of freedom
AIC: 222.76
Number of Fisher Scoring iterations: 6
> ## Re-arrange data for a mosaic plot
> ttt <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
> o <- with(esoph, order(tobgp, alcgp, agegp))
> ttt[ttt == 1] <- esoph$ncases[o]
> tt1 <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
> tt1[tt1 == 1] <- esoph$ncontrols[o]
> tt <- array(c(ttt, tt1), c(dim(ttt),2),
+ c(dimnames(ttt), list(c("Cancer", "control"))))
> mosaicplot(tt, main = "esoph data set", color = TRUE)
>
>
>
> cleanEx()
> nameEx("euro")
> ### * euro
>
> flush(stderr()); flush(stdout())
>
> ### Name: euro
> ### Title: Conversion Rates of Euro Currencies
> ### Aliases: euro euro.cross
> ### Keywords: datasets
>
> ### ** Examples
>
> cbind(euro)
euro
ATS 13.760300
BEF 40.339900
DEM 1.955830
ESP 166.386000
FIM 5.945730
FRF 6.559570
IEP 0.787564
ITL 1936.270000
LUF 40.339900
NLG 2.203710
PTE 200.482000
>
> ## These relations hold:
> euro == signif(euro,6) # [6 digit precision in Euro's definition]
ATS BEF DEM ESP FIM FRF IEP ITL LUF NLG PTE
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> all(euro.cross == outer(1/euro, euro))
[1] TRUE
>
> ## Convert 20 Euro to Belgian Franc
> 20 * euro["BEF"]
BEF
806.798
> ## Convert 20 Austrian Schilling to Euro
> 20 / euro["ATS"]
ATS
1.453457
> ## Convert 20 Spanish Pesetas to Italian Lira
> 20 * euro.cross["ESP", "ITL"]
[1] 232.7443
>
> require(graphics)
> dotchart(euro,
+ main = "euro data: 1 Euro in currency unit")
> dotchart(1/euro,
+ main = "euro data: 1 currency unit in Euros")
> dotchart(log(euro, 10),
+ main = "euro data: log10(1 Euro in currency unit)")
>
>
>
> cleanEx()
> nameEx("faithful")
> ### * faithful
>
> flush(stderr()); flush(stdout())
>
> ### Name: faithful
> ### Title: Old Faithful Geyser Data
> ### Aliases: faithful
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> f.tit <- "faithful data: Eruptions of Old Faithful"
>
> ne60 <- round(e60 <- 60 * faithful$eruptions)
> all.equal(e60, ne60) # relative diff. ~ 1/10000
[1] "Mean relative difference: 9.515332e-05"
> table(zapsmall(abs(e60 - ne60))) # 0, 0.02 or 0.04
0 0.02 0.04
106 163 3
> faithful$better.eruptions <- ne60 / 60
> te <- table(ne60)
> te[te >= 4] # (too) many multiples of 5 !
ne60
105 108 110 112 113 120 216 230 240 245 249 250 255 260 261 262 265 270 272 275
6 4 7 8 4 4 4 5 6 5 4 4 4 5 4 4 4 8 5 4
276 282 288
4 6 6
> plot(names(te), te, type="h", main = f.tit, xlab = "Eruption time (sec)")
>
> plot(faithful[, -3], main = f.tit,
+ xlab = "Eruption time (min)",
+ ylab = "Waiting time to next eruption (min)")
> lines(lowess(faithful$eruptions, faithful$waiting, f = 2/3, iter = 3),
+ col = "red")
>
>
>
> cleanEx()
> nameEx("freeny")
> ### * freeny
>
> flush(stderr()); flush(stdout())
>
> ### Name: freeny
> ### Title: Freeny's Revenue Data
> ### Aliases: freeny freeny.x freeny.y
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> summary(freeny)
y lag.quarterly.revenue price.index income.level
Min. :8.791 Min. :8.791 Min. :4.278 Min. :5.821
1st Qu.:9.045 1st Qu.:9.020 1st Qu.:4.392 1st Qu.:5.948
Median :9.314 Median :9.284 Median :4.510 Median :6.061
Mean :9.306 Mean :9.281 Mean :4.496 Mean :6.039
3rd Qu.:9.591 3rd Qu.:9.561 3rd Qu.:4.605 3rd Qu.:6.139
Max. :9.794 Max. :9.775 Max. :4.710 Max. :6.200
market.potential
Min. :12.97
1st Qu.:13.01
Median :13.07
Mean :13.07
3rd Qu.:13.12
Max. :13.17
> pairs(freeny, main = "freeny data")
> # gives warning: freeny$y has class "ts"
>
> summary(fm1 <- lm(y ~ ., data = freeny))
Call:
lm(formula = y ~ ., data = freeny)
Residuals:
Min 1Q Median 3Q Max
-0.0259426 -0.0101033 0.0003824 0.0103236 0.0267124
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10.4726 6.0217 -1.739 0.0911 .
lag.quarterly.revenue 0.1239 0.1424 0.870 0.3904
price.index -0.7542 0.1607 -4.693 4.28e-05 ***
income.level 0.7675 0.1339 5.730 1.93e-06 ***
market.potential 1.3306 0.5093 2.613 0.0133 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01473 on 34 degrees of freedom
Multiple R-squared: 0.9981, Adjusted R-squared: 0.9978
F-statistic: 4354 on 4 and 34 DF, p-value: < 2.2e-16
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+ mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("infert")
> ### * infert
>
> flush(stderr()); flush(stdout())
>
> ### Name: infert
> ### Title: Infertility after Spontaneous and Induced Abortion
> ### Aliases: infert
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats)
> model1 <- glm(case ~ spontaneous+induced, data=infert,family=binomial())
> summary(model1)
Call:
glm(formula = case ~ spontaneous + induced, family = binomial(),
data = infert)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6678 -0.8360 -0.5772 0.9030 1.9362
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.7079 0.2677 -6.380 1.78e-10 ***
spontaneous 1.1972 0.2116 5.657 1.54e-08 ***
induced 0.4181 0.2056 2.033 0.042 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 316.17 on 247 degrees of freedom
Residual deviance: 279.61 on 245 degrees of freedom
AIC: 285.61
Number of Fisher Scoring iterations: 4
> ## adjusted for other potential confounders:
> summary(model2 <- glm(case ~ age+parity+education+spontaneous+induced,
+ data=infert,family=binomial()))
Call:
glm(formula = case ~ age + parity + education + spontaneous +
induced, family = binomial(), data = infert)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7603 -0.8162 -0.4956 0.8349 2.6536
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.14924 1.41220 -0.814 0.4158
age 0.03958 0.03120 1.269 0.2046
parity -0.82828 0.19649 -4.215 2.49e-05 ***
education6-11yrs -1.04424 0.79255 -1.318 0.1876
education12+ yrs -1.40321 0.83416 -1.682 0.0925 .
spontaneous 2.04591 0.31016 6.596 4.21e-11 ***
induced 1.28876 0.30146 4.275 1.91e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 316.17 on 247 degrees of freedom
Residual deviance: 257.80 on 241 degrees of freedom
AIC: 271.8
Number of Fisher Scoring iterations: 4
> ## Really should be analysed by conditional logistic regression
> ## which is in the survival package
> if(require(survival)){
+ model3 <- clogit(case~spontaneous+induced+strata(stratum),data=infert)
+ print(summary(model3))
+ detach()# survival (conflicts)
+ }
Loading required package: survival
Loading required package: splines
Call:
coxph(formula = Surv(rep(1, 248L), case) ~ spontaneous + induced +
strata(stratum), data = infert, method = "exact")
n= 248, number of events= 83
coef exp(coef) se(coef) z Pr(>|z|)
spontaneous 1.9859 7.2854 0.3524 5.635 1.75e-08 ***
induced 1.4090 4.0919 0.3607 3.906 9.38e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
spontaneous 7.285 0.1373 3.651 14.536
induced 4.092 0.2444 2.018 8.298
Rsquare= 0.193 (max possible= 0.519 )
Likelihood ratio test= 53.15 on 2 df, p=2.869e-12
Wald test = 31.84 on 2 df, p=1.221e-07
Score (logrank) test = 48.44 on 2 df, p=3.032e-11
>
>
>
> cleanEx()
detaching ‘package:splines’
> nameEx("iris")
> ### * iris
>
> flush(stderr()); flush(stdout())
>
> ### Name: iris
> ### Title: Edgar Anderson's Iris Data
> ### Aliases: iris iris3
> ### Keywords: datasets
>
> ### ** Examples
>
> dni3 <- dimnames(iris3)
> ii <- data.frame(matrix(aperm(iris3, c(1,3,2)), ncol=4,
+ dimnames = list(NULL, sub(" L.",".Length",
+ sub(" W.",".Width", dni3[[2]])))),
+ Species = gl(3, 50, labels=sub("S", "s", sub("V", "v", dni3[[3]]))))
> all.equal(ii, iris) # TRUE
[1] TRUE
>
>
>
> cleanEx()
> nameEx("islands")
> ### * islands
>
> flush(stderr()); flush(stdout())
>
> ### Name: islands
> ### Title: Areas of the World's Major Landmasses
> ### Aliases: islands
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> dotchart(log(islands, 10),
+ main = "islands data: log10(area) (log10(sq. miles))")
> dotchart(log(islands[order(islands)], 10),
+ main = "islands data: log10(area) (log10(sq. miles))")
>
>
>
> cleanEx()
> nameEx("longley")
> ### * longley
>
> flush(stderr()); flush(stdout())
>
> ### Name: longley
> ### Title: Longley's Economic Regression Data
> ### Aliases: longley
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> ## give the data set in the form it is used in S-PLUS:
> longley.x <- data.matrix(longley[, 1:6])
> longley.y <- longley[, "Employed"]
> pairs(longley, main = "longley data")
> summary(fm1 <- lm(Employed ~ ., data = longley))
Call:
lm(formula = Employed ~ ., data = longley)
Residuals:
Min 1Q Median 3Q Max
-0.41011 -0.15767 -0.02816 0.10155 0.45539
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.482e+03 8.904e+02 -3.911 0.003560 **
GNP.deflator 1.506e-02 8.492e-02 0.177 0.863141
GNP -3.582e-02 3.349e-02 -1.070 0.312681
Unemployed -2.020e-02 4.884e-03 -4.136 0.002535 **
Armed.Forces -1.033e-02 2.143e-03 -4.822 0.000944 ***
Population -5.110e-02 2.261e-01 -0.226 0.826212
Year 1.829e+00 4.555e-01 4.016 0.003037 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3049 on 9 degrees of freedom
Multiple R-squared: 0.9955, Adjusted R-squared: 0.9925
F-statistic: 330.3 on 6 and 9 DF, p-value: 4.984e-10
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+ mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("morley")
> ### * morley
>
> flush(stderr()); flush(stdout())
>
> ### Name: morley
> ### Title: Michaelson-Morley Speed of Light Data
> ### Aliases: morley
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> morley$Expt <- factor(morley$Expt)
> morley$Run <- factor(morley$Run)
>
> xtabs(~ Expt + Run, data = morley)# 5 x 20 balanced (two-way)
Run
Expt 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> plot(Speed ~ Expt, data = morley,
+ main = "Speed of Light Data", xlab = "Experiment No.")
> fm <- aov(Speed ~ Run + Expt, data = morley)
> summary(fm)
Df Sum Sq Mean Sq F value Pr(>F)
Run 19 113344 5965 1.105 0.36321
Expt 4 94514 23629 4.378 0.00307 **
Residuals 76 410166 5397
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> fm0 <- update(fm, . ~ . - Run)
> anova(fm0, fm)
Analysis of Variance Table
Model 1: Speed ~ Expt
Model 2: Speed ~ Run + Expt
Res.Df RSS Df Sum of Sq F Pr(>F)
1 95 523510
2 76 410166 19 113344 1.1053 0.3632
>
>
>
> cleanEx()
> nameEx("mtcars")
> ### * mtcars
>
> flush(stderr()); flush(stdout())
>
> ### Name: mtcars
> ### Title: Motor Trend Car Road Tests
> ### Aliases: mtcars
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> pairs(mtcars, main = "mtcars data")
> coplot(mpg ~ disp | as.factor(cyl), data = mtcars,
+ panel = panel.smooth, rows = 1)
>
>
>
> cleanEx()
> nameEx("nhtemp")
> ### * nhtemp
>
> flush(stderr()); flush(stdout())
>
> ### Name: nhtemp
> ### Title: Average Yearly Temperatures in New Haven
> ### Aliases: nhtemp
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> plot(nhtemp, main = "nhtemp data",
+ ylab = "Mean annual temperature in New Haven, CT (deg. F)")
>
>
>
> cleanEx()
> nameEx("nottem")
> ### * nottem
>
> flush(stderr()); flush(stdout())
>
> ### Name: nottem
> ### Title: Average Monthly Temperatures at Nottingham, 1920-1939
> ### Aliases: nottem
> ### Keywords: datasets
>
> ### ** Examples
>
> ## Not run:
> ##D require(stats); require(graphics)
> ##D nott <- window(nottem, end=c(1936,12))
> ##D fit <- arima(nott,order=c(1,0,0), list(order=c(2,1,0), period=12))
> ##D nott.fore <- predict(fit, n.ahead=36)
> ##D ts.plot(nott, nott.fore$pred, nott.fore$pred+2*nott.fore$se,
> ##D nott.fore$pred-2*nott.fore$se, gpars=list(col=c(1,1,4,4)))
> ## End(Not run)
>
>
> cleanEx()
> nameEx("occupationalStatus")
> ### * occupationalStatus
>
> flush(stderr()); flush(stdout())
>
> ### Name: occupationalStatus
> ### Title: Occupational Status of Fathers and their Sons
> ### Aliases: occupationalStatus
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
>
> plot(occupationalStatus)
>
> ## Fit a uniform association model separating diagonal effects
> Diag <- as.factor(diag(1:8))
> Rscore <- scale(as.numeric(row(occupationalStatus)), scale = FALSE)
> Cscore <- scale(as.numeric(col(occupationalStatus)), scale = FALSE)
> modUnif <- glm(Freq ~ origin + destination + Diag + Rscore:Cscore,
+ family = poisson, data = occupationalStatus)
>
> summary(modUnif)
Call:
glm(formula = Freq ~ origin + destination + Diag + Rscore:Cscore,
family = poisson, data = occupationalStatus)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6521 -0.6267 0.0000 0.5913 2.0964
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.568592 0.183358 3.101 0.001929 **
origin2 0.431314 0.149415 2.887 0.003893 **
origin3 1.461862 0.131141 11.147 < 2e-16 ***
origin4 1.788731 0.126588 14.130 < 2e-16 ***
origin5 0.441011 0.144844 3.045 0.002329 **
origin6 2.491735 0.121219 20.556 < 2e-16 ***
origin7 1.127536 0.129032 8.738 < 2e-16 ***
origin8 0.796445 0.131863 6.040 1.54e-09 ***
destination2 0.873202 0.166902 5.232 1.68e-07 ***
destination3 1.813992 0.153861 11.790 < 2e-16 ***
destination4 2.082515 0.150887 13.802 < 2e-16 ***
destination5 1.366383 0.155590 8.782 < 2e-16 ***
destination6 2.816369 0.146054 19.283 < 2e-16 ***
destination7 1.903918 0.147810 12.881 < 2e-16 ***
destination8 1.398585 0.151658 9.222 < 2e-16 ***
Diag1 1.665495 0.237383 7.016 2.28e-12 ***
Diag2 0.959681 0.212122 4.524 6.06e-06 ***
Diag3 0.021750 0.156530 0.139 0.889490
Diag4 0.226399 0.124364 1.820 0.068689 .
Diag5 0.808646 0.229754 3.520 0.000432 ***
Diag6 0.132277 0.077191 1.714 0.086597 .
Diag7 0.506709 0.115936 4.371 1.24e-05 ***
Diag8 0.221880 0.134803 1.646 0.099771 .
Rscore:Cscore 0.136974 0.007489 18.289 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 4679.004 on 63 degrees of freedom
Residual deviance: 58.436 on 40 degrees of freedom
AIC: 428.78
Number of Fisher Scoring iterations: 4
> plot(modUnif) # 4 plots, with warning about h_ii ~= 1
Warning: Not plotting observations with leverage one:
1, 10, 19, 28, 37, 46, 55, 64
Warning: Not plotting observations with leverage one:
1, 10, 19, 28, 37, 46, 55, 64
>
>
>
> cleanEx()
> nameEx("precip")
> ### * precip
>
> flush(stderr()); flush(stdout())
>
> ### Name: precip
> ### Title: Annual Precipitation in US Cities
> ### Aliases: precip
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> dotchart(precip[order(precip)], main = "precip data")
> title(sub = "Average annual precipitation (in.)")
>
>
>
> cleanEx()
> nameEx("presidents")
> ### * presidents
>
> flush(stderr()); flush(stdout())
>
> ### Name: presidents
> ### Title: Quarterly Approval Ratings of US Presidents
> ### Aliases: presidents
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> plot(presidents, las = 1, ylab = "Approval rating (%)",
+ main = "presidents data")
>
>
>
> cleanEx()
> nameEx("pressure")
> ### * pressure
>
> flush(stderr()); flush(stdout())
>
> ### Name: pressure
> ### Title: Vapor Pressure of Mercury as a Function of Temperature
> ### Aliases: pressure
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> plot(pressure, xlab = "Temperature (deg C)",
+ ylab = "Pressure (mm of Hg)",
+ main = "pressure data: Vapor Pressure of Mercury")
> plot(pressure, xlab = "Temperature (deg C)", log = "y",
+ ylab = "Pressure (mm of Hg)",
+ main = "pressure data: Vapor Pressure of Mercury")
>
>
>
> cleanEx()
> nameEx("quakes")
> ### * quakes
>
> flush(stderr()); flush(stdout())
>
> ### Name: quakes
> ### Title: Locations of Earthquakes off Fiji
> ### Aliases: quakes
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> pairs(quakes, main = "Fiji Earthquakes, N = 1000", cex.main=1.2, pch=".")
>
>
>
> cleanEx()
> nameEx("randu")
> ### * randu
>
> flush(stderr()); flush(stdout())
>
> ### Name: randu
> ### Title: Random Numbers from Congruential Generator RANDU
> ### Aliases: randu
> ### Keywords: datasets
>
> ### ** Examples
>
> ## Not run:
> ##D ## We could re-generate the dataset by the following R code
> ##D seed <- as.double(1)
> ##D RANDU <- function() {
> ##D seed <<- ((2^16 + 3) * seed) %% (2^31)
> ##D seed/(2^31)
> ##D }
> ##D for(i in 1:400) {
> ##D U <- c(RANDU(), RANDU(), RANDU(), RANDU(), RANDU())
> ##D print(round(U[1:3], 6))
> ##D }
> ## End(Not run)
>
>
>
> cleanEx()
> nameEx("sleep")
> ### * sleep
>
> flush(stderr()); flush(stdout())
>
> ### Name: sleep
> ### Title: Student's Sleep Data
> ### Aliases: sleep
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats)
> ## Student's paired t-test
> with(sleep,
+ t.test(extra[group == 1],
+ extra[group == 2], paired = TRUE))
Paired t-test
data: extra[group == 1] and extra[group == 2]
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.4598858 -0.7001142
sample estimates:
mean of the differences
-1.58
>
>
>
> cleanEx()
> nameEx("stackloss")
> ### * stackloss
>
> flush(stderr()); flush(stdout())
>
> ### Name: stackloss
> ### Title: Brownlee's Stack Loss Plant Data
> ### Aliases: stackloss stack.loss stack.x
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats)
> summary(lm.stack <- lm(stack.loss ~ stack.x))
Call:
lm(formula = stack.loss ~ stack.x)
Residuals:
Min 1Q Median 3Q Max
-7.2377 -1.7117 -0.4551 2.3614 5.6978
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -39.9197 11.8960 -3.356 0.00375 **
stack.xAir.Flow 0.7156 0.1349 5.307 5.8e-05 ***
stack.xWater.Temp 1.2953 0.3680 3.520 0.00263 **
stack.xAcid.Conc. -0.1521 0.1563 -0.973 0.34405
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.243 on 17 degrees of freedom
Multiple R-squared: 0.9136, Adjusted R-squared: 0.8983
F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09
>
>
>
> cleanEx()
> nameEx("sunspot.month")
> ### * sunspot.month
>
> flush(stderr()); flush(stdout())
>
> ### Name: sunspot.month
> ### Title: Monthly Sunspot Data, 1749-1997
> ### Aliases: sunspot.month
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> ## Compare the monthly series
> plot (sunspot.month, main = "sunspot.month [stats]", col = 2)
> lines(sunspots) # "very barely" see something
>
> ## Now look at the difference :
> all(tsp(sunspots) [c(1,3)] ==
+ tsp(sunspot.month)[c(1,3)]) ## Start & Periodicity are the same
[1] TRUE
> n1 <- length(sunspots)
> table(eq <- sunspots == sunspot.month[1:n1]) #> 132 are different !
FALSE TRUE
132 2688
> i <- which(!eq)
> rug(time(eq)[i])
> s1 <- sunspots[i] ; s2 <- sunspot.month[i]
> cbind(i = i, sunspots = s1, ss.month = s2,
+ perc.diff = round(100*2*abs(s1-s2)/(s1+s2), 1))
i sunspots ss.month perc.diff
[1,] 55 22.2 22.0 0.9
[2,] 838 31.7 31.6 0.3
[3,] 841 32.5 32.8 0.9
[4,] 862 9.0 8.9 1.1
[5,] 864 9.7 9.1 6.4
[6,] 866 4.3 4.2 2.4
[7,] 876 0.0 0.2 200.0
[8,] 901 21.6 21.7 0.5
[9,] 917 15.4 15.5 0.6
[10,] 920 25.4 25.7 1.2
[11,] 943 42.9 42.3 1.4
[12,] 946 57.2 56.1 1.9
[13,] 955 54.3 54.2 0.2
[14,] 960 46.6 46.9 0.6
[15,] 965 67.5 67.4 0.1
[16,] 968 78.3 77.6 0.9
[17,] 976 107.1 106.3 0.7
[18,] 988 54.6 54.5 0.2
[19,] 992 54.9 55.0 0.2
[20,] 994 46.2 46.3 0.2
[21,] 998 55.5 55.6 0.2
[22,] 1003 13.9 14.0 0.7
[23,] 1047 98.1 98.2 0.1
[24,] 1061 111.3 111.7 0.4
[25,] 1081 107.6 105.6 1.9
[26,] 1087 84.7 84.8 0.1
[27,] 1090 90.8 90.9 0.1
[28,] 1092 63.6 63.7 0.2
[29,] 1095 55.5 67.8 20.0
[30,] 1102 49.8 55.0 9.9
[31,] 1105 24.0 24.1 0.4
[32,] 1108 42.6 40.2 5.8
[33,] 1109 67.4 67.5 0.1
[34,] 1113 35.1 36.5 3.9
[35,] 1124 26.5 26.6 0.4
[36,] 1125 18.5 18.4 0.5
[37,] 1132 8.8 9.5 7.7
[38,] 1145 12.0 11.6 3.4
[39,] 1149 6.9 7.0 1.4
[40,] 1156 56.9 57.0 0.2
[41,] 1168 69.2 69.3 0.1
[42,] 1185 161.2 160.9 0.2
[43,] 1191 108.9 108.6 0.3
[44,] 1194 123.8 129.0 4.1
[45,] 1196 132.5 132.6 0.1
[46,] 1200 159.9 159.5 0.3
[47,] 1201 156.7 157.0 0.2
[48,] 1203 96.5 96.2 0.3
[49,] 1206 81.2 81.1 0.1
[50,] 1208 61.3 67.7 9.9
[51,] 1211 99.7 99.0 0.7
[52,] 1224 60.0 61.0 1.7
[53,] 1235 50.9 51.0 0.2
[54,] 1238 67.5 66.4 1.6
[55,] 1243 42.0 42.1 0.2
[56,] 1264 26.4 26.5 0.4
[57,] 1270 12.7 12.6 0.8
[58,] 1272 21.4 21.6 0.9
[59,] 1282 9.7 9.6 1.0
[60,] 1283 4.3 4.2 2.4
[61,] 1290 5.0 5.2 3.9
[62,] 1301 29.2 28.6 2.1
[63,] 1333 81.5 82.4 1.1
[64,] 1334 88.0 88.3 0.3
[65,] 1346 77.8 77.7 0.1
[66,] 1350 87.8 88.1 0.3
[67,] 1366 42.0 41.9 0.2
[68,] 1424 4.9 4.8 2.1
[69,] 1427 9.3 9.6 3.2
[70,] 1430 15.8 15.7 0.6
[71,] 1435 28.6 29.0 1.4
[72,] 1437 43.8 47.2 7.5
[73,] 1438 61.7 61.6 0.2
[74,] 1442 59.3 59.9 1.0
[75,] 1445 104.0 103.9 0.1
[76,] 1450 59.4 59.3 0.2
[77,] 1451 77.4 78.1 0.9
[78,] 1455 159.4 157.6 1.1
[79,] 1472 110.0 110.1 0.1
[80,] 1476 90.3 90.4 0.1
[81,] 1486 103.5 102.6 0.9
[82,] 1497 47.5 47.1 0.8
[83,] 1498 47.4 47.1 0.6
[84,] 1514 22.2 21.5 3.2
[85,] 1527 31.2 30.6 1.9
[86,] 1539 11.7 11.9 1.7
[87,] 1541 21.2 21.6 1.9
[88,] 1542 13.4 14.2 5.8
[89,] 1543 5.9 6.0 1.7
[90,] 1545 16.4 16.9 3.0
[91,] 1547 14.5 14.2 2.1
[92,] 1548 2.3 2.2 4.4
[93,] 1550 6.0 6.6 9.5
[94,] 1553 5.8 5.9 1.7
[95,] 1561 0.8 1.0 22.2
[96,] 1571 12.9 13.1 1.5
[97,] 1572 7.2 7.3 1.4
[98,] 1574 27.5 27.2 1.1
[99,] 1575 19.5 19.3 1.0
[100,] 1576 19.3 19.5 1.0
[101,] 1588 51.7 51.6 0.2
[102,] 1592 58.0 58.4 0.7
[103,] 1594 64.0 64.4 0.6
[104,] 1598 69.3 69.5 0.3
[105,] 1599 67.5 66.8 1.0
[106,] 1613 32.1 31.5 1.9
[107,] 1614 76.5 76.3 0.3
[108,] 1623 86.8 87.5 0.8
[109,] 1643 33.3 30.9 7.5
[110,] 1656 12.4 13.0 4.7
[111,] 1683 7.0 6.7 4.4
[112,] 1687 9.7 9.4 3.1
[113,] 1712 33.2 33.0 0.6
[114,] 1716 32.3 32.5 0.6
[115,] 1723 76.8 76.5 0.4
[116,] 1734 88.2 89.9 1.9
[117,] 1735 88.8 88.6 0.2
[118,] 1738 79.7 80.0 0.4
[119,] 1774 28.4 28.7 1.1
[120,] 1837 5.2 5.5 5.6
[121,] 2126 70.0 69.9 0.1
[122,] 2153 76.9 77.0 0.1
[123,] 2162 64.1 62.8 2.0
[124,] 2174 49.2 49.9 1.4
[125,] 2233 18.9 18.6 1.6
[126,] 2315 38.3 38.4 0.3
[127,] 2378 113.4 133.4 16.2
[128,] 2427 59.9 55.9 6.9
[129,] 2552 55.9 55.8 0.2
[130,] 2556 40.0 39.9 0.3
[131,] 2694 39.5 37.5 5.2
[132,] 2790 90.0 90.9 1.0
>
>
>
> cleanEx()
> nameEx("sunspots")
> ### * sunspots
>
> flush(stderr()); flush(stdout())
>
> ### Name: sunspots
> ### Title: Monthly Sunspot Numbers, 1749-1983
> ### Aliases: sunspots
> ### Keywords: datasets
>
> ### ** Examples
>
> require(graphics)
> plot(sunspots, main = "sunspots data", xlab = "Year",
+ ylab = "Monthly sunspot numbers")
>
>
>
> cleanEx()
> nameEx("swiss")
> ### * swiss
>
> flush(stderr()); flush(stdout())
>
> ### Name: swiss
> ### Title: Swiss Fertility and Socioeconomic Indicators (1888) Data
> ### Aliases: swiss
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> pairs(swiss, panel = panel.smooth, main = "swiss data",
+ col = 3 + (swiss$Catholic > 50))
> summary(lm(Fertility ~ . , data = swiss))
Call:
lm(formula = Fertility ~ ., data = swiss)
Residuals:
Min 1Q Median 3Q Max
-15.2743 -5.2617 0.5032 4.1198 15.3213
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
Agriculture -0.17211 0.07030 -2.448 0.01873 *
Examination -0.25801 0.25388 -1.016 0.31546
Education -0.87094 0.18303 -4.758 2.43e-05 ***
Catholic 0.10412 0.03526 2.953 0.00519 **
Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.165 on 41 degrees of freedom
Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10
>
>
>
> cleanEx()
> nameEx("trees")
> ### * trees
>
> flush(stderr()); flush(stdout())
>
> ### Name: trees
> ### Title: Girth, Height and Volume for Black Cherry Trees
> ### Aliases: trees
> ### Keywords: datasets
>
> ### ** Examples
>
> require(stats); require(graphics)
> pairs(trees, panel = panel.smooth, main = "trees data")
> plot(Volume ~ Girth, data = trees, log = "xy")
> coplot(log(Volume) ~ log(Girth) | Height, data = trees,
+ panel = panel.smooth)
> summary(fm1 <- lm(log(Volume) ~ log(Girth), data = trees))
Call:
lm(formula = log(Volume) ~ log(Girth), data = trees)
Residuals:
Min 1Q Median 3Q Max
-0.205999 -0.068702 0.001011 0.072585 0.247963
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.35332 0.23066 -10.20 4.18e-11 ***
log(Girth) 2.19997 0.08983 24.49 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.115 on 29 degrees of freedom
Multiple R-squared: 0.9539, Adjusted R-squared: 0.9523
F-statistic: 599.7 on 1 and 29 DF, p-value: < 2.2e-16
> summary(fm2 <- update(fm1, ~ . + log(Height), data = trees))
Call:
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
Residuals:
Min 1Q Median 3Q Max
-0.168561 -0.048488 0.002431 0.063637 0.129223
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.63162 0.79979 -8.292 5.06e-09 ***
log(Girth) 1.98265 0.07501 26.432 < 2e-16 ***
log(Height) 1.11712 0.20444 5.464 7.81e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.08139 on 28 degrees of freedom
Multiple R-squared: 0.9777, Adjusted R-squared: 0.9761
F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16
> step(fm2)
Start: AIC=-152.69
log(Volume) ~ log(Girth) + log(Height)
Df Sum of Sq RSS AIC
<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>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 2.645 0.094 2.778 0 0
> grDevices::dev.off()
null device
1
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')