| 56186 |
murdoch |
1 |
|
| 89316 |
smeyer |
2 |
R Under development (unstable) (2026-01-21 r89314) -- "Unsuffered Consequences"
|
| 89306 |
smeyer |
3 |
Copyright (C) 2026 The R Foundation for Statistical Computing
|
| 87723 |
maechler |
4 |
Platform: x86_64-pc-linux-gnu
|
| 56186 |
murdoch |
5 |
|
|
|
6 |
R is free software and comes with ABSOLUTELY NO WARRANTY.
|
|
|
7 |
You are welcome to redistribute it under certain conditions.
|
|
|
8 |
Type 'license()' or 'licence()' for distribution details.
|
|
|
9 |
|
|
|
10 |
R is a collaborative project with many contributors.
|
|
|
11 |
Type 'contributors()' for more information and
|
|
|
12 |
'citation()' on how to cite R or R packages in publications.
|
|
|
13 |
|
|
|
14 |
Type 'demo()' for some demos, 'help()' for on-line help, or
|
|
|
15 |
'help.start()' for an HTML browser interface to help.
|
|
|
16 |
Type 'q()' to quit R.
|
|
|
17 |
|
|
|
18 |
> pkgname <- "datasets"
|
|
|
19 |
> source(file.path(R.home("share"), "R", "examples-header.R"))
|
|
|
20 |
> options(warn = 1)
|
|
|
21 |
> library('datasets')
|
|
|
22 |
>
|
| 61787 |
ripley |
23 |
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
|
| 73819 |
hornik |
24 |
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
|
| 56186 |
murdoch |
25 |
> cleanEx()
|
|
|
26 |
> nameEx("AirPassengers")
|
|
|
27 |
> ### * AirPassengers
|
|
|
28 |
>
|
|
|
29 |
> flush(stderr()); flush(stdout())
|
|
|
30 |
>
|
|
|
31 |
> ### Name: AirPassengers
|
|
|
32 |
> ### Title: Monthly Airline Passenger Numbers 1949-1960
|
|
|
33 |
> ### Aliases: AirPassengers
|
|
|
34 |
> ### Keywords: datasets
|
|
|
35 |
>
|
|
|
36 |
> ### ** Examples
|
|
|
37 |
>
|
| 86049 |
smeyer |
38 |
> ## The classic 'airline model', by full ML
|
|
|
39 |
> (fit <- arima(log10(AirPassengers), c(0, 1, 1),
|
|
|
40 |
+ seasonal = list(order = c(0, 1, 1), period = 12)))
|
|
|
41 |
|
|
|
42 |
Call:
|
|
|
43 |
arima(x = log10(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0,
|
|
|
44 |
1, 1), period = 12))
|
|
|
45 |
|
|
|
46 |
Coefficients:
|
|
|
47 |
ma1 sma1
|
|
|
48 |
-0.4018 -0.5569
|
|
|
49 |
s.e. 0.0896 0.0731
|
|
|
50 |
|
|
|
51 |
sigma^2 estimated as 0.0002543: log likelihood = 353.96, aic = -701.92
|
|
|
52 |
> update(fit, method = "CSS")
|
|
|
53 |
|
|
|
54 |
Call:
|
|
|
55 |
arima(x = log10(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0,
|
|
|
56 |
1, 1), period = 12), method = "CSS")
|
|
|
57 |
|
|
|
58 |
Coefficients:
|
|
|
59 |
ma1 sma1
|
|
|
60 |
-0.3772 -0.5724
|
|
|
61 |
s.e. 0.0883 0.0704
|
|
|
62 |
|
|
|
63 |
sigma^2 estimated as 0.0002619: part log likelihood = 354.32
|
|
|
64 |
> update(fit, x = window(log10(AirPassengers), start = 1954))
|
|
|
65 |
|
|
|
66 |
Call:
|
|
|
67 |
arima(x = window(log10(AirPassengers), start = 1954), order = c(0, 1, 1), seasonal = list(order = c(0,
|
|
|
68 |
1, 1), period = 12))
|
|
|
69 |
|
|
|
70 |
Coefficients:
|
|
|
71 |
ma1 sma1
|
|
|
72 |
-0.4797 -0.4460
|
|
|
73 |
s.e. 0.1000 0.1514
|
|
|
74 |
|
|
|
75 |
sigma^2 estimated as 0.0001603: log likelihood = 208.02, aic = -410.04
|
|
|
76 |
> pred <- predict(fit, n.ahead = 24)
|
|
|
77 |
> tl <- pred$pred - 1.96 * pred$se
|
|
|
78 |
> tu <- pred$pred + 1.96 * pred$se
|
|
|
79 |
> ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty = c(1, 2, 2))
|
| 56186 |
murdoch |
80 |
>
|
| 86049 |
smeyer |
81 |
> ## full ML fit is the same if the series is reversed, CSS fit is not
|
|
|
82 |
> ap0 <- rev(log10(AirPassengers))
|
|
|
83 |
> attributes(ap0) <- attributes(AirPassengers)
|
|
|
84 |
> arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
|
|
|
85 |
|
|
|
86 |
Call:
|
|
|
87 |
arima(x = ap0, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
|
|
|
88 |
|
|
|
89 |
Coefficients:
|
|
|
90 |
ma1 sma1
|
|
|
91 |
-0.4018 -0.5569
|
|
|
92 |
s.e. 0.0896 0.0731
|
|
|
93 |
|
|
|
94 |
sigma^2 estimated as 0.0002543: log likelihood = 353.95, aic = -701.91
|
|
|
95 |
> arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),
|
|
|
96 |
+ method = "CSS")
|
|
|
97 |
|
|
|
98 |
Call:
|
|
|
99 |
arima(x = ap0, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),
|
|
|
100 |
method = "CSS")
|
|
|
101 |
|
|
|
102 |
Coefficients:
|
|
|
103 |
ma1 sma1
|
|
|
104 |
-0.4188 -0.5224
|
|
|
105 |
s.e. 0.0875 0.0695
|
|
|
106 |
|
|
|
107 |
sigma^2 estimated as 0.0002673: part log likelihood = 353
|
| 56186 |
murdoch |
108 |
>
|
| 86049 |
smeyer |
109 |
> ## Structural Time Series
|
|
|
110 |
> ap <- log10(AirPassengers) - 2
|
|
|
111 |
> (fit <- StructTS(ap, type = "BSM"))
|
|
|
112 |
|
|
|
113 |
Call:
|
|
|
114 |
StructTS(x = ap, type = "BSM")
|
|
|
115 |
|
|
|
116 |
Variances:
|
|
|
117 |
level slope seas epsilon
|
|
|
118 |
0.0001456 0.0000000 0.0002635 0.0000000
|
|
|
119 |
> par(mfrow = c(1, 2))
|
|
|
120 |
> plot(cbind(ap, fitted(fit)), plot.type = "single")
|
|
|
121 |
> plot(cbind(ap, tsSmooth(fit)), plot.type = "single")
|
|
|
122 |
>
|
|
|
123 |
>
|
|
|
124 |
>
|
|
|
125 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
| 56186 |
murdoch |
126 |
> cleanEx()
|
|
|
127 |
> nameEx("BOD")
|
|
|
128 |
> ### * BOD
|
|
|
129 |
>
|
|
|
130 |
> flush(stderr()); flush(stdout())
|
|
|
131 |
>
|
|
|
132 |
> ### Name: BOD
|
|
|
133 |
> ### Title: Biochemical Oxygen Demand
|
|
|
134 |
> ### Aliases: BOD
|
|
|
135 |
> ### Keywords: datasets
|
|
|
136 |
>
|
|
|
137 |
> ### ** Examples
|
|
|
138 |
>
|
|
|
139 |
> ## Don't show:
|
|
|
140 |
> options(show.nls.convergence=FALSE)
|
|
|
141 |
> old <- options(digits = 5)
|
| 66943 |
maechler |
142 |
> ## End(Don't show)
|
| 56186 |
murdoch |
143 |
> require(stats)
|
|
|
144 |
> # simplest form of fitting a first-order model to these data
|
|
|
145 |
> fm1 <- nls(demand ~ A*(1-exp(-exp(lrc)*Time)), data = BOD,
|
| 79454 |
maechler |
146 |
+ start = c(A = 20, lrc = log(.35)))
|
| 56186 |
murdoch |
147 |
> coef(fm1)
|
|
|
148 |
A lrc
|
|
|
149 |
19.14258 -0.63282
|
|
|
150 |
> fm1
|
|
|
151 |
Nonlinear regression model
|
| 61435 |
ripley |
152 |
model: demand ~ A * (1 - exp(-exp(lrc) * Time))
|
|
|
153 |
data: BOD
|
| 56186 |
murdoch |
154 |
A lrc
|
|
|
155 |
19.143 -0.633
|
|
|
156 |
residual sum-of-squares: 26
|
| 85990 |
smeyer |
157 |
> ## IGNORE_RDIFF_BEGIN
|
| 79454 |
maechler |
158 |
> # using the plinear algorithm (trace o/p differs by platform)
|
|
|
159 |
> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,
|
|
|
160 |
+ start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)
|
|
|
161 |
32.946 (7.09e-01): par = (-1.0498 22.126)
|
|
|
162 |
25.992 (1.29e-02): par = (-0.62572 19.103)
|
|
|
163 |
25.990 (2.13e-04): par = (-0.6327 19.142)
|
|
|
164 |
25.990 (3.55e-06): par = (-0.63282 19.143)
|
|
|
165 |
> ## IGNORE_RDIFF_END
|
| 56186 |
murdoch |
166 |
> # using a self-starting model
|
|
|
167 |
> fm3 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD)
|
|
|
168 |
> summary(fm3)
|
|
|
169 |
|
|
|
170 |
Formula: demand ~ SSasympOrig(Time, A, lrc)
|
|
|
171 |
|
|
|
172 |
Parameters:
|
|
|
173 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
174 |
A 19.143 2.496 7.67 0.0016 **
|
|
|
175 |
lrc -0.633 0.382 -1.65 0.1733
|
|
|
176 |
---
|
| 61435 |
ripley |
177 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
178 |
|
| 62439 |
ripley |
179 |
Residual standard error: 2.55 on 4 degrees of freedom
|
| 56186 |
murdoch |
180 |
|
|
|
181 |
> ## Don't show:
|
|
|
182 |
> options(old)
|
| 66943 |
maechler |
183 |
> ## End(Don't show)
|
| 56186 |
murdoch |
184 |
>
|
|
|
185 |
>
|
|
|
186 |
>
|
|
|
187 |
> cleanEx()
|
|
|
188 |
> nameEx("ChickWeight")
|
|
|
189 |
> ### * ChickWeight
|
|
|
190 |
>
|
|
|
191 |
> flush(stderr()); flush(stdout())
|
|
|
192 |
>
|
|
|
193 |
> ### Name: ChickWeight
|
| 89605 |
smeyer |
194 |
> ### Title: Weight Versus Age of Chicks on Different Diets
|
| 56186 |
murdoch |
195 |
> ### Aliases: ChickWeight
|
|
|
196 |
> ### Keywords: datasets
|
|
|
197 |
>
|
|
|
198 |
> ### ** Examples
|
|
|
199 |
>
|
| 86027 |
smeyer |
200 |
> require(graphics)
|
|
|
201 |
> coplot(weight ~ Time | Chick, data = ChickWeight,
|
|
|
202 |
+ type = "b", show.given = FALSE)
|
| 56186 |
murdoch |
203 |
>
|
| 86027 |
smeyer |
204 |
>
|
|
|
205 |
>
|
| 56186 |
murdoch |
206 |
> cleanEx()
|
|
|
207 |
> nameEx("DNase")
|
|
|
208 |
> ### * DNase
|
|
|
209 |
>
|
|
|
210 |
> flush(stderr()); flush(stdout())
|
|
|
211 |
>
|
|
|
212 |
> ### Name: DNase
|
|
|
213 |
> ### Title: Elisa assay of DNase
|
|
|
214 |
> ### Aliases: DNase
|
|
|
215 |
> ### Keywords: datasets
|
|
|
216 |
>
|
|
|
217 |
> ### ** Examples
|
|
|
218 |
>
|
|
|
219 |
> require(stats); require(graphics)
|
|
|
220 |
> ## Don't show:
|
|
|
221 |
> options(show.nls.convergence=FALSE)
|
| 66943 |
maechler |
222 |
> ## End(Don't show)
|
| 56186 |
murdoch |
223 |
> coplot(density ~ conc | Run, data = DNase,
|
|
|
224 |
+ show.given = FALSE, type = "b")
|
|
|
225 |
> coplot(density ~ log(conc) | Run, data = DNase,
|
|
|
226 |
+ show.given = FALSE, type = "b")
|
|
|
227 |
> ## fit a representative run
|
|
|
228 |
> fm1 <- nls(density ~ SSlogis( log(conc), Asym, xmid, scal ),
|
|
|
229 |
+ data = DNase, subset = Run == 1)
|
|
|
230 |
> ## compare with a four-parameter logistic
|
|
|
231 |
> fm2 <- nls(density ~ SSfpl( log(conc), A, B, xmid, scal ),
|
|
|
232 |
+ data = DNase, subset = Run == 1)
|
|
|
233 |
> summary(fm2)
|
|
|
234 |
|
|
|
235 |
Formula: density ~ SSfpl(log(conc), A, B, xmid, scal)
|
|
|
236 |
|
|
|
237 |
Parameters:
|
|
|
238 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
239 |
A -0.007897 0.017200 -0.459 0.654
|
|
|
240 |
B 2.377239 0.109516 21.707 5.35e-11 ***
|
|
|
241 |
xmid 1.507403 0.102080 14.767 4.65e-09 ***
|
|
|
242 |
scal 1.062579 0.056996 18.643 3.16e-10 ***
|
|
|
243 |
---
|
| 61435 |
ripley |
244 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
245 |
|
| 62439 |
ripley |
246 |
Residual standard error: 0.01981 on 12 degrees of freedom
|
| 56186 |
murdoch |
247 |
|
|
|
248 |
> anova(fm1, fm2)
|
|
|
249 |
Analysis of Variance Table
|
|
|
250 |
|
|
|
251 |
Model 1: density ~ SSlogis(log(conc), Asym, xmid, scal)
|
|
|
252 |
Model 2: density ~ SSfpl(log(conc), A, B, xmid, scal)
|
|
|
253 |
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
|
|
|
254 |
1 13 0.0047896
|
|
|
255 |
2 12 0.0047073 1 8.2314e-05 0.2098 0.6551
|
|
|
256 |
>
|
|
|
257 |
>
|
|
|
258 |
>
|
|
|
259 |
> cleanEx()
|
|
|
260 |
> nameEx("Formaldehyde")
|
|
|
261 |
> ### * Formaldehyde
|
|
|
262 |
>
|
|
|
263 |
> flush(stderr()); flush(stdout())
|
|
|
264 |
>
|
|
|
265 |
> ### Name: Formaldehyde
|
|
|
266 |
> ### Title: Determination of Formaldehyde
|
|
|
267 |
> ### Aliases: Formaldehyde
|
|
|
268 |
> ### Keywords: datasets
|
|
|
269 |
>
|
|
|
270 |
> ### ** Examples
|
|
|
271 |
>
|
|
|
272 |
> require(stats); require(graphics)
|
|
|
273 |
> plot(optden ~ carb, data = Formaldehyde,
|
|
|
274 |
+ xlab = "Carbohydrate (ml)", ylab = "Optical Density",
|
|
|
275 |
+ main = "Formaldehyde data", col = 4, las = 1)
|
|
|
276 |
> abline(fm1 <- lm(optden ~ carb, data = Formaldehyde))
|
|
|
277 |
> summary(fm1)
|
|
|
278 |
|
|
|
279 |
Call:
|
|
|
280 |
lm(formula = optden ~ carb, data = Formaldehyde)
|
|
|
281 |
|
|
|
282 |
Residuals:
|
|
|
283 |
1 2 3 4 5 6
|
|
|
284 |
-0.006714 0.001029 0.002771 0.007143 0.007514 -0.011743
|
|
|
285 |
|
|
|
286 |
Coefficients:
|
|
|
287 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
288 |
(Intercept) 0.005086 0.007834 0.649 0.552
|
|
|
289 |
carb 0.876286 0.013535 64.744 3.41e-07 ***
|
|
|
290 |
---
|
| 61435 |
ripley |
291 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
292 |
|
|
|
293 |
Residual standard error: 0.008649 on 4 degrees of freedom
|
| 61435 |
ripley |
294 |
Multiple R-squared: 0.999, Adjusted R-squared: 0.9988
|
|
|
295 |
F-statistic: 4192 on 1 and 4 DF, p-value: 3.409e-07
|
| 56186 |
murdoch |
296 |
|
| 61169 |
ripley |
297 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
298 |
> plot(fm1)
|
|
|
299 |
> par(opar)
|
|
|
300 |
>
|
|
|
301 |
>
|
|
|
302 |
>
|
|
|
303 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
304 |
> cleanEx()
|
|
|
305 |
> nameEx("HairEyeColor")
|
|
|
306 |
> ### * HairEyeColor
|
|
|
307 |
>
|
|
|
308 |
> flush(stderr()); flush(stdout())
|
|
|
309 |
>
|
|
|
310 |
> ### Name: HairEyeColor
|
|
|
311 |
> ### Title: Hair and Eye Color of Statistics Students
|
|
|
312 |
> ### Aliases: HairEyeColor
|
|
|
313 |
> ### Keywords: datasets
|
|
|
314 |
>
|
|
|
315 |
> ### ** Examples
|
|
|
316 |
>
|
|
|
317 |
> require(graphics)
|
|
|
318 |
> ## Full mosaic
|
|
|
319 |
> mosaicplot(HairEyeColor)
|
|
|
320 |
> ## Aggregate over sex (as in Snee's original data)
|
|
|
321 |
> x <- apply(HairEyeColor, c(1, 2), sum)
|
|
|
322 |
> x
|
|
|
323 |
Eye
|
|
|
324 |
Hair Brown Blue Hazel Green
|
|
|
325 |
Black 68 20 15 5
|
|
|
326 |
Brown 119 84 54 29
|
|
|
327 |
Red 26 17 14 14
|
|
|
328 |
Blond 7 94 10 16
|
|
|
329 |
> mosaicplot(x, main = "Relation between hair and eye color")
|
|
|
330 |
>
|
|
|
331 |
>
|
|
|
332 |
>
|
|
|
333 |
> cleanEx()
|
|
|
334 |
> nameEx("Harman23.cor")
|
|
|
335 |
> ### * Harman23.cor
|
|
|
336 |
>
|
|
|
337 |
> flush(stderr()); flush(stdout())
|
|
|
338 |
>
|
|
|
339 |
> ### Name: Harman23.cor
|
|
|
340 |
> ### Title: Harman Example 2.3
|
|
|
341 |
> ### Aliases: Harman23.cor
|
|
|
342 |
> ### Keywords: datasets
|
|
|
343 |
>
|
|
|
344 |
> ### ** Examples
|
|
|
345 |
>
|
|
|
346 |
> require(stats)
|
|
|
347 |
> (Harman23.FA <- factanal(factors = 1, covmat = Harman23.cor))
|
|
|
348 |
|
|
|
349 |
Call:
|
|
|
350 |
factanal(factors = 1, covmat = Harman23.cor)
|
|
|
351 |
|
|
|
352 |
Uniquenesses:
|
|
|
353 |
height arm.span forearm lower.leg weight
|
|
|
354 |
0.158 0.135 0.190 0.187 0.760
|
|
|
355 |
bitro.diameter chest.girth chest.width
|
|
|
356 |
0.829 0.877 0.801
|
|
|
357 |
|
|
|
358 |
Loadings:
|
|
|
359 |
Factor1
|
|
|
360 |
height 0.918
|
|
|
361 |
arm.span 0.930
|
|
|
362 |
forearm 0.900
|
|
|
363 |
lower.leg 0.902
|
|
|
364 |
weight 0.490
|
|
|
365 |
bitro.diameter 0.413
|
|
|
366 |
chest.girth 0.351
|
|
|
367 |
chest.width 0.446
|
|
|
368 |
|
|
|
369 |
Factor1
|
|
|
370 |
SS loadings 4.064
|
|
|
371 |
Proportion Var 0.508
|
|
|
372 |
|
|
|
373 |
Test of the hypothesis that 1 factor is sufficient.
|
|
|
374 |
The chi square statistic is 611.44 on 20 degrees of freedom.
|
|
|
375 |
The p-value is 1.12e-116
|
|
|
376 |
> for(factors in 2:4) print(update(Harman23.FA, factors = factors))
|
|
|
377 |
|
|
|
378 |
Call:
|
|
|
379 |
factanal(factors = factors, covmat = Harman23.cor)
|
|
|
380 |
|
|
|
381 |
Uniquenesses:
|
|
|
382 |
height arm.span forearm lower.leg weight
|
|
|
383 |
0.170 0.107 0.166 0.199 0.089
|
|
|
384 |
bitro.diameter chest.girth chest.width
|
|
|
385 |
0.364 0.416 0.537
|
|
|
386 |
|
|
|
387 |
Loadings:
|
|
|
388 |
Factor1 Factor2
|
|
|
389 |
height 0.865 0.287
|
|
|
390 |
arm.span 0.927 0.181
|
|
|
391 |
forearm 0.895 0.179
|
|
|
392 |
lower.leg 0.859 0.252
|
|
|
393 |
weight 0.233 0.925
|
|
|
394 |
bitro.diameter 0.194 0.774
|
|
|
395 |
chest.girth 0.134 0.752
|
|
|
396 |
chest.width 0.278 0.621
|
|
|
397 |
|
|
|
398 |
Factor1 Factor2
|
|
|
399 |
SS loadings 3.335 2.617
|
|
|
400 |
Proportion Var 0.417 0.327
|
|
|
401 |
Cumulative Var 0.417 0.744
|
|
|
402 |
|
|
|
403 |
Test of the hypothesis that 2 factors are sufficient.
|
|
|
404 |
The chi square statistic is 75.74 on 13 degrees of freedom.
|
|
|
405 |
The p-value is 6.94e-11
|
|
|
406 |
|
|
|
407 |
Call:
|
|
|
408 |
factanal(factors = factors, covmat = Harman23.cor)
|
|
|
409 |
|
|
|
410 |
Uniquenesses:
|
|
|
411 |
height arm.span forearm lower.leg weight
|
|
|
412 |
0.127 0.005 0.193 0.157 0.090
|
|
|
413 |
bitro.diameter chest.girth chest.width
|
|
|
414 |
0.359 0.411 0.490
|
|
|
415 |
|
|
|
416 |
Loadings:
|
|
|
417 |
Factor1 Factor2 Factor3
|
|
|
418 |
height 0.886 0.267 -0.130
|
|
|
419 |
arm.span 0.937 0.195 0.280
|
|
|
420 |
forearm 0.874 0.188
|
|
|
421 |
lower.leg 0.877 0.230 -0.145
|
|
|
422 |
weight 0.242 0.916 -0.106
|
|
|
423 |
bitro.diameter 0.193 0.777
|
|
|
424 |
chest.girth 0.137 0.755
|
|
|
425 |
chest.width 0.261 0.646 0.159
|
|
|
426 |
|
|
|
427 |
Factor1 Factor2 Factor3
|
|
|
428 |
SS loadings 3.379 2.628 0.162
|
|
|
429 |
Proportion Var 0.422 0.329 0.020
|
|
|
430 |
Cumulative Var 0.422 0.751 0.771
|
|
|
431 |
|
|
|
432 |
Test of the hypothesis that 3 factors are sufficient.
|
|
|
433 |
The chi square statistic is 22.81 on 7 degrees of freedom.
|
|
|
434 |
The p-value is 0.00184
|
|
|
435 |
|
|
|
436 |
Call:
|
|
|
437 |
factanal(factors = factors, covmat = Harman23.cor)
|
|
|
438 |
|
|
|
439 |
Uniquenesses:
|
|
|
440 |
height arm.span forearm lower.leg weight
|
|
|
441 |
0.137 0.005 0.191 0.116 0.138
|
|
|
442 |
bitro.diameter chest.girth chest.width
|
|
|
443 |
0.283 0.178 0.488
|
|
|
444 |
|
|
|
445 |
Loadings:
|
|
|
446 |
Factor1 Factor2 Factor3 Factor4
|
|
|
447 |
height 0.879 0.277 -0.115
|
|
|
448 |
arm.span 0.937 0.194 0.277
|
|
|
449 |
forearm 0.875 0.191
|
|
|
450 |
lower.leg 0.887 0.209 0.135 -0.188
|
|
|
451 |
weight 0.246 0.882 0.111 -0.109
|
|
|
452 |
bitro.diameter 0.187 0.822
|
|
|
453 |
chest.girth 0.117 0.729 0.526
|
|
|
454 |
chest.width 0.263 0.644 0.141
|
|
|
455 |
|
|
|
456 |
Factor1 Factor2 Factor3 Factor4
|
|
|
457 |
SS loadings 3.382 2.595 0.323 0.165
|
|
|
458 |
Proportion Var 0.423 0.324 0.040 0.021
|
|
|
459 |
Cumulative Var 0.423 0.747 0.787 0.808
|
|
|
460 |
|
|
|
461 |
Test of the hypothesis that 4 factors are sufficient.
|
|
|
462 |
The chi square statistic is 4.63 on 2 degrees of freedom.
|
|
|
463 |
The p-value is 0.0988
|
|
|
464 |
>
|
|
|
465 |
>
|
|
|
466 |
>
|
|
|
467 |
> cleanEx()
|
|
|
468 |
> nameEx("Harman74.cor")
|
|
|
469 |
> ### * Harman74.cor
|
|
|
470 |
>
|
|
|
471 |
> flush(stderr()); flush(stdout())
|
|
|
472 |
>
|
|
|
473 |
> ### Name: Harman74.cor
|
|
|
474 |
> ### Title: Harman Example 7.4
|
|
|
475 |
> ### Aliases: Harman74.cor
|
|
|
476 |
> ### Keywords: datasets
|
|
|
477 |
>
|
|
|
478 |
> ### ** Examples
|
|
|
479 |
>
|
|
|
480 |
> require(stats)
|
|
|
481 |
> (Harman74.FA <- factanal(factors = 1, covmat = Harman74.cor))
|
|
|
482 |
|
|
|
483 |
Call:
|
|
|
484 |
factanal(factors = 1, covmat = Harman74.cor)
|
|
|
485 |
|
|
|
486 |
Uniquenesses:
|
|
|
487 |
VisualPerception Cubes PaperFormBoard
|
|
|
488 |
0.677 0.866 0.830
|
|
|
489 |
Flags GeneralInformation PargraphComprehension
|
|
|
490 |
0.768 0.487 0.491
|
|
|
491 |
SentenceCompletion WordClassification WordMeaning
|
|
|
492 |
0.500 0.514 0.474
|
|
|
493 |
Addition Code CountingDots
|
|
|
494 |
0.818 0.731 0.824
|
|
|
495 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
496 |
0.681 0.833 0.863
|
|
|
497 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
498 |
0.775 0.812 0.778
|
|
|
499 |
FigureWord Deduction NumericalPuzzles
|
|
|
500 |
0.816 0.612 0.676
|
|
|
501 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
502 |
0.619 0.524 0.593
|
|
|
503 |
|
|
|
504 |
Loadings:
|
|
|
505 |
Factor1
|
|
|
506 |
VisualPerception 0.569
|
|
|
507 |
Cubes 0.366
|
|
|
508 |
PaperFormBoard 0.412
|
|
|
509 |
Flags 0.482
|
|
|
510 |
GeneralInformation 0.716
|
|
|
511 |
PargraphComprehension 0.713
|
|
|
512 |
SentenceCompletion 0.707
|
|
|
513 |
WordClassification 0.697
|
|
|
514 |
WordMeaning 0.725
|
|
|
515 |
Addition 0.426
|
|
|
516 |
Code 0.519
|
|
|
517 |
CountingDots 0.419
|
|
|
518 |
StraightCurvedCapitals 0.565
|
|
|
519 |
WordRecognition 0.408
|
|
|
520 |
NumberRecognition 0.370
|
|
|
521 |
FigureRecognition 0.474
|
|
|
522 |
ObjectNumber 0.434
|
|
|
523 |
NumberFigure 0.471
|
|
|
524 |
FigureWord 0.429
|
|
|
525 |
Deduction 0.623
|
|
|
526 |
NumericalPuzzles 0.569
|
|
|
527 |
ProblemReasoning 0.617
|
|
|
528 |
SeriesCompletion 0.690
|
|
|
529 |
ArithmeticProblems 0.638
|
|
|
530 |
|
|
|
531 |
Factor1
|
|
|
532 |
SS loadings 7.438
|
|
|
533 |
Proportion Var 0.310
|
|
|
534 |
|
|
|
535 |
Test of the hypothesis that 1 factor is sufficient.
|
|
|
536 |
The chi square statistic is 622.91 on 252 degrees of freedom.
|
|
|
537 |
The p-value is 2.28e-33
|
|
|
538 |
> for(factors in 2:5) print(update(Harman74.FA, factors = factors))
|
|
|
539 |
|
|
|
540 |
Call:
|
|
|
541 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
542 |
|
|
|
543 |
Uniquenesses:
|
|
|
544 |
VisualPerception Cubes PaperFormBoard
|
|
|
545 |
0.650 0.864 0.844
|
|
|
546 |
Flags GeneralInformation PargraphComprehension
|
|
|
547 |
0.778 0.375 0.316
|
|
|
548 |
SentenceCompletion WordClassification WordMeaning
|
|
|
549 |
0.319 0.503 0.258
|
|
|
550 |
Addition Code CountingDots
|
|
|
551 |
0.670 0.608 0.581
|
|
|
552 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
553 |
0.567 0.832 0.850
|
|
|
554 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
555 |
0.743 0.770 0.625
|
|
|
556 |
FigureWord Deduction NumericalPuzzles
|
|
|
557 |
0.792 0.629 0.579
|
|
|
558 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
559 |
0.634 0.539 0.553
|
|
|
560 |
|
|
|
561 |
Loadings:
|
|
|
562 |
Factor1 Factor2
|
|
|
563 |
VisualPerception 0.506 0.306
|
|
|
564 |
Cubes 0.304 0.209
|
|
|
565 |
PaperFormBoard 0.297 0.260
|
|
|
566 |
Flags 0.327 0.339
|
|
|
567 |
GeneralInformation 0.240 0.753
|
|
|
568 |
PargraphComprehension 0.171 0.809
|
|
|
569 |
SentenceCompletion 0.163 0.809
|
|
|
570 |
WordClassification 0.344 0.615
|
|
|
571 |
WordMeaning 0.148 0.849
|
|
|
572 |
Addition 0.563 0.115
|
|
|
573 |
Code 0.591 0.207
|
|
|
574 |
CountingDots 0.647
|
|
|
575 |
StraightCurvedCapitals 0.612 0.241
|
|
|
576 |
WordRecognition 0.315 0.263
|
|
|
577 |
NumberRecognition 0.328 0.205
|
|
|
578 |
FigureRecognition 0.457 0.218
|
|
|
579 |
ObjectNumber 0.431 0.209
|
|
|
580 |
NumberFigure 0.601 0.116
|
|
|
581 |
FigureWord 0.399 0.222
|
|
|
582 |
Deduction 0.379 0.477
|
|
|
583 |
NumericalPuzzles 0.604 0.237
|
|
|
584 |
ProblemReasoning 0.390 0.462
|
|
|
585 |
SeriesCompletion 0.486 0.474
|
|
|
586 |
ArithmeticProblems 0.544 0.389
|
|
|
587 |
|
|
|
588 |
Factor1 Factor2
|
|
|
589 |
SS loadings 4.573 4.548
|
|
|
590 |
Proportion Var 0.191 0.190
|
|
|
591 |
Cumulative Var 0.191 0.380
|
|
|
592 |
|
|
|
593 |
Test of the hypothesis that 2 factors are sufficient.
|
|
|
594 |
The chi square statistic is 420.24 on 229 degrees of freedom.
|
|
|
595 |
The p-value is 2.01e-13
|
|
|
596 |
|
|
|
597 |
Call:
|
|
|
598 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
599 |
|
|
|
600 |
Uniquenesses:
|
|
|
601 |
VisualPerception Cubes PaperFormBoard
|
|
|
602 |
0.500 0.793 0.662
|
|
|
603 |
Flags GeneralInformation PargraphComprehension
|
|
|
604 |
0.694 0.352 0.316
|
|
|
605 |
SentenceCompletion WordClassification WordMeaning
|
|
|
606 |
0.300 0.502 0.256
|
|
|
607 |
Addition Code CountingDots
|
|
|
608 |
0.200 0.586 0.494
|
|
|
609 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
610 |
0.569 0.838 0.848
|
|
|
611 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
612 |
0.643 0.780 0.635
|
|
|
613 |
FigureWord Deduction NumericalPuzzles
|
|
|
614 |
0.788 0.590 0.580
|
|
|
615 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
616 |
0.597 0.498 0.500
|
|
|
617 |
|
|
|
618 |
Loadings:
|
|
|
619 |
Factor1 Factor2 Factor3
|
|
|
620 |
VisualPerception 0.176 0.656 0.198
|
|
|
621 |
Cubes 0.122 0.428
|
|
|
622 |
PaperFormBoard 0.145 0.563
|
|
|
623 |
Flags 0.239 0.487 0.107
|
|
|
624 |
GeneralInformation 0.745 0.191 0.237
|
|
|
625 |
PargraphComprehension 0.780 0.249 0.118
|
|
|
626 |
SentenceCompletion 0.802 0.175 0.160
|
|
|
627 |
WordClassification 0.571 0.327 0.256
|
|
|
628 |
WordMeaning 0.821 0.248
|
|
|
629 |
Addition 0.162 -0.118 0.871
|
|
|
630 |
Code 0.198 0.219 0.572
|
|
|
631 |
CountingDots 0.179 0.688
|
|
|
632 |
StraightCurvedCapitals 0.190 0.381 0.499
|
|
|
633 |
WordRecognition 0.231 0.253 0.210
|
|
|
634 |
NumberRecognition 0.158 0.299 0.195
|
|
|
635 |
FigureRecognition 0.108 0.557 0.186
|
|
|
636 |
ObjectNumber 0.178 0.267 0.342
|
|
|
637 |
NumberFigure 0.427 0.424
|
|
|
638 |
FigureWord 0.167 0.355 0.240
|
|
|
639 |
Deduction 0.392 0.472 0.181
|
|
|
640 |
NumericalPuzzles 0.178 0.406 0.473
|
|
|
641 |
ProblemReasoning 0.382 0.473 0.182
|
|
|
642 |
SeriesCompletion 0.379 0.528 0.283
|
|
|
643 |
ArithmeticProblems 0.377 0.226 0.554
|
|
|
644 |
|
|
|
645 |
Factor1 Factor2 Factor3
|
|
|
646 |
SS loadings 3.802 3.488 3.186
|
|
|
647 |
Proportion Var 0.158 0.145 0.133
|
|
|
648 |
Cumulative Var 0.158 0.304 0.436
|
|
|
649 |
|
|
|
650 |
Test of the hypothesis that 3 factors are sufficient.
|
|
|
651 |
The chi square statistic is 295.59 on 207 degrees of freedom.
|
|
|
652 |
The p-value is 5.12e-05
|
|
|
653 |
|
|
|
654 |
Call:
|
|
|
655 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
656 |
|
|
|
657 |
Uniquenesses:
|
|
|
658 |
VisualPerception Cubes PaperFormBoard
|
|
|
659 |
0.438 0.780 0.644
|
|
|
660 |
Flags GeneralInformation PargraphComprehension
|
|
|
661 |
0.651 0.352 0.312
|
|
|
662 |
SentenceCompletion WordClassification WordMeaning
|
|
|
663 |
0.283 0.485 0.257
|
|
|
664 |
Addition Code CountingDots
|
|
|
665 |
0.240 0.551 0.435
|
|
|
666 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
667 |
0.491 0.646 0.696
|
|
|
668 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
669 |
0.549 0.598 0.593
|
|
|
670 |
FigureWord Deduction NumericalPuzzles
|
|
|
671 |
0.762 0.592 0.583
|
|
|
672 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
673 |
0.601 0.497 0.500
|
|
|
674 |
|
|
|
675 |
Loadings:
|
|
|
676 |
Factor1 Factor2 Factor3 Factor4
|
|
|
677 |
VisualPerception 0.160 0.689 0.187 0.160
|
|
|
678 |
Cubes 0.117 0.436
|
|
|
679 |
PaperFormBoard 0.137 0.570 0.110
|
|
|
680 |
Flags 0.233 0.527
|
|
|
681 |
GeneralInformation 0.739 0.185 0.213 0.150
|
|
|
682 |
PargraphComprehension 0.767 0.205 0.233
|
|
|
683 |
SentenceCompletion 0.806 0.197 0.153
|
|
|
684 |
WordClassification 0.569 0.339 0.242 0.132
|
|
|
685 |
WordMeaning 0.806 0.201 0.227
|
|
|
686 |
Addition 0.167 -0.118 0.831 0.166
|
|
|
687 |
Code 0.180 0.120 0.512 0.374
|
|
|
688 |
CountingDots 0.210 0.716
|
|
|
689 |
StraightCurvedCapitals 0.188 0.438 0.525
|
|
|
690 |
WordRecognition 0.197 0.553
|
|
|
691 |
NumberRecognition 0.122 0.116 0.520
|
|
|
692 |
FigureRecognition 0.408 0.525
|
|
|
693 |
ObjectNumber 0.142 0.219 0.574
|
|
|
694 |
NumberFigure 0.293 0.336 0.456
|
|
|
695 |
FigureWord 0.148 0.239 0.161 0.365
|
|
|
696 |
Deduction 0.378 0.402 0.118 0.301
|
|
|
697 |
NumericalPuzzles 0.175 0.381 0.438 0.223
|
|
|
698 |
ProblemReasoning 0.366 0.399 0.123 0.301
|
|
|
699 |
SeriesCompletion 0.369 0.500 0.244 0.239
|
|
|
700 |
ArithmeticProblems 0.370 0.158 0.496 0.304
|
|
|
701 |
|
|
|
702 |
Factor1 Factor2 Factor3 Factor4
|
|
|
703 |
SS loadings 3.647 2.872 2.657 2.290
|
|
|
704 |
Proportion Var 0.152 0.120 0.111 0.095
|
|
|
705 |
Cumulative Var 0.152 0.272 0.382 0.478
|
|
|
706 |
|
|
|
707 |
Test of the hypothesis that 4 factors are sufficient.
|
|
|
708 |
The chi square statistic is 226.68 on 186 degrees of freedom.
|
|
|
709 |
The p-value is 0.0224
|
|
|
710 |
|
|
|
711 |
Call:
|
|
|
712 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
713 |
|
|
|
714 |
Uniquenesses:
|
|
|
715 |
VisualPerception Cubes PaperFormBoard
|
|
|
716 |
0.450 0.781 0.639
|
|
|
717 |
Flags GeneralInformation PargraphComprehension
|
|
|
718 |
0.649 0.357 0.288
|
|
|
719 |
SentenceCompletion WordClassification WordMeaning
|
|
|
720 |
0.277 0.485 0.262
|
|
|
721 |
Addition Code CountingDots
|
|
|
722 |
0.215 0.386 0.444
|
|
|
723 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
724 |
0.256 0.639 0.706
|
|
|
725 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
726 |
0.550 0.614 0.596
|
|
|
727 |
FigureWord Deduction NumericalPuzzles
|
|
|
728 |
0.764 0.521 0.564
|
|
|
729 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
730 |
0.580 0.442 0.478
|
|
|
731 |
|
|
|
732 |
Loadings:
|
|
|
733 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
734 |
VisualPerception 0.161 0.658 0.136 0.182 0.199
|
|
|
735 |
Cubes 0.113 0.435 0.107
|
|
|
736 |
PaperFormBoard 0.135 0.562 0.107 0.116
|
|
|
737 |
Flags 0.231 0.533
|
|
|
738 |
GeneralInformation 0.736 0.188 0.192 0.162
|
|
|
739 |
PargraphComprehension 0.775 0.187 0.251 0.113
|
|
|
740 |
SentenceCompletion 0.809 0.208 0.136
|
|
|
741 |
WordClassification 0.568 0.348 0.223 0.131
|
|
|
742 |
WordMeaning 0.800 0.215 0.224
|
|
|
743 |
Addition 0.175 -0.100 0.844 0.176
|
|
|
744 |
Code 0.185 0.438 0.451 0.426
|
|
|
745 |
CountingDots 0.222 0.690 0.101 0.140
|
|
|
746 |
StraightCurvedCapitals 0.186 0.425 0.458 0.559
|
|
|
747 |
WordRecognition 0.197 0.557
|
|
|
748 |
NumberRecognition 0.121 0.130 0.508
|
|
|
749 |
FigureRecognition 0.400 0.529
|
|
|
750 |
ObjectNumber 0.145 0.208 0.562
|
|
|
751 |
NumberFigure 0.306 0.325 0.452
|
|
|
752 |
FigureWord 0.147 0.242 0.145 0.364
|
|
|
753 |
Deduction 0.370 0.452 0.139 0.287 -0.190
|
|
|
754 |
NumericalPuzzles 0.170 0.402 0.439 0.230
|
|
|
755 |
ProblemReasoning 0.358 0.423 0.126 0.302
|
|
|
756 |
SeriesCompletion 0.360 0.549 0.256 0.223 -0.107
|
|
|
757 |
ArithmeticProblems 0.371 0.185 0.502 0.307
|
|
|
758 |
|
|
|
759 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
760 |
SS loadings 3.632 2.964 2.456 2.345 0.663
|
|
|
761 |
Proportion Var 0.151 0.124 0.102 0.098 0.028
|
|
|
762 |
Cumulative Var 0.151 0.275 0.377 0.475 0.503
|
|
|
763 |
|
|
|
764 |
Test of the hypothesis that 5 factors are sufficient.
|
|
|
765 |
The chi square statistic is 186.82 on 166 degrees of freedom.
|
|
|
766 |
The p-value is 0.128
|
|
|
767 |
> Harman74.FA <- factanal(factors = 5, covmat = Harman74.cor,
|
| 61157 |
ripley |
768 |
+ rotation = "promax")
|
| 56186 |
murdoch |
769 |
> print(Harman74.FA$loadings, sort = TRUE)
|
|
|
770 |
|
|
|
771 |
Loadings:
|
|
|
772 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
773 |
VisualPerception 0.831 -0.127 0.230
|
|
|
774 |
Cubes 0.534
|
|
|
775 |
PaperFormBoard 0.736 -0.290 0.136
|
|
|
776 |
Flags 0.647 -0.104
|
|
|
777 |
SeriesCompletion 0.555 0.126 0.127
|
|
|
778 |
GeneralInformation 0.764
|
|
|
779 |
PargraphComprehension 0.845 -0.140 0.140
|
|
|
780 |
SentenceCompletion 0.872 -0.140
|
|
|
781 |
WordClassification 0.277 0.505 0.104
|
|
|
782 |
WordMeaning 0.846 -0.108
|
|
|
783 |
Addition -0.334 1.012
|
|
|
784 |
CountingDots 0.206 -0.200 0.722 0.185
|
|
|
785 |
ArithmeticProblems 0.197 0.500 0.139
|
|
|
786 |
WordRecognition -0.126 0.127 -0.103 0.657
|
|
|
787 |
NumberRecognition 0.568
|
|
|
788 |
FigureRecognition 0.399 -0.142 -0.207 0.562
|
|
|
789 |
ObjectNumber -0.108 0.107 0.613
|
|
|
790 |
StraightCurvedCapitals 0.542 0.247 0.618
|
|
|
791 |
Code 0.112 0.288 0.486 0.424
|
|
|
792 |
NumberFigure 0.255 -0.230 0.211 0.413
|
|
|
793 |
FigureWord 0.187 0.347
|
|
|
794 |
Deduction 0.404 0.169 0.117 -0.203
|
|
|
795 |
NumericalPuzzles 0.393 0.368
|
|
|
796 |
ProblemReasoning 0.381 0.188 0.169
|
|
|
797 |
|
|
|
798 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
799 |
SS loadings 3.529 3.311 2.367 2.109 0.762
|
|
|
800 |
Proportion Var 0.147 0.138 0.099 0.088 0.032
|
|
|
801 |
Cumulative Var 0.147 0.285 0.384 0.471 0.503
|
|
|
802 |
>
|
|
|
803 |
>
|
|
|
804 |
>
|
|
|
805 |
> cleanEx()
|
|
|
806 |
> nameEx("InsectSprays")
|
|
|
807 |
> ### * InsectSprays
|
|
|
808 |
>
|
|
|
809 |
> flush(stderr()); flush(stdout())
|
|
|
810 |
>
|
|
|
811 |
> ### Name: InsectSprays
|
|
|
812 |
> ### Title: Effectiveness of Insect Sprays
|
|
|
813 |
> ### Aliases: InsectSprays
|
|
|
814 |
> ### Keywords: datasets
|
|
|
815 |
>
|
|
|
816 |
> ### ** Examples
|
|
|
817 |
>
|
|
|
818 |
> require(stats); require(graphics)
|
|
|
819 |
> boxplot(count ~ spray, data = InsectSprays,
|
|
|
820 |
+ xlab = "Type of spray", ylab = "Insect count",
|
|
|
821 |
+ main = "InsectSprays data", varwidth = TRUE, col = "lightgray")
|
|
|
822 |
> fm1 <- aov(count ~ spray, data = InsectSprays)
|
|
|
823 |
> summary(fm1)
|
| 57044 |
ripley |
824 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
825 |
spray 5 2669 533.8 34.7 <2e-16 ***
|
|
|
826 |
Residuals 66 1015 15.4
|
| 56186 |
murdoch |
827 |
---
|
| 61435 |
ripley |
828 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 61169 |
ripley |
829 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
830 |
> plot(fm1)
|
|
|
831 |
> fm2 <- aov(sqrt(count) ~ spray, data = InsectSprays)
|
|
|
832 |
> summary(fm2)
|
| 57044 |
ripley |
833 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
834 |
spray 5 88.44 17.688 44.8 <2e-16 ***
|
|
|
835 |
Residuals 66 26.06 0.395
|
| 56186 |
murdoch |
836 |
---
|
| 61435 |
ripley |
837 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
838 |
> plot(fm2)
|
|
|
839 |
> par(opar)
|
|
|
840 |
>
|
|
|
841 |
>
|
|
|
842 |
>
|
|
|
843 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
844 |
> cleanEx()
|
|
|
845 |
> nameEx("JohnsonJohnson")
|
|
|
846 |
> ### * JohnsonJohnson
|
|
|
847 |
>
|
|
|
848 |
> flush(stderr()); flush(stdout())
|
|
|
849 |
>
|
|
|
850 |
> ### Name: JohnsonJohnson
|
|
|
851 |
> ### Title: Quarterly Earnings per Johnson & Johnson Share
|
|
|
852 |
> ### Aliases: JohnsonJohnson
|
|
|
853 |
> ### Keywords: datasets
|
|
|
854 |
>
|
|
|
855 |
> ### ** Examples
|
|
|
856 |
>
|
|
|
857 |
>
|
|
|
858 |
> cleanEx()
|
|
|
859 |
> nameEx("LifeCycleSavings")
|
|
|
860 |
> ### * LifeCycleSavings
|
|
|
861 |
>
|
|
|
862 |
> flush(stderr()); flush(stdout())
|
|
|
863 |
>
|
|
|
864 |
> ### Name: LifeCycleSavings
|
|
|
865 |
> ### Title: Intercountry Life-Cycle Savings Data
|
|
|
866 |
> ### Aliases: LifeCycleSavings
|
|
|
867 |
> ### Keywords: datasets
|
|
|
868 |
>
|
|
|
869 |
> ### ** Examples
|
|
|
870 |
>
|
|
|
871 |
> require(stats); require(graphics)
|
|
|
872 |
> pairs(LifeCycleSavings, panel = panel.smooth,
|
|
|
873 |
+ main = "LifeCycleSavings data")
|
|
|
874 |
> fm1 <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
|
|
|
875 |
> summary(fm1)
|
|
|
876 |
|
|
|
877 |
Call:
|
|
|
878 |
lm(formula = sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
|
|
|
879 |
|
|
|
880 |
Residuals:
|
|
|
881 |
Min 1Q Median 3Q Max
|
|
|
882 |
-8.2422 -2.6857 -0.2488 2.4280 9.7509
|
|
|
883 |
|
|
|
884 |
Coefficients:
|
|
|
885 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
886 |
(Intercept) 28.5660865 7.3545161 3.884 0.000334 ***
|
|
|
887 |
pop15 -0.4611931 0.1446422 -3.189 0.002603 **
|
|
|
888 |
pop75 -1.6914977 1.0835989 -1.561 0.125530
|
|
|
889 |
dpi -0.0003369 0.0009311 -0.362 0.719173
|
|
|
890 |
ddpi 0.4096949 0.1961971 2.088 0.042471 *
|
|
|
891 |
---
|
| 61435 |
ripley |
892 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
893 |
|
|
|
894 |
Residual standard error: 3.803 on 45 degrees of freedom
|
| 61435 |
ripley |
895 |
Multiple R-squared: 0.3385, Adjusted R-squared: 0.2797
|
|
|
896 |
F-statistic: 5.756 on 4 and 45 DF, p-value: 0.0007904
|
| 56186 |
murdoch |
897 |
|
|
|
898 |
>
|
| 88216 |
smeyer |
899 |
> ## recreate pre-4.6 dataset with a misspelled row name
|
|
|
900 |
> idx <- which(row.names(LifeCycleSavings) == "Guatemala")
|
|
|
901 |
> stopifnot(length(idx) == 1)
|
|
|
902 |
> LCS <- LifeCycleSavings; row.names(LCS)[idx] <- "Guatamala" # sic!
|
| 56186 |
murdoch |
903 |
>
|
|
|
904 |
>
|
| 88216 |
smeyer |
905 |
>
|
| 56186 |
murdoch |
906 |
> cleanEx()
|
|
|
907 |
> nameEx("Loblolly")
|
|
|
908 |
> ### * Loblolly
|
|
|
909 |
>
|
|
|
910 |
> flush(stderr()); flush(stdout())
|
|
|
911 |
>
|
|
|
912 |
> ### Name: Loblolly
|
| 85956 |
ripley |
913 |
> ### Title: Growth of Loblolly Pine Trees
|
| 56186 |
murdoch |
914 |
> ### Aliases: Loblolly
|
|
|
915 |
> ### Keywords: datasets
|
|
|
916 |
>
|
|
|
917 |
> ### ** Examples
|
|
|
918 |
>
|
|
|
919 |
> require(stats); require(graphics)
|
|
|
920 |
> plot(height ~ age, data = Loblolly, subset = Seed == 329,
|
|
|
921 |
+ xlab = "Tree age (yr)", las = 1,
|
|
|
922 |
+ ylab = "Tree height (ft)",
|
|
|
923 |
+ main = "Loblolly data and fitted curve (Seed 329 only)")
|
|
|
924 |
> fm1 <- nls(height ~ SSasymp(age, Asym, R0, lrc),
|
|
|
925 |
+ data = Loblolly, subset = Seed == 329)
|
|
|
926 |
> age <- seq(0, 30, length.out = 101)
|
|
|
927 |
> lines(age, predict(fm1, list(age = age)))
|
|
|
928 |
>
|
|
|
929 |
>
|
|
|
930 |
>
|
|
|
931 |
> cleanEx()
|
|
|
932 |
> nameEx("Nile")
|
|
|
933 |
> ### * Nile
|
|
|
934 |
>
|
|
|
935 |
> flush(stderr()); flush(stdout())
|
|
|
936 |
>
|
|
|
937 |
> ### Name: Nile
|
|
|
938 |
> ### Title: Flow of the River Nile
|
|
|
939 |
> ### Aliases: Nile
|
|
|
940 |
> ### Keywords: datasets
|
|
|
941 |
>
|
|
|
942 |
> ### ** Examples
|
|
|
943 |
>
|
|
|
944 |
> require(stats); require(graphics)
|
| 61169 |
ripley |
945 |
> par(mfrow = c(2, 2))
|
| 56186 |
murdoch |
946 |
> plot(Nile)
|
|
|
947 |
> acf(Nile)
|
|
|
948 |
> pacf(Nile)
|
|
|
949 |
> ar(Nile) # selects order 2
|
|
|
950 |
|
|
|
951 |
Call:
|
|
|
952 |
ar(x = Nile)
|
|
|
953 |
|
|
|
954 |
Coefficients:
|
|
|
955 |
1 2
|
|
|
956 |
0.4081 0.1812
|
|
|
957 |
|
| 61435 |
ripley |
958 |
Order selected 2 sigma^2 estimated as 21247
|
| 56186 |
murdoch |
959 |
> cpgram(ar(Nile)$resid)
|
| 61169 |
ripley |
960 |
> par(mfrow = c(1, 1))
|
| 56186 |
murdoch |
961 |
> arima(Nile, c(2, 0, 0))
|
|
|
962 |
|
|
|
963 |
Call:
|
|
|
964 |
arima(x = Nile, order = c(2, 0, 0))
|
|
|
965 |
|
|
|
966 |
Coefficients:
|
|
|
967 |
ar1 ar2 intercept
|
|
|
968 |
0.4096 0.1987 919.8397
|
|
|
969 |
s.e. 0.0974 0.0990 35.6410
|
|
|
970 |
|
|
|
971 |
sigma^2 estimated as 20291: log likelihood = -637.98, aic = 1283.96
|
|
|
972 |
>
|
|
|
973 |
> ## Now consider missing values, following Durbin & Koopman
|
|
|
974 |
> NileNA <- Nile
|
|
|
975 |
> NileNA[c(21:40, 61:80)] <- NA
|
|
|
976 |
> arima(NileNA, c(2, 0, 0))
|
|
|
977 |
|
|
|
978 |
Call:
|
|
|
979 |
arima(x = NileNA, order = c(2, 0, 0))
|
|
|
980 |
|
|
|
981 |
Coefficients:
|
|
|
982 |
ar1 ar2 intercept
|
|
|
983 |
0.3622 0.1678 918.3103
|
|
|
984 |
s.e. 0.1273 0.1323 39.5037
|
|
|
985 |
|
|
|
986 |
sigma^2 estimated as 23676: log likelihood = -387.7, aic = 783.41
|
|
|
987 |
> plot(NileNA)
|
|
|
988 |
> pred <-
|
| 61169 |
ripley |
989 |
+ predict(arima(window(NileNA, 1871, 1890), c(2, 0, 0)), n.ahead = 20)
|
| 56186 |
murdoch |
990 |
> lines(pred$pred, lty = 3, col = "red")
|
| 61157 |
ripley |
991 |
> lines(pred$pred + 2*pred$se, lty = 2, col = "blue")
|
|
|
992 |
> lines(pred$pred - 2*pred$se, lty = 2, col = "blue")
|
| 56186 |
murdoch |
993 |
> pred <-
|
| 61169 |
ripley |
994 |
+ predict(arima(window(NileNA, 1871, 1930), c(2, 0, 0)), n.ahead = 20)
|
| 56186 |
murdoch |
995 |
> lines(pred$pred, lty = 3, col = "red")
|
| 61157 |
ripley |
996 |
> lines(pred$pred + 2*pred$se, lty = 2, col = "blue")
|
|
|
997 |
> lines(pred$pred - 2*pred$se, lty = 2, col = "blue")
|
| 56186 |
murdoch |
998 |
>
|
|
|
999 |
> ## Structural time series models
|
|
|
1000 |
> par(mfrow = c(3, 1))
|
|
|
1001 |
> plot(Nile)
|
|
|
1002 |
> ## local level model
|
|
|
1003 |
> (fit <- StructTS(Nile, type = "level"))
|
|
|
1004 |
|
|
|
1005 |
Call:
|
|
|
1006 |
StructTS(x = Nile, type = "level")
|
|
|
1007 |
|
|
|
1008 |
Variances:
|
|
|
1009 |
level epsilon
|
|
|
1010 |
1469 15099
|
|
|
1011 |
> lines(fitted(fit), lty = 2) # contemporaneous smoothing
|
|
|
1012 |
> lines(tsSmooth(fit), lty = 2, col = 4) # fixed-interval smoothing
|
|
|
1013 |
> plot(residuals(fit)); abline(h = 0, lty = 3)
|
|
|
1014 |
> ## local trend model
|
|
|
1015 |
> (fit2 <- StructTS(Nile, type = "trend")) ## constant trend fitted
|
|
|
1016 |
|
|
|
1017 |
Call:
|
|
|
1018 |
StructTS(x = Nile, type = "trend")
|
|
|
1019 |
|
|
|
1020 |
Variances:
|
|
|
1021 |
level slope epsilon
|
|
|
1022 |
1427 0 15047
|
|
|
1023 |
> pred <- predict(fit, n.ahead = 30)
|
|
|
1024 |
> ## with 50% confidence interval
|
|
|
1025 |
> ts.plot(Nile, pred$pred,
|
|
|
1026 |
+ pred$pred + 0.67*pred$se, pred$pred -0.67*pred$se)
|
|
|
1027 |
>
|
|
|
1028 |
> ## Now consider missing values
|
|
|
1029 |
> plot(NileNA)
|
|
|
1030 |
> (fit3 <- StructTS(NileNA, type = "level"))
|
|
|
1031 |
|
|
|
1032 |
Call:
|
|
|
1033 |
StructTS(x = NileNA, type = "level")
|
|
|
1034 |
|
|
|
1035 |
Variances:
|
|
|
1036 |
level epsilon
|
|
|
1037 |
685.8 17899.8
|
|
|
1038 |
> lines(fitted(fit3), lty = 2)
|
|
|
1039 |
> lines(tsSmooth(fit3), lty = 3)
|
|
|
1040 |
> plot(residuals(fit3)); abline(h = 0, lty = 3)
|
|
|
1041 |
>
|
|
|
1042 |
>
|
|
|
1043 |
>
|
|
|
1044 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1045 |
> cleanEx()
|
|
|
1046 |
> nameEx("Orange")
|
|
|
1047 |
> ### * Orange
|
|
|
1048 |
>
|
|
|
1049 |
> flush(stderr()); flush(stdout())
|
|
|
1050 |
>
|
|
|
1051 |
> ### Name: Orange
|
|
|
1052 |
> ### Title: Growth of Orange Trees
|
|
|
1053 |
> ### Aliases: Orange
|
|
|
1054 |
> ### Keywords: datasets
|
|
|
1055 |
>
|
|
|
1056 |
> ### ** Examples
|
|
|
1057 |
>
|
|
|
1058 |
> require(stats); require(graphics)
|
|
|
1059 |
> coplot(circumference ~ age | Tree, data = Orange, show.given = FALSE)
|
|
|
1060 |
> fm1 <- nls(circumference ~ SSlogis(age, Asym, xmid, scal),
|
|
|
1061 |
+ data = Orange, subset = Tree == 3)
|
|
|
1062 |
> plot(circumference ~ age, data = Orange, subset = Tree == 3,
|
|
|
1063 |
+ xlab = "Tree age (days since 1968/12/31)",
|
|
|
1064 |
+ ylab = "Tree circumference (mm)", las = 1,
|
|
|
1065 |
+ main = "Orange tree data and fitted model (Tree 3 only)")
|
|
|
1066 |
> age <- seq(0, 1600, length.out = 101)
|
|
|
1067 |
> lines(age, predict(fm1, list(age = age)))
|
|
|
1068 |
>
|
|
|
1069 |
>
|
|
|
1070 |
>
|
|
|
1071 |
> cleanEx()
|
|
|
1072 |
> nameEx("OrchardSprays")
|
|
|
1073 |
> ### * OrchardSprays
|
|
|
1074 |
>
|
|
|
1075 |
> flush(stderr()); flush(stdout())
|
|
|
1076 |
>
|
|
|
1077 |
> ### Name: OrchardSprays
|
|
|
1078 |
> ### Title: Potency of Orchard Sprays
|
|
|
1079 |
> ### Aliases: OrchardSprays
|
|
|
1080 |
> ### Keywords: datasets
|
|
|
1081 |
>
|
|
|
1082 |
> ### ** Examples
|
|
|
1083 |
>
|
|
|
1084 |
> require(graphics)
|
|
|
1085 |
> pairs(OrchardSprays, main = "OrchardSprays data")
|
|
|
1086 |
>
|
|
|
1087 |
>
|
|
|
1088 |
>
|
|
|
1089 |
> cleanEx()
|
|
|
1090 |
> nameEx("PlantGrowth")
|
|
|
1091 |
> ### * PlantGrowth
|
|
|
1092 |
>
|
|
|
1093 |
> flush(stderr()); flush(stdout())
|
|
|
1094 |
>
|
|
|
1095 |
> ### Name: PlantGrowth
|
|
|
1096 |
> ### Title: Results from an Experiment on Plant Growth
|
|
|
1097 |
> ### Aliases: PlantGrowth
|
|
|
1098 |
> ### Keywords: datasets
|
|
|
1099 |
>
|
|
|
1100 |
> ### ** Examples
|
|
|
1101 |
>
|
|
|
1102 |
> ## One factor ANOVA example from Dobson's book, cf. Table 7.4:
|
|
|
1103 |
> require(stats); require(graphics)
|
|
|
1104 |
> boxplot(weight ~ group, data = PlantGrowth, main = "PlantGrowth data",
|
|
|
1105 |
+ ylab = "Dried weight of plants", col = "lightgray",
|
|
|
1106 |
+ notch = TRUE, varwidth = TRUE)
|
| 78669 |
ripley |
1107 |
Warning in (function (z, notch = FALSE, width = NULL, varwidth = FALSE, :
|
| 56186 |
murdoch |
1108 |
some notches went outside hinges ('box'): maybe set notch=FALSE
|
|
|
1109 |
> anova(lm(weight ~ group, data = PlantGrowth))
|
|
|
1110 |
Analysis of Variance Table
|
|
|
1111 |
|
|
|
1112 |
Response: weight
|
|
|
1113 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1114 |
group 2 3.7663 1.8832 4.8461 0.01591 *
|
|
|
1115 |
Residuals 27 10.4921 0.3886
|
|
|
1116 |
---
|
| 61435 |
ripley |
1117 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1118 |
>
|
|
|
1119 |
>
|
|
|
1120 |
>
|
|
|
1121 |
> cleanEx()
|
|
|
1122 |
> nameEx("Puromycin")
|
|
|
1123 |
> ### * Puromycin
|
|
|
1124 |
>
|
|
|
1125 |
> flush(stderr()); flush(stdout())
|
|
|
1126 |
>
|
|
|
1127 |
> ### Name: Puromycin
|
|
|
1128 |
> ### Title: Reaction Velocity of an Enzymatic Reaction
|
|
|
1129 |
> ### Aliases: Puromycin
|
|
|
1130 |
> ### Keywords: datasets
|
|
|
1131 |
>
|
|
|
1132 |
> ### ** Examples
|
|
|
1133 |
>
|
|
|
1134 |
> require(stats); require(graphics)
|
|
|
1135 |
> ## Don't show:
|
|
|
1136 |
> options(show.nls.convergence=FALSE)
|
| 66943 |
maechler |
1137 |
> ## End(Don't show)
|
| 56186 |
murdoch |
1138 |
> plot(rate ~ conc, data = Puromycin, las = 1,
|
|
|
1139 |
+ xlab = "Substrate concentration (ppm)",
|
|
|
1140 |
+ ylab = "Reaction velocity (counts/min/min)",
|
|
|
1141 |
+ pch = as.integer(Puromycin$state),
|
|
|
1142 |
+ col = as.integer(Puromycin$state),
|
|
|
1143 |
+ main = "Puromycin data and fitted Michaelis-Menten curves")
|
|
|
1144 |
> ## simplest form of fitting the Michaelis-Menten model to these data
|
|
|
1145 |
> fm1 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
|
|
|
1146 |
+ subset = state == "treated",
|
|
|
1147 |
+ start = c(Vm = 200, K = 0.05))
|
|
|
1148 |
> fm2 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
|
|
|
1149 |
+ subset = state == "untreated",
|
|
|
1150 |
+ start = c(Vm = 160, K = 0.05))
|
|
|
1151 |
> summary(fm1)
|
|
|
1152 |
|
|
|
1153 |
Formula: rate ~ Vm * conc/(K + conc)
|
|
|
1154 |
|
|
|
1155 |
Parameters:
|
|
|
1156 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1157 |
Vm 2.127e+02 6.947e+00 30.615 3.24e-11 ***
|
|
|
1158 |
K 6.412e-02 8.281e-03 7.743 1.57e-05 ***
|
|
|
1159 |
---
|
| 61435 |
ripley |
1160 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1161 |
|
| 62439 |
ripley |
1162 |
Residual standard error: 10.93 on 10 degrees of freedom
|
| 56186 |
murdoch |
1163 |
|
|
|
1164 |
> summary(fm2)
|
|
|
1165 |
|
|
|
1166 |
Formula: rate ~ Vm * conc/(K + conc)
|
|
|
1167 |
|
|
|
1168 |
Parameters:
|
|
|
1169 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1170 |
Vm 1.603e+02 6.480e+00 24.734 1.38e-09 ***
|
|
|
1171 |
K 4.771e-02 7.782e-03 6.131 0.000173 ***
|
|
|
1172 |
---
|
| 61435 |
ripley |
1173 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1174 |
|
| 62439 |
ripley |
1175 |
Residual standard error: 9.773 on 9 degrees of freedom
|
| 56186 |
murdoch |
1176 |
|
|
|
1177 |
> ## add fitted lines to the plot
|
|
|
1178 |
> conc <- seq(0, 1.2, length.out = 101)
|
|
|
1179 |
> lines(conc, predict(fm1, list(conc = conc)), lty = 1, col = 1)
|
|
|
1180 |
> lines(conc, predict(fm2, list(conc = conc)), lty = 2, col = 2)
|
|
|
1181 |
> legend(0.8, 120, levels(Puromycin$state),
|
|
|
1182 |
+ col = 1:2, lty = 1:2, pch = 1:2)
|
|
|
1183 |
>
|
|
|
1184 |
> ## using partial linearity
|
|
|
1185 |
> fm3 <- nls(rate ~ conc/(K + conc), data = Puromycin,
|
|
|
1186 |
+ subset = state == "treated", start = c(K = 0.05),
|
|
|
1187 |
+ algorithm = "plinear")
|
|
|
1188 |
>
|
|
|
1189 |
>
|
|
|
1190 |
>
|
|
|
1191 |
> cleanEx()
|
|
|
1192 |
> nameEx("Theoph")
|
|
|
1193 |
> ### * Theoph
|
|
|
1194 |
>
|
|
|
1195 |
> flush(stderr()); flush(stdout())
|
|
|
1196 |
>
|
|
|
1197 |
> ### Name: Theoph
|
|
|
1198 |
> ### Title: Pharmacokinetics of Theophylline
|
|
|
1199 |
> ### Aliases: Theoph
|
|
|
1200 |
> ### Keywords: datasets
|
|
|
1201 |
>
|
|
|
1202 |
> ### ** Examples
|
|
|
1203 |
>
|
|
|
1204 |
> require(stats); require(graphics)
|
|
|
1205 |
> ## Don't show:
|
|
|
1206 |
> options(show.nls.convergence=FALSE)
|
| 66943 |
maechler |
1207 |
> ## End(Don't show)
|
| 56186 |
murdoch |
1208 |
> coplot(conc ~ Time | Subject, data = Theoph, show.given = FALSE)
|
|
|
1209 |
> Theoph.4 <- subset(Theoph, Subject == 4)
|
|
|
1210 |
> fm1 <- nls(conc ~ SSfol(Dose, Time, lKe, lKa, lCl),
|
|
|
1211 |
+ data = Theoph.4)
|
|
|
1212 |
> summary(fm1)
|
|
|
1213 |
|
|
|
1214 |
Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl)
|
|
|
1215 |
|
|
|
1216 |
Parameters:
|
|
|
1217 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1218 |
lKe -2.4365 0.2257 -10.797 4.77e-06 ***
|
|
|
1219 |
lKa 0.1583 0.2297 0.689 0.51
|
|
|
1220 |
lCl -3.2861 0.1448 -22.695 1.51e-08 ***
|
|
|
1221 |
---
|
| 61435 |
ripley |
1222 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1223 |
|
| 62439 |
ripley |
1224 |
Residual standard error: 0.8465 on 8 degrees of freedom
|
| 56186 |
murdoch |
1225 |
|
|
|
1226 |
> plot(conc ~ Time, data = Theoph.4,
|
|
|
1227 |
+ xlab = "Time since drug administration (hr)",
|
|
|
1228 |
+ ylab = "Theophylline concentration (mg/L)",
|
|
|
1229 |
+ main = "Observed concentrations and fitted model",
|
|
|
1230 |
+ sub = "Theophylline data - Subject 4 only",
|
|
|
1231 |
+ las = 1, col = 4)
|
|
|
1232 |
> xvals <- seq(0, par("usr")[2], length.out = 55)
|
|
|
1233 |
> lines(xvals, predict(fm1, newdata = list(Time = xvals)),
|
|
|
1234 |
+ col = 4)
|
|
|
1235 |
>
|
|
|
1236 |
>
|
|
|
1237 |
>
|
|
|
1238 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1239 |
> cleanEx()
|
|
|
1240 |
> nameEx("Titanic")
|
|
|
1241 |
> ### * Titanic
|
|
|
1242 |
>
|
|
|
1243 |
> flush(stderr()); flush(stdout())
|
|
|
1244 |
>
|
|
|
1245 |
> ### Name: Titanic
|
|
|
1246 |
> ### Title: Survival of passengers on the Titanic
|
|
|
1247 |
> ### Aliases: Titanic
|
|
|
1248 |
> ### Keywords: datasets
|
|
|
1249 |
>
|
|
|
1250 |
> ### ** Examples
|
|
|
1251 |
>
|
|
|
1252 |
> require(graphics)
|
|
|
1253 |
> mosaicplot(Titanic, main = "Survival on the Titanic")
|
|
|
1254 |
> ## Higher survival rates in children?
|
|
|
1255 |
> apply(Titanic, c(3, 4), sum)
|
|
|
1256 |
Survived
|
|
|
1257 |
Age No Yes
|
|
|
1258 |
Child 52 57
|
|
|
1259 |
Adult 1438 654
|
|
|
1260 |
> ## Higher survival rates in females?
|
|
|
1261 |
> apply(Titanic, c(2, 4), sum)
|
|
|
1262 |
Survived
|
|
|
1263 |
Sex No Yes
|
|
|
1264 |
Male 1364 367
|
|
|
1265 |
Female 126 344
|
|
|
1266 |
> ## Use loglm() in package 'MASS' for further analysis ...
|
|
|
1267 |
>
|
|
|
1268 |
>
|
|
|
1269 |
>
|
|
|
1270 |
> cleanEx()
|
|
|
1271 |
> nameEx("ToothGrowth")
|
|
|
1272 |
> ### * ToothGrowth
|
|
|
1273 |
>
|
|
|
1274 |
> flush(stderr()); flush(stdout())
|
|
|
1275 |
>
|
|
|
1276 |
> ### Name: ToothGrowth
|
|
|
1277 |
> ### Title: The Effect of Vitamin C on Tooth Growth in Guinea Pigs
|
|
|
1278 |
> ### Aliases: ToothGrowth
|
|
|
1279 |
> ### Keywords: datasets
|
|
|
1280 |
>
|
|
|
1281 |
> ### ** Examples
|
|
|
1282 |
>
|
|
|
1283 |
> require(graphics)
|
|
|
1284 |
> coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
|
|
|
1285 |
+ xlab = "ToothGrowth data: length vs dose, given type of supplement")
|
|
|
1286 |
>
|
|
|
1287 |
>
|
|
|
1288 |
>
|
|
|
1289 |
> cleanEx()
|
|
|
1290 |
> nameEx("UCBAdmissions")
|
|
|
1291 |
> ### * UCBAdmissions
|
|
|
1292 |
>
|
|
|
1293 |
> flush(stderr()); flush(stdout())
|
|
|
1294 |
>
|
|
|
1295 |
> ### Name: UCBAdmissions
|
|
|
1296 |
> ### Title: Student Admissions at UC Berkeley
|
|
|
1297 |
> ### Aliases: UCBAdmissions
|
|
|
1298 |
> ### Keywords: datasets
|
|
|
1299 |
>
|
|
|
1300 |
> ### ** Examples
|
|
|
1301 |
>
|
|
|
1302 |
> require(graphics)
|
|
|
1303 |
> ## Data aggregated over departments
|
|
|
1304 |
> apply(UCBAdmissions, c(1, 2), sum)
|
|
|
1305 |
Gender
|
|
|
1306 |
Admit Male Female
|
|
|
1307 |
Admitted 1198 557
|
|
|
1308 |
Rejected 1493 1278
|
|
|
1309 |
> mosaicplot(apply(UCBAdmissions, c(1, 2), sum),
|
|
|
1310 |
+ main = "Student admissions at UC Berkeley")
|
|
|
1311 |
> ## Data for individual departments
|
|
|
1312 |
> opar <- par(mfrow = c(2, 3), oma = c(0, 0, 2, 0))
|
|
|
1313 |
> for(i in 1:6)
|
|
|
1314 |
+ mosaicplot(UCBAdmissions[,,i],
|
|
|
1315 |
+ xlab = "Admit", ylab = "Sex",
|
|
|
1316 |
+ main = paste("Department", LETTERS[i]))
|
|
|
1317 |
> mtext(expression(bold("Student admissions at UC Berkeley")),
|
|
|
1318 |
+ outer = TRUE, cex = 1.5)
|
|
|
1319 |
> par(opar)
|
|
|
1320 |
>
|
|
|
1321 |
>
|
|
|
1322 |
>
|
|
|
1323 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1324 |
> cleanEx()
|
|
|
1325 |
> nameEx("UKDriverDeaths")
|
|
|
1326 |
> ### * UKDriverDeaths
|
|
|
1327 |
>
|
|
|
1328 |
> flush(stderr()); flush(stdout())
|
|
|
1329 |
>
|
|
|
1330 |
> ### Name: UKDriverDeaths
|
|
|
1331 |
> ### Title: Road Casualties in Great Britain 1969-84
|
|
|
1332 |
> ### Aliases: UKDriverDeaths Seatbelts
|
|
|
1333 |
> ### Keywords: datasets
|
|
|
1334 |
>
|
|
|
1335 |
> ### ** Examples
|
|
|
1336 |
>
|
|
|
1337 |
> require(stats); require(graphics)
|
|
|
1338 |
> ## work with pre-seatbelt period to identify a model, use logs
|
|
|
1339 |
> work <- window(log10(UKDriverDeaths), end = 1982+11/12)
|
| 61169 |
ripley |
1340 |
> par(mfrow = c(3, 1))
|
| 56186 |
murdoch |
1341 |
> plot(work); acf(work); pacf(work)
|
| 61169 |
ripley |
1342 |
> par(mfrow = c(1, 1))
|
|
|
1343 |
> (fit <- arima(work, c(1, 0, 0), seasonal = list(order = c(1, 0, 0))))
|
| 56186 |
murdoch |
1344 |
|
|
|
1345 |
Call:
|
|
|
1346 |
arima(x = work, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0)))
|
|
|
1347 |
|
|
|
1348 |
Coefficients:
|
|
|
1349 |
ar1 sar1 intercept
|
|
|
1350 |
0.4378 0.6281 3.2274
|
|
|
1351 |
s.e. 0.0764 0.0637 0.0131
|
|
|
1352 |
|
|
|
1353 |
sigma^2 estimated as 0.00157: log likelihood = 300.85, aic = -593.7
|
|
|
1354 |
> z <- predict(fit, n.ahead = 24)
|
|
|
1355 |
> ts.plot(log10(UKDriverDeaths), z$pred, z$pred+2*z$se, z$pred-2*z$se,
|
| 61169 |
ripley |
1356 |
+ lty = c(1, 3, 2, 2), col = c("black", "red", "blue", "blue"))
|
| 56186 |
murdoch |
1357 |
>
|
|
|
1358 |
> ## now see the effect of the explanatory variables
|
|
|
1359 |
> X <- Seatbelts[, c("kms", "PetrolPrice", "law")]
|
|
|
1360 |
> X[, 1] <- log10(X[, 1]) - 4
|
| 61169 |
ripley |
1361 |
> arima(log10(Seatbelts[, "drivers"]), c(1, 0, 0),
|
|
|
1362 |
+ seasonal = list(order = c(1, 0, 0)), xreg = X)
|
| 56186 |
murdoch |
1363 |
|
|
|
1364 |
Call:
|
|
|
1365 |
arima(x = log10(Seatbelts[, "drivers"]), order = c(1, 0, 0), seasonal = list(order = c(1,
|
|
|
1366 |
0, 0)), xreg = X)
|
|
|
1367 |
|
|
|
1368 |
Coefficients:
|
|
|
1369 |
ar1 sar1 intercept kms PetrolPrice law
|
|
|
1370 |
0.3348 0.6672 3.3539 0.0082 -1.2224 -0.0963
|
|
|
1371 |
s.e. 0.0775 0.0612 0.0441 0.0902 0.3839 0.0166
|
|
|
1372 |
|
|
|
1373 |
sigma^2 estimated as 0.001476: log likelihood = 349.73, aic = -685.46
|
|
|
1374 |
>
|
|
|
1375 |
>
|
|
|
1376 |
>
|
|
|
1377 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1378 |
> cleanEx()
|
|
|
1379 |
> nameEx("UKLungDeaths")
|
|
|
1380 |
> ### * UKLungDeaths
|
|
|
1381 |
>
|
|
|
1382 |
> flush(stderr()); flush(stdout())
|
|
|
1383 |
>
|
|
|
1384 |
> ### Name: UKLungDeaths
|
|
|
1385 |
> ### Title: Monthly Deaths from Lung Diseases in the UK
|
|
|
1386 |
> ### Aliases: UKLungDeaths ldeaths fdeaths mdeaths
|
|
|
1387 |
> ### Keywords: datasets
|
|
|
1388 |
>
|
|
|
1389 |
> ### ** Examples
|
|
|
1390 |
>
|
|
|
1391 |
> require(stats); require(graphics) # for time
|
| 61435 |
ripley |
1392 |
> plot(ldeaths)
|
|
|
1393 |
> plot(mdeaths, fdeaths)
|
| 56186 |
murdoch |
1394 |
> ## Better labels:
|
|
|
1395 |
> yr <- floor(tt <- time(mdeaths))
|
|
|
1396 |
> plot(mdeaths, fdeaths,
|
| 61157 |
ripley |
1397 |
+ xy.labels = paste(month.abb[12*(tt - yr)], yr-1900, sep = "'"))
|
| 56186 |
murdoch |
1398 |
>
|
|
|
1399 |
>
|
|
|
1400 |
>
|
|
|
1401 |
> cleanEx()
|
|
|
1402 |
> nameEx("UKgas")
|
|
|
1403 |
> ### * UKgas
|
|
|
1404 |
>
|
|
|
1405 |
> flush(stderr()); flush(stdout())
|
|
|
1406 |
>
|
|
|
1407 |
> ### Name: UKgas
|
|
|
1408 |
> ### Title: UK Quarterly Gas Consumption
|
|
|
1409 |
> ### Aliases: UKgas
|
|
|
1410 |
> ### Keywords: datasets
|
|
|
1411 |
>
|
|
|
1412 |
> ### ** Examples
|
|
|
1413 |
>
|
|
|
1414 |
> ## maybe str(UKgas) ; plot(UKgas) ...
|
|
|
1415 |
>
|
|
|
1416 |
>
|
|
|
1417 |
>
|
|
|
1418 |
> cleanEx()
|
|
|
1419 |
> nameEx("USArrests")
|
|
|
1420 |
> ### * USArrests
|
|
|
1421 |
>
|
|
|
1422 |
> flush(stderr()); flush(stdout())
|
|
|
1423 |
>
|
|
|
1424 |
> ### Name: USArrests
|
|
|
1425 |
> ### Title: Violent Crime Rates by US State
|
|
|
1426 |
> ### Aliases: USArrests
|
|
|
1427 |
> ### Keywords: datasets
|
|
|
1428 |
>
|
|
|
1429 |
> ### ** Examples
|
|
|
1430 |
>
|
| 71632 |
ripley |
1431 |
> summary(USArrests)
|
|
|
1432 |
Murder Assault UrbanPop Rape
|
|
|
1433 |
Min. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.30
|
|
|
1434 |
1st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07
|
|
|
1435 |
Median : 7.250 Median :159.0 Median :66.00 Median :20.10
|
|
|
1436 |
Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.23
|
|
|
1437 |
3rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18
|
|
|
1438 |
Max. :17.400 Max. :337.0 Max. :91.00 Max. :46.00
|
|
|
1439 |
>
|
| 56186 |
murdoch |
1440 |
> require(graphics)
|
|
|
1441 |
> pairs(USArrests, panel = panel.smooth, main = "USArrests data")
|
|
|
1442 |
>
|
| 71632 |
ripley |
1443 |
> ## Difference between 'USArrests' and its correction
|
|
|
1444 |
> USArrests["Maryland", "UrbanPop"] # 67 -- the transcription error
|
|
|
1445 |
[1] 67
|
|
|
1446 |
> UA.C <- USArrests
|
|
|
1447 |
> UA.C["Maryland", "UrbanPop"] <- 76.6
|
| 56186 |
murdoch |
1448 |
>
|
| 71632 |
ripley |
1449 |
> ## also +/- 0.5 to restore the original <n>.5 percentages
|
|
|
1450 |
> s5u <- c("Colorado", "Florida", "Mississippi", "Wyoming")
|
|
|
1451 |
> s5d <- c("Nebraska", "Pennsylvania")
|
|
|
1452 |
> UA.C[s5u, "UrbanPop"] <- UA.C[s5u, "UrbanPop"] + 0.5
|
|
|
1453 |
> UA.C[s5d, "UrbanPop"] <- UA.C[s5d, "UrbanPop"] - 0.5
|
| 56186 |
murdoch |
1454 |
>
|
| 71632 |
ripley |
1455 |
> ## ==> UA.C is now a *C*orrected version of USArrests
|
|
|
1456 |
>
|
|
|
1457 |
>
|
|
|
1458 |
>
|
| 56186 |
murdoch |
1459 |
> cleanEx()
|
|
|
1460 |
> nameEx("USJudgeRatings")
|
|
|
1461 |
> ### * USJudgeRatings
|
|
|
1462 |
>
|
|
|
1463 |
> flush(stderr()); flush(stdout())
|
|
|
1464 |
>
|
|
|
1465 |
> ### Name: USJudgeRatings
|
|
|
1466 |
> ### Title: Lawyers' Ratings of State Judges in the US Superior Court
|
|
|
1467 |
> ### Aliases: USJudgeRatings
|
|
|
1468 |
> ### Keywords: datasets
|
|
|
1469 |
>
|
|
|
1470 |
> ### ** Examples
|
|
|
1471 |
>
|
|
|
1472 |
> require(graphics)
|
|
|
1473 |
> pairs(USJudgeRatings, main = "USJudgeRatings data")
|
|
|
1474 |
>
|
|
|
1475 |
>
|
|
|
1476 |
>
|
|
|
1477 |
> cleanEx()
|
|
|
1478 |
> nameEx("USPersonalExpenditure")
|
|
|
1479 |
> ### * USPersonalExpenditure
|
|
|
1480 |
>
|
|
|
1481 |
> flush(stderr()); flush(stdout())
|
|
|
1482 |
>
|
|
|
1483 |
> ### Name: USPersonalExpenditure
|
|
|
1484 |
> ### Title: Personal Expenditure Data
|
|
|
1485 |
> ### Aliases: USPersonalExpenditure
|
|
|
1486 |
> ### Keywords: datasets
|
|
|
1487 |
>
|
|
|
1488 |
> ### ** Examples
|
|
|
1489 |
>
|
|
|
1490 |
> require(stats) # for medpolish
|
|
|
1491 |
> USPersonalExpenditure
|
|
|
1492 |
1940 1945 1950 1955 1960
|
|
|
1493 |
Food and Tobacco 22.200 44.500 59.60 73.2 86.80
|
|
|
1494 |
Household Operation 10.500 15.500 29.00 36.5 46.20
|
|
|
1495 |
Medical and Health 3.530 5.760 9.71 14.0 21.10
|
|
|
1496 |
Personal Care 1.040 1.980 2.45 3.4 5.40
|
|
|
1497 |
Private Education 0.341 0.974 1.80 2.6 3.64
|
|
|
1498 |
> medpolish(log10(USPersonalExpenditure))
|
| 61435 |
ripley |
1499 |
1: 1.126317
|
|
|
1500 |
2: 1.032421
|
|
|
1501 |
Final: 1.032421
|
| 56186 |
murdoch |
1502 |
|
|
|
1503 |
Median Polish Results (Dataset: "log10(USPersonalExpenditure)")
|
|
|
1504 |
|
| 61435 |
ripley |
1505 |
Overall: 0.9872192
|
| 56186 |
murdoch |
1506 |
|
|
|
1507 |
Row Effects:
|
|
|
1508 |
Food and Tobacco Household Operation Medical and Health Personal Care
|
|
|
1509 |
0.7880270 0.4327608 0.0000000 -0.5606543
|
|
|
1510 |
Private Education
|
|
|
1511 |
-0.7319467
|
|
|
1512 |
|
|
|
1513 |
Column Effects:
|
|
|
1514 |
1940 1945 1950 1955 1960
|
|
|
1515 |
-0.4288933 -0.2267967 0.0000000 0.1423128 0.3058289
|
|
|
1516 |
|
|
|
1517 |
Residuals:
|
|
|
1518 |
1940 1945 1950 1955 1960
|
|
|
1519 |
Food and Tobacco 0.000000 0.0999105 0.000000 -0.053048 -0.142555
|
|
|
1520 |
Household Operation 0.030103 -0.0028516 0.042418 0.000000 -0.061167
|
|
|
1521 |
Medical and Health -0.010551 0.0000000 0.000000 0.016596 0.031234
|
|
|
1522 |
Personal Care 0.019362 0.0968971 -0.037399 -0.037399 0.000000
|
|
|
1523 |
Private Education -0.293625 -0.0399168 0.000000 0.017388 0.000000
|
|
|
1524 |
|
|
|
1525 |
>
|
|
|
1526 |
>
|
|
|
1527 |
>
|
|
|
1528 |
> cleanEx()
|
|
|
1529 |
> nameEx("VADeaths")
|
|
|
1530 |
> ### * VADeaths
|
|
|
1531 |
>
|
|
|
1532 |
> flush(stderr()); flush(stdout())
|
|
|
1533 |
>
|
|
|
1534 |
> ### Name: VADeaths
|
|
|
1535 |
> ### Title: Death Rates in Virginia (1940)
|
|
|
1536 |
> ### Aliases: VADeaths
|
|
|
1537 |
> ### Keywords: datasets
|
|
|
1538 |
>
|
|
|
1539 |
> ### ** Examples
|
|
|
1540 |
>
|
|
|
1541 |
> require(stats); require(graphics)
|
|
|
1542 |
> n <- length(dr <- c(VADeaths))
|
|
|
1543 |
> nam <- names(VADeaths)
|
|
|
1544 |
> d.VAD <- data.frame(
|
|
|
1545 |
+ Drate = dr,
|
| 61157 |
ripley |
1546 |
+ age = rep(ordered(rownames(VADeaths)), length.out = n),
|
| 61169 |
ripley |
1547 |
+ gender = gl(2, 5, n, labels = c("M", "F")),
|
|
|
1548 |
+ site = gl(2, 10, labels = c("rural", "urban")))
|
| 56186 |
murdoch |
1549 |
> coplot(Drate ~ as.numeric(age) | gender * site, data = d.VAD,
|
|
|
1550 |
+ panel = panel.smooth, xlab = "VADeaths data - Given: gender")
|
|
|
1551 |
> summary(aov.VAD <- aov(Drate ~ .^2, data = d.VAD))
|
| 57044 |
ripley |
1552 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1553 |
age 4 6288 1572.1 590.858 8.55e-06 ***
|
|
|
1554 |
gender 1 648 647.5 243.361 9.86e-05 ***
|
|
|
1555 |
site 1 77 76.8 28.876 0.00579 **
|
|
|
1556 |
age:gender 4 86 21.6 8.100 0.03358 *
|
|
|
1557 |
age:site 4 43 10.6 3.996 0.10414
|
|
|
1558 |
gender:site 1 73 73.0 27.422 0.00636 **
|
|
|
1559 |
Residuals 4 11 2.7
|
| 56186 |
murdoch |
1560 |
---
|
| 61435 |
ripley |
1561 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 61169 |
ripley |
1562 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
1563 |
> plot(aov.VAD)
|
|
|
1564 |
> par(opar)
|
|
|
1565 |
>
|
|
|
1566 |
>
|
|
|
1567 |
>
|
|
|
1568 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1569 |
> cleanEx()
|
|
|
1570 |
> nameEx("WWWusage")
|
|
|
1571 |
> ### * WWWusage
|
|
|
1572 |
>
|
|
|
1573 |
> flush(stderr()); flush(stdout())
|
|
|
1574 |
>
|
|
|
1575 |
> ### Name: WWWusage
|
|
|
1576 |
> ### Title: Internet Usage per Minute
|
|
|
1577 |
> ### Aliases: WWWusage
|
|
|
1578 |
> ### Keywords: datasets
|
|
|
1579 |
>
|
|
|
1580 |
> ### ** Examples
|
|
|
1581 |
>
|
|
|
1582 |
> require(graphics)
|
|
|
1583 |
> work <- diff(WWWusage)
|
| 61169 |
ripley |
1584 |
> par(mfrow = c(2, 1)); plot(WWWusage); plot(work)
|
| 56186 |
murdoch |
1585 |
> ## Not run:
|
|
|
1586 |
> ##D require(stats)
|
| 61157 |
ripley |
1587 |
> ##D aics <- matrix(, 6, 6, dimnames = list(p = 0:5, q = 0:5))
|
| 61169 |
ripley |
1588 |
> ##D for(q in 1:5) aics[1, 1+q] <- arima(WWWusage, c(0, 1, q),
|
| 56186 |
murdoch |
1589 |
> ##D optim.control = list(maxit = 500))$aic
|
|
|
1590 |
> ##D for(p in 1:5)
|
| 61169 |
ripley |
1591 |
> ##D for(q in 0:5) aics[1+p, 1+q] <- arima(WWWusage, c(p, 1, q),
|
| 56186 |
murdoch |
1592 |
> ##D optim.control = list(maxit = 500))$aic
|
| 61157 |
ripley |
1593 |
> ##D round(aics - min(aics, na.rm = TRUE), 2)
|
| 56186 |
murdoch |
1594 |
> ## End(Not run)
|
|
|
1595 |
>
|
|
|
1596 |
>
|
|
|
1597 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1598 |
> cleanEx()
|
|
|
1599 |
> nameEx("WorldPhones")
|
|
|
1600 |
> ### * WorldPhones
|
|
|
1601 |
>
|
|
|
1602 |
> flush(stderr()); flush(stdout())
|
|
|
1603 |
>
|
|
|
1604 |
> ### Name: WorldPhones
|
|
|
1605 |
> ### Title: The World's Telephones
|
|
|
1606 |
> ### Aliases: WorldPhones
|
|
|
1607 |
> ### Keywords: datasets
|
|
|
1608 |
>
|
|
|
1609 |
> ### ** Examples
|
|
|
1610 |
>
|
|
|
1611 |
> require(graphics)
|
|
|
1612 |
> matplot(rownames(WorldPhones), WorldPhones, type = "b", log = "y",
|
|
|
1613 |
+ xlab = "Year", ylab = "Number of telephones (1000's)")
|
| 61435 |
ripley |
1614 |
> legend(1951.5, 80000, colnames(WorldPhones), col = 1:6, lty = 1:5,
|
| 56186 |
murdoch |
1615 |
+ pch = rep(21, 7))
|
|
|
1616 |
> title(main = "World phones data: log scale for response")
|
|
|
1617 |
>
|
|
|
1618 |
>
|
|
|
1619 |
>
|
|
|
1620 |
> cleanEx()
|
|
|
1621 |
> nameEx("ability.cov")
|
|
|
1622 |
> ### * ability.cov
|
|
|
1623 |
>
|
|
|
1624 |
> flush(stderr()); flush(stdout())
|
|
|
1625 |
>
|
|
|
1626 |
> ### Name: ability.cov
|
|
|
1627 |
> ### Title: Ability and Intelligence Tests
|
|
|
1628 |
> ### Aliases: ability.cov
|
|
|
1629 |
> ### Keywords: datasets
|
|
|
1630 |
>
|
|
|
1631 |
> ### ** Examples
|
|
|
1632 |
>
|
|
|
1633 |
>
|
|
|
1634 |
> cleanEx()
|
|
|
1635 |
> nameEx("airmiles")
|
|
|
1636 |
> ### * airmiles
|
|
|
1637 |
>
|
|
|
1638 |
> flush(stderr()); flush(stdout())
|
|
|
1639 |
>
|
|
|
1640 |
> ### Name: airmiles
|
|
|
1641 |
> ### Title: Passenger Miles on Commercial US Airlines, 1937-1960
|
|
|
1642 |
> ### Aliases: airmiles
|
|
|
1643 |
> ### Keywords: datasets
|
|
|
1644 |
>
|
|
|
1645 |
> ### ** Examples
|
|
|
1646 |
>
|
|
|
1647 |
> require(graphics)
|
|
|
1648 |
> plot(airmiles, main = "airmiles data",
|
|
|
1649 |
+ xlab = "Passenger-miles flown by U.S. commercial airlines", col = 4)
|
|
|
1650 |
>
|
|
|
1651 |
>
|
|
|
1652 |
>
|
|
|
1653 |
> cleanEx()
|
|
|
1654 |
> nameEx("airquality")
|
|
|
1655 |
> ### * airquality
|
|
|
1656 |
>
|
|
|
1657 |
> flush(stderr()); flush(stdout())
|
|
|
1658 |
>
|
|
|
1659 |
> ### Name: airquality
|
|
|
1660 |
> ### Title: New York Air Quality Measurements
|
|
|
1661 |
> ### Aliases: airquality
|
|
|
1662 |
> ### Keywords: datasets
|
|
|
1663 |
>
|
|
|
1664 |
> ### ** Examples
|
|
|
1665 |
>
|
|
|
1666 |
> require(graphics)
|
|
|
1667 |
> pairs(airquality, panel = panel.smooth, main = "airquality data")
|
|
|
1668 |
>
|
|
|
1669 |
>
|
|
|
1670 |
>
|
|
|
1671 |
> cleanEx()
|
|
|
1672 |
> nameEx("anscombe")
|
|
|
1673 |
> ### * anscombe
|
|
|
1674 |
>
|
|
|
1675 |
> flush(stderr()); flush(stdout())
|
|
|
1676 |
>
|
|
|
1677 |
> ### Name: anscombe
|
|
|
1678 |
> ### Title: Anscombe's Quartet of 'Identical' Simple Linear Regressions
|
|
|
1679 |
> ### Aliases: anscombe
|
|
|
1680 |
> ### Keywords: datasets
|
|
|
1681 |
>
|
|
|
1682 |
> ### ** Examples
|
|
|
1683 |
>
|
|
|
1684 |
> require(stats); require(graphics)
|
|
|
1685 |
> summary(anscombe)
|
|
|
1686 |
x1 x2 x3 x4 y1
|
|
|
1687 |
Min. : 4.0 Min. : 4.0 Min. : 4.0 Min. : 8 Min. : 4.260
|
|
|
1688 |
1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 8 1st Qu.: 6.315
|
|
|
1689 |
Median : 9.0 Median : 9.0 Median : 9.0 Median : 8 Median : 7.580
|
|
|
1690 |
Mean : 9.0 Mean : 9.0 Mean : 9.0 Mean : 9 Mean : 7.501
|
|
|
1691 |
3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.: 8 3rd Qu.: 8.570
|
|
|
1692 |
Max. :14.0 Max. :14.0 Max. :14.0 Max. :19 Max. :10.840
|
|
|
1693 |
y2 y3 y4
|
|
|
1694 |
Min. :3.100 Min. : 5.39 Min. : 5.250
|
|
|
1695 |
1st Qu.:6.695 1st Qu.: 6.25 1st Qu.: 6.170
|
|
|
1696 |
Median :8.140 Median : 7.11 Median : 7.040
|
|
|
1697 |
Mean :7.501 Mean : 7.50 Mean : 7.501
|
|
|
1698 |
3rd Qu.:8.950 3rd Qu.: 7.98 3rd Qu.: 8.190
|
|
|
1699 |
Max. :9.260 Max. :12.74 Max. :12.500
|
|
|
1700 |
>
|
|
|
1701 |
> ##-- now some "magic" to do the 4 regressions in a loop:
|
|
|
1702 |
> ff <- y ~ x
|
| 57821 |
maechler |
1703 |
> mods <- setNames(as.list(1:4), paste0("lm", 1:4))
|
| 56186 |
murdoch |
1704 |
> for(i in 1:4) {
|
| 57821 |
maechler |
1705 |
+ ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
|
|
|
1706 |
+ ## or ff[[2]] <- as.name(paste0("y", i))
|
|
|
1707 |
+ ## ff[[3]] <- as.name(paste0("x", i))
|
| 61157 |
ripley |
1708 |
+ mods[[i]] <- lmi <- lm(ff, data = anscombe)
|
| 56186 |
murdoch |
1709 |
+ print(anova(lmi))
|
|
|
1710 |
+ }
|
|
|
1711 |
Analysis of Variance Table
|
|
|
1712 |
|
|
|
1713 |
Response: y1
|
|
|
1714 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1715 |
x1 1 27.510 27.5100 17.99 0.00217 **
|
|
|
1716 |
Residuals 9 13.763 1.5292
|
|
|
1717 |
---
|
| 61435 |
ripley |
1718 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1719 |
Analysis of Variance Table
|
|
|
1720 |
|
|
|
1721 |
Response: y2
|
|
|
1722 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1723 |
x2 1 27.500 27.5000 17.966 0.002179 **
|
|
|
1724 |
Residuals 9 13.776 1.5307
|
|
|
1725 |
---
|
| 61435 |
ripley |
1726 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1727 |
Analysis of Variance Table
|
|
|
1728 |
|
|
|
1729 |
Response: y3
|
|
|
1730 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1731 |
x3 1 27.470 27.4700 17.972 0.002176 **
|
|
|
1732 |
Residuals 9 13.756 1.5285
|
|
|
1733 |
---
|
| 61435 |
ripley |
1734 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1735 |
Analysis of Variance Table
|
|
|
1736 |
|
|
|
1737 |
Response: y4
|
|
|
1738 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1739 |
x4 1 27.490 27.4900 18.003 0.002165 **
|
|
|
1740 |
Residuals 9 13.742 1.5269
|
|
|
1741 |
---
|
| 61435 |
ripley |
1742 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1743 |
>
|
|
|
1744 |
> ## See how close they are (numerically!)
|
| 57821 |
maechler |
1745 |
> sapply(mods, coef)
|
|
|
1746 |
lm1 lm2 lm3 lm4
|
| 56186 |
murdoch |
1747 |
(Intercept) 3.0000909 3.000909 3.0024545 3.0017273
|
|
|
1748 |
x1 0.5000909 0.500000 0.4997273 0.4999091
|
| 57821 |
maechler |
1749 |
> lapply(mods, function(fm) coef(summary(fm)))
|
|
|
1750 |
$lm1
|
| 56186 |
murdoch |
1751 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1752 |
(Intercept) 3.0000909 1.1247468 2.667348 0.025734051
|
|
|
1753 |
x1 0.5000909 0.1179055 4.241455 0.002169629
|
|
|
1754 |
|
| 57821 |
maechler |
1755 |
$lm2
|
| 56186 |
murdoch |
1756 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1757 |
(Intercept) 3.000909 1.1253024 2.666758 0.025758941
|
|
|
1758 |
x2 0.500000 0.1179637 4.238590 0.002178816
|
|
|
1759 |
|
| 57821 |
maechler |
1760 |
$lm3
|
| 56186 |
murdoch |
1761 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1762 |
(Intercept) 3.0024545 1.1244812 2.670080 0.025619109
|
|
|
1763 |
x3 0.4997273 0.1178777 4.239372 0.002176305
|
|
|
1764 |
|
| 57821 |
maechler |
1765 |
$lm4
|
| 56186 |
murdoch |
1766 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1767 |
(Intercept) 3.0017273 1.1239211 2.670763 0.025590425
|
|
|
1768 |
x4 0.4999091 0.1178189 4.243028 0.002164602
|
|
|
1769 |
|
|
|
1770 |
>
|
|
|
1771 |
> ## Now, do what you should have done in the first place: PLOTS
|
| 61169 |
ripley |
1772 |
> op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma = c(0, 0, 2, 0))
|
| 56186 |
murdoch |
1773 |
> for(i in 1:4) {
|
| 57821 |
maechler |
1774 |
+ ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
|
| 61157 |
ripley |
1775 |
+ plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
|
| 61169 |
ripley |
1776 |
+ xlim = c(3, 19), ylim = c(3, 13))
|
| 61157 |
ripley |
1777 |
+ abline(mods[[i]], col = "blue")
|
| 56186 |
murdoch |
1778 |
+ }
|
| 61157 |
ripley |
1779 |
> mtext("Anscombe's 4 Regression data sets", outer = TRUE, cex = 1.5)
|
| 56186 |
murdoch |
1780 |
> par(op)
|
|
|
1781 |
>
|
|
|
1782 |
>
|
|
|
1783 |
>
|
|
|
1784 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1785 |
> cleanEx()
|
|
|
1786 |
> nameEx("attenu")
|
|
|
1787 |
> ### * attenu
|
|
|
1788 |
>
|
|
|
1789 |
> flush(stderr()); flush(stdout())
|
|
|
1790 |
>
|
|
|
1791 |
> ### Name: attenu
|
|
|
1792 |
> ### Title: The Joyner-Boore Attenuation Data
|
|
|
1793 |
> ### Aliases: attenu
|
|
|
1794 |
> ### Keywords: datasets
|
|
|
1795 |
>
|
|
|
1796 |
> ### ** Examples
|
|
|
1797 |
>
|
|
|
1798 |
> require(graphics)
|
|
|
1799 |
> ## check the data class of the variables
|
|
|
1800 |
> sapply(attenu, data.class)
|
|
|
1801 |
event mag station dist accel
|
|
|
1802 |
"numeric" "numeric" "factor" "numeric" "numeric"
|
|
|
1803 |
> summary(attenu)
|
|
|
1804 |
event mag station dist
|
|
|
1805 |
Min. : 1.00 Min. :5.000 117 : 5 Min. : 0.50
|
|
|
1806 |
1st Qu.: 9.00 1st Qu.:5.300 1028 : 4 1st Qu.: 11.32
|
|
|
1807 |
Median :18.00 Median :6.100 113 : 4 Median : 23.40
|
|
|
1808 |
Mean :14.74 Mean :6.084 112 : 3 Mean : 45.60
|
|
|
1809 |
3rd Qu.:20.00 3rd Qu.:6.600 135 : 3 3rd Qu.: 47.55
|
|
|
1810 |
Max. :23.00 Max. :7.700 (Other):147 Max. :370.00
|
| 89306 |
smeyer |
1811 |
NAs : 16
|
| 56186 |
murdoch |
1812 |
accel
|
|
|
1813 |
Min. :0.00300
|
|
|
1814 |
1st Qu.:0.04425
|
|
|
1815 |
Median :0.11300
|
|
|
1816 |
Mean :0.15422
|
|
|
1817 |
3rd Qu.:0.21925
|
|
|
1818 |
Max. :0.81000
|
|
|
1819 |
|
|
|
1820 |
> pairs(attenu, main = "attenu data")
|
|
|
1821 |
> coplot(accel ~ dist | as.factor(event), data = attenu, show.given = FALSE)
|
|
|
1822 |
> coplot(log(accel) ~ log(dist) | as.factor(event),
|
|
|
1823 |
+ data = attenu, panel = panel.smooth, show.given = FALSE)
|
|
|
1824 |
>
|
|
|
1825 |
>
|
|
|
1826 |
>
|
|
|
1827 |
> cleanEx()
|
|
|
1828 |
> nameEx("attitude")
|
|
|
1829 |
> ### * attitude
|
|
|
1830 |
>
|
|
|
1831 |
> flush(stderr()); flush(stdout())
|
|
|
1832 |
>
|
|
|
1833 |
> ### Name: attitude
|
|
|
1834 |
> ### Title: The Chatterjee-Price Attitude Data
|
|
|
1835 |
> ### Aliases: attitude
|
|
|
1836 |
> ### Keywords: datasets
|
|
|
1837 |
>
|
|
|
1838 |
> ### ** Examples
|
|
|
1839 |
>
|
|
|
1840 |
> require(stats); require(graphics)
|
|
|
1841 |
> pairs(attitude, main = "attitude data")
|
|
|
1842 |
> summary(attitude)
|
|
|
1843 |
rating complaints privileges learning raises
|
|
|
1844 |
Min. :40.00 Min. :37.0 Min. :30.00 Min. :34.00 Min. :43.00
|
|
|
1845 |
1st Qu.:58.75 1st Qu.:58.5 1st Qu.:45.00 1st Qu.:47.00 1st Qu.:58.25
|
|
|
1846 |
Median :65.50 Median :65.0 Median :51.50 Median :56.50 Median :63.50
|
|
|
1847 |
Mean :64.63 Mean :66.6 Mean :53.13 Mean :56.37 Mean :64.63
|
|
|
1848 |
3rd Qu.:71.75 3rd Qu.:77.0 3rd Qu.:62.50 3rd Qu.:66.75 3rd Qu.:71.00
|
|
|
1849 |
Max. :85.00 Max. :90.0 Max. :83.00 Max. :75.00 Max. :88.00
|
|
|
1850 |
critical advance
|
|
|
1851 |
Min. :49.00 Min. :25.00
|
|
|
1852 |
1st Qu.:69.25 1st Qu.:35.00
|
|
|
1853 |
Median :77.50 Median :41.00
|
|
|
1854 |
Mean :74.77 Mean :42.93
|
|
|
1855 |
3rd Qu.:80.00 3rd Qu.:47.75
|
|
|
1856 |
Max. :92.00 Max. :72.00
|
|
|
1857 |
> summary(fm1 <- lm(rating ~ ., data = attitude))
|
|
|
1858 |
|
|
|
1859 |
Call:
|
|
|
1860 |
lm(formula = rating ~ ., data = attitude)
|
|
|
1861 |
|
|
|
1862 |
Residuals:
|
|
|
1863 |
Min 1Q Median 3Q Max
|
|
|
1864 |
-10.9418 -4.3555 0.3158 5.5425 11.5990
|
|
|
1865 |
|
|
|
1866 |
Coefficients:
|
|
|
1867 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1868 |
(Intercept) 10.78708 11.58926 0.931 0.361634
|
|
|
1869 |
complaints 0.61319 0.16098 3.809 0.000903 ***
|
|
|
1870 |
privileges -0.07305 0.13572 -0.538 0.595594
|
|
|
1871 |
learning 0.32033 0.16852 1.901 0.069925 .
|
|
|
1872 |
raises 0.08173 0.22148 0.369 0.715480
|
|
|
1873 |
critical 0.03838 0.14700 0.261 0.796334
|
|
|
1874 |
advance -0.21706 0.17821 -1.218 0.235577
|
|
|
1875 |
---
|
| 61435 |
ripley |
1876 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1877 |
|
|
|
1878 |
Residual standard error: 7.068 on 23 degrees of freedom
|
| 61435 |
ripley |
1879 |
Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
|
|
|
1880 |
F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05
|
| 56186 |
murdoch |
1881 |
|
|
|
1882 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
1883 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
1884 |
> plot(fm1)
|
|
|
1885 |
> summary(fm2 <- lm(rating ~ complaints, data = attitude))
|
|
|
1886 |
|
|
|
1887 |
Call:
|
|
|
1888 |
lm(formula = rating ~ complaints, data = attitude)
|
|
|
1889 |
|
|
|
1890 |
Residuals:
|
|
|
1891 |
Min 1Q Median 3Q Max
|
|
|
1892 |
-12.8799 -5.9905 0.1783 6.2978 9.6294
|
|
|
1893 |
|
|
|
1894 |
Coefficients:
|
|
|
1895 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1896 |
(Intercept) 14.37632 6.61999 2.172 0.0385 *
|
|
|
1897 |
complaints 0.75461 0.09753 7.737 1.99e-08 ***
|
|
|
1898 |
---
|
| 61435 |
ripley |
1899 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1900 |
|
|
|
1901 |
Residual standard error: 6.993 on 28 degrees of freedom
|
| 61435 |
ripley |
1902 |
Multiple R-squared: 0.6813, Adjusted R-squared: 0.6699
|
|
|
1903 |
F-statistic: 59.86 on 1 and 28 DF, p-value: 1.988e-08
|
| 56186 |
murdoch |
1904 |
|
|
|
1905 |
> plot(fm2)
|
|
|
1906 |
> par(opar)
|
|
|
1907 |
>
|
|
|
1908 |
>
|
|
|
1909 |
>
|
|
|
1910 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1911 |
> cleanEx()
|
|
|
1912 |
> nameEx("beavers")
|
|
|
1913 |
> ### * beavers
|
|
|
1914 |
>
|
|
|
1915 |
> flush(stderr()); flush(stdout())
|
|
|
1916 |
>
|
|
|
1917 |
> ### Name: beavers
|
|
|
1918 |
> ### Title: Body Temperature Series of Two Beavers
|
|
|
1919 |
> ### Aliases: beavers beaver1 beaver2
|
|
|
1920 |
> ### Keywords: datasets
|
|
|
1921 |
>
|
|
|
1922 |
> ### ** Examples
|
|
|
1923 |
>
|
|
|
1924 |
> require(graphics)
|
|
|
1925 |
> (yl <- range(beaver1$temp, beaver2$temp))
|
|
|
1926 |
[1] 36.33 38.35
|
|
|
1927 |
>
|
|
|
1928 |
> beaver.plot <- function(bdat, ...) {
|
|
|
1929 |
+ nam <- deparse(substitute(bdat))
|
|
|
1930 |
+ with(bdat, {
|
|
|
1931 |
+ # Hours since start of day:
|
|
|
1932 |
+ hours <- time %/% 100 + 24*(day - day[1]) + (time %% 100)/60
|
|
|
1933 |
+ plot (hours, temp, type = "l", ...,
|
|
|
1934 |
+ main = paste(nam, "body temperature"))
|
|
|
1935 |
+ abline(h = 37.5, col = "gray", lty = 2)
|
|
|
1936 |
+ is.act <- activ == 1
|
|
|
1937 |
+ points(hours[is.act], temp[is.act], col = 2, cex = .8)
|
|
|
1938 |
+ })
|
|
|
1939 |
+ }
|
| 61169 |
ripley |
1940 |
> op <- par(mfrow = c(2, 1), mar = c(3, 3, 4, 2), mgp = 0.9 * 2:0)
|
| 56186 |
murdoch |
1941 |
> beaver.plot(beaver1, ylim = yl)
|
|
|
1942 |
> beaver.plot(beaver2, ylim = yl)
|
|
|
1943 |
> par(op)
|
|
|
1944 |
>
|
|
|
1945 |
>
|
|
|
1946 |
>
|
|
|
1947 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1948 |
> cleanEx()
|
|
|
1949 |
> nameEx("cars")
|
|
|
1950 |
> ### * cars
|
|
|
1951 |
>
|
|
|
1952 |
> flush(stderr()); flush(stdout())
|
|
|
1953 |
>
|
|
|
1954 |
> ### Name: cars
|
|
|
1955 |
> ### Title: Speed and Stopping Distances of Cars
|
|
|
1956 |
> ### Aliases: cars
|
|
|
1957 |
> ### Keywords: datasets
|
|
|
1958 |
>
|
|
|
1959 |
> ### ** Examples
|
|
|
1960 |
>
|
|
|
1961 |
> require(stats); require(graphics)
|
|
|
1962 |
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
|
|
|
1963 |
+ las = 1)
|
|
|
1964 |
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
|
|
|
1965 |
> title(main = "cars data")
|
|
|
1966 |
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
|
|
|
1967 |
+ las = 1, log = "xy")
|
|
|
1968 |
> title(main = "cars data (logarithmic scales)")
|
|
|
1969 |
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
|
|
|
1970 |
> summary(fm1 <- lm(log(dist) ~ log(speed), data = cars))
|
|
|
1971 |
|
|
|
1972 |
Call:
|
|
|
1973 |
lm(formula = log(dist) ~ log(speed), data = cars)
|
|
|
1974 |
|
|
|
1975 |
Residuals:
|
|
|
1976 |
Min 1Q Median 3Q Max
|
|
|
1977 |
-1.00215 -0.24578 -0.02898 0.20717 0.88289
|
|
|
1978 |
|
|
|
1979 |
Coefficients:
|
|
|
1980 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1981 |
(Intercept) -0.7297 0.3758 -1.941 0.0581 .
|
|
|
1982 |
log(speed) 1.6024 0.1395 11.484 2.26e-15 ***
|
|
|
1983 |
---
|
| 61435 |
ripley |
1984 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1985 |
|
|
|
1986 |
Residual standard error: 0.4053 on 48 degrees of freedom
|
| 61435 |
ripley |
1987 |
Multiple R-squared: 0.7331, Adjusted R-squared: 0.7276
|
|
|
1988 |
F-statistic: 131.9 on 1 and 48 DF, p-value: 2.259e-15
|
| 56186 |
murdoch |
1989 |
|
|
|
1990 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
1991 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
1992 |
> plot(fm1)
|
|
|
1993 |
> par(opar)
|
|
|
1994 |
>
|
|
|
1995 |
> ## An example of polynomial regression
|
|
|
1996 |
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
|
|
|
1997 |
+ las = 1, xlim = c(0, 25))
|
|
|
1998 |
> d <- seq(0, 25, length.out = 200)
|
|
|
1999 |
> for(degree in 1:4) {
|
|
|
2000 |
+ fm <- lm(dist ~ poly(speed, degree), data = cars)
|
| 61157 |
ripley |
2001 |
+ assign(paste("cars", degree, sep = "."), fm)
|
|
|
2002 |
+ lines(d, predict(fm, data.frame(speed = d)), col = degree)
|
| 56186 |
murdoch |
2003 |
+ }
|
|
|
2004 |
> anova(cars.1, cars.2, cars.3, cars.4)
|
|
|
2005 |
Analysis of Variance Table
|
|
|
2006 |
|
|
|
2007 |
Model 1: dist ~ poly(speed, degree)
|
|
|
2008 |
Model 2: dist ~ poly(speed, degree)
|
|
|
2009 |
Model 3: dist ~ poly(speed, degree)
|
|
|
2010 |
Model 4: dist ~ poly(speed, degree)
|
|
|
2011 |
Res.Df RSS Df Sum of Sq F Pr(>F)
|
|
|
2012 |
1 48 11354
|
|
|
2013 |
2 47 10825 1 528.81 2.3108 0.1355
|
|
|
2014 |
3 46 10634 1 190.35 0.8318 0.3666
|
|
|
2015 |
4 45 10298 1 336.55 1.4707 0.2316
|
|
|
2016 |
>
|
|
|
2017 |
>
|
|
|
2018 |
>
|
|
|
2019 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
2020 |
> cleanEx()
|
|
|
2021 |
> nameEx("chickwts")
|
|
|
2022 |
> ### * chickwts
|
|
|
2023 |
>
|
|
|
2024 |
> flush(stderr()); flush(stdout())
|
|
|
2025 |
>
|
|
|
2026 |
> ### Name: chickwts
|
|
|
2027 |
> ### Title: Chicken Weights by Feed Type
|
|
|
2028 |
> ### Aliases: chickwts
|
|
|
2029 |
> ### Keywords: datasets
|
|
|
2030 |
>
|
|
|
2031 |
> ### ** Examples
|
|
|
2032 |
>
|
|
|
2033 |
> require(stats); require(graphics)
|
|
|
2034 |
> boxplot(weight ~ feed, data = chickwts, col = "lightgray",
|
|
|
2035 |
+ varwidth = TRUE, notch = TRUE, main = "chickwt data",
|
|
|
2036 |
+ ylab = "Weight at six weeks (gm)")
|
| 78669 |
ripley |
2037 |
Warning in (function (z, notch = FALSE, width = NULL, varwidth = FALSE, :
|
| 56186 |
murdoch |
2038 |
some notches went outside hinges ('box'): maybe set notch=FALSE
|
|
|
2039 |
> anova(fm1 <- lm(weight ~ feed, data = chickwts))
|
|
|
2040 |
Analysis of Variance Table
|
|
|
2041 |
|
|
|
2042 |
Response: weight
|
|
|
2043 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
2044 |
feed 5 231129 46226 15.365 5.936e-10 ***
|
|
|
2045 |
Residuals 65 195556 3009
|
|
|
2046 |
---
|
| 61435 |
ripley |
2047 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2048 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
2049 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
2050 |
> plot(fm1)
|
|
|
2051 |
> par(opar)
|
|
|
2052 |
>
|
|
|
2053 |
>
|
|
|
2054 |
>
|
|
|
2055 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
2056 |
> cleanEx()
|
|
|
2057 |
> nameEx("co2")
|
|
|
2058 |
> ### * co2
|
|
|
2059 |
>
|
|
|
2060 |
> flush(stderr()); flush(stdout())
|
|
|
2061 |
>
|
|
|
2062 |
> ### Name: co2
|
|
|
2063 |
> ### Title: Mauna Loa Atmospheric CO2 Concentration
|
|
|
2064 |
> ### Aliases: co2
|
|
|
2065 |
> ### Keywords: datasets
|
|
|
2066 |
>
|
|
|
2067 |
> ### ** Examples
|
|
|
2068 |
>
|
|
|
2069 |
> require(graphics)
|
|
|
2070 |
> plot(co2, ylab = expression("Atmospheric concentration of CO"[2]),
|
|
|
2071 |
+ las = 1)
|
|
|
2072 |
> title(main = "co2 data set")
|
|
|
2073 |
>
|
|
|
2074 |
>
|
|
|
2075 |
>
|
|
|
2076 |
> cleanEx()
|
|
|
2077 |
> nameEx("crimtab")
|
|
|
2078 |
> ### * crimtab
|
|
|
2079 |
>
|
|
|
2080 |
> flush(stderr()); flush(stdout())
|
|
|
2081 |
>
|
|
|
2082 |
> ### Name: crimtab
|
|
|
2083 |
> ### Title: Student's 3000 Criminals Data
|
|
|
2084 |
> ### Aliases: crimtab
|
|
|
2085 |
> ### Keywords: datasets
|
|
|
2086 |
>
|
|
|
2087 |
> ### ** Examples
|
|
|
2088 |
>
|
|
|
2089 |
> require(stats)
|
|
|
2090 |
> dim(crimtab)
|
|
|
2091 |
[1] 42 22
|
|
|
2092 |
> utils::str(crimtab)
|
|
|
2093 |
'table' int [1:42, 1:22] 0 0 0 0 0 0 1 0 0 0 ...
|
|
|
2094 |
- attr(*, "dimnames")=List of 2
|
|
|
2095 |
..$ : chr [1:42] "9.4" "9.5" "9.6" "9.7" ...
|
|
|
2096 |
..$ : chr [1:22] "142.24" "144.78" "147.32" "149.86" ...
|
|
|
2097 |
> ## for nicer printing:
|
|
|
2098 |
> local({cT <- crimtab
|
| 61169 |
ripley |
2099 |
+ colnames(cT) <- substring(colnames(cT), 2, 3)
|
| 56186 |
murdoch |
2100 |
+ print(cT, zero.print = " ")
|
|
|
2101 |
+ })
|
|
|
2102 |
42 44 47 49 52 54 57 60 62 65 67 70 72 75 77 80 82 85 87 90 93 95
|
|
|
2103 |
9.4
|
|
|
2104 |
9.5 1
|
|
|
2105 |
9.6
|
|
|
2106 |
9.7
|
|
|
2107 |
9.8 1
|
|
|
2108 |
9.9 1 1 1
|
|
|
2109 |
10 1 1 2 2 1
|
|
|
2110 |
10.1 1 3 1 1 1
|
|
|
2111 |
10.2 2 2 2 1 2 1
|
|
|
2112 |
10.3 1 1 3 2 2 3 5
|
|
|
2113 |
10.4 1 1 2 3 3 4 3 3
|
|
|
2114 |
10.5 1 3 7 6 4 3 1 3 1 1
|
|
|
2115 |
10.6 1 4 5 9 14 6 3 1 1
|
|
|
2116 |
10.7 1 2 4 9 14 16 15 7 3 1 2
|
|
|
2117 |
10.8 2 5 6 14 27 10 7 1 2 1
|
|
|
2118 |
10.9 2 6 14 24 27 14 10 4 1
|
|
|
2119 |
11 2 6 12 15 31 37 27 17 10 6
|
|
|
2120 |
11.1 3 3 12 22 26 24 26 24 7 4 1
|
|
|
2121 |
11.2 3 2 7 21 30 38 29 27 20 4 1 1
|
|
|
2122 |
11.3 1 5 10 24 26 39 26 24 7 2
|
|
|
2123 |
11.4 3 4 9 29 56 58 26 22 10 11
|
|
|
2124 |
11.5 5 11 17 33 57 38 34 25 11 2
|
|
|
2125 |
11.6 2 1 4 13 37 39 48 38 27 12 2 2 1
|
|
|
2126 |
11.7 2 9 17 30 37 48 45 24 9 9 2
|
|
|
2127 |
11.8 1 2 11 15 35 41 34 29 10 5 1
|
|
|
2128 |
11.9 1 1 2 12 10 27 32 35 19 10 9 3 1
|
|
|
2129 |
12 1 4 8 19 42 39 22 16 8 2 2
|
|
|
2130 |
12.1 2 4 13 22 28 15 27 10 4 1
|
|
|
2131 |
12.2 1 2 5 6 23 17 16 11 8 1 1
|
|
|
2132 |
12.3 4 8 10 13 20 23 6 5
|
|
|
2133 |
12.4 1 1 1 2 7 12 4 7 7 1 1
|
|
|
2134 |
12.5 1 1 3 12 11 8 6 8 2
|
|
|
2135 |
12.6 1 3 5 7 8 6 3 1 1
|
|
|
2136 |
12.7 1 1 7 5 5 8 2 2
|
|
|
2137 |
12.8 1 2 3 1 8 5 3 1 1
|
|
|
2138 |
12.9 1 2 2 1 1
|
|
|
2139 |
13 3 1 1 2 1
|
|
|
2140 |
13.1 1 1
|
|
|
2141 |
13.2 1 1 1 3
|
|
|
2142 |
13.3 1 1
|
|
|
2143 |
13.4
|
|
|
2144 |
13.5 1
|
|
|
2145 |
>
|
|
|
2146 |
> ## Repeat Student's experiment:
|
|
|
2147 |
>
|
|
|
2148 |
> # 1) Reconstitute 3000 raw data for heights in inches and rounded to
|
|
|
2149 |
> # nearest integer as in Student's paper:
|
|
|
2150 |
>
|
|
|
2151 |
> (heIn <- round(as.numeric(colnames(crimtab)) / 2.54))
|
|
|
2152 |
[1] 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
|
|
|
2153 |
> d.hei <- data.frame(height = rep(heIn, colSums(crimtab)))
|
|
|
2154 |
>
|
|
|
2155 |
> # 2) shuffle the data:
|
|
|
2156 |
>
|
|
|
2157 |
> set.seed(1)
|
|
|
2158 |
> d.hei <- d.hei[sample(1:3000), , drop = FALSE]
|
|
|
2159 |
>
|
|
|
2160 |
> # 3) Make 750 samples each of size 4:
|
|
|
2161 |
>
|
|
|
2162 |
> d.hei$sample <- as.factor(rep(1:750, each = 4))
|
|
|
2163 |
>
|
|
|
2164 |
> # 4) Compute the means and standard deviations (n) for the 750 samples:
|
|
|
2165 |
>
|
|
|
2166 |
> h.mean <- with(d.hei, tapply(height, sample, FUN = mean))
|
|
|
2167 |
> h.sd <- with(d.hei, tapply(height, sample, FUN = sd)) * sqrt(3/4)
|
|
|
2168 |
>
|
|
|
2169 |
> # 5) Compute the difference between the mean of each sample and
|
|
|
2170 |
> # the mean of the population and then divide by the
|
|
|
2171 |
> # standard deviation of the sample:
|
|
|
2172 |
>
|
|
|
2173 |
> zobs <- (h.mean - mean(d.hei[,"height"]))/h.sd
|
|
|
2174 |
>
|
|
|
2175 |
> # 6) Replace infinite values by +/- 6 as in Student's paper:
|
|
|
2176 |
>
|
| 76196 |
ripley |
2177 |
> zobs[infZ <- is.infinite(zobs)] # none of them
|
| 76160 |
luke |
2178 |
named numeric(0)
|
| 56186 |
murdoch |
2179 |
> zobs[infZ] <- 6 * sign(zobs[infZ])
|
|
|
2180 |
>
|
|
|
2181 |
> # 7) Plot the distribution:
|
|
|
2182 |
>
|
|
|
2183 |
> require(grDevices); require(graphics)
|
|
|
2184 |
> hist(x = zobs, probability = TRUE, xlab = "Student's z",
|
|
|
2185 |
+ col = grey(0.8), border = grey(0.5),
|
|
|
2186 |
+ main = "Distribution of Student's z score for 'crimtab' data")
|
|
|
2187 |
>
|
|
|
2188 |
>
|
|
|
2189 |
>
|
|
|
2190 |
> cleanEx()
|
|
|
2191 |
> nameEx("discoveries")
|
|
|
2192 |
> ### * discoveries
|
|
|
2193 |
>
|
|
|
2194 |
> flush(stderr()); flush(stdout())
|
|
|
2195 |
>
|
|
|
2196 |
> ### Name: discoveries
|
|
|
2197 |
> ### Title: Yearly Numbers of Important Discoveries
|
|
|
2198 |
> ### Aliases: discoveries
|
|
|
2199 |
> ### Keywords: datasets
|
|
|
2200 |
>
|
|
|
2201 |
> ### ** Examples
|
|
|
2202 |
>
|
|
|
2203 |
> require(graphics)
|
|
|
2204 |
> plot(discoveries, ylab = "Number of important discoveries",
|
|
|
2205 |
+ las = 1)
|
|
|
2206 |
> title(main = "discoveries data set")
|
|
|
2207 |
>
|
|
|
2208 |
>
|
|
|
2209 |
>
|
|
|
2210 |
> cleanEx()
|
|
|
2211 |
> nameEx("esoph")
|
|
|
2212 |
> ### * esoph
|
|
|
2213 |
>
|
|
|
2214 |
> flush(stderr()); flush(stdout())
|
|
|
2215 |
>
|
|
|
2216 |
> ### Name: esoph
|
|
|
2217 |
> ### Title: Smoking, Alcohol and (O)esophageal Cancer
|
|
|
2218 |
> ### Aliases: esoph
|
|
|
2219 |
> ### Keywords: datasets
|
|
|
2220 |
>
|
|
|
2221 |
> ### ** Examples
|
|
|
2222 |
>
|
|
|
2223 |
> require(stats)
|
|
|
2224 |
> require(graphics) # for mosaicplot
|
|
|
2225 |
> summary(esoph)
|
| 79510 |
hornik |
2226 |
agegp alcgp tobgp ncases ncontrols
|
|
|
2227 |
25-34:15 0-39g/day:23 0-9g/day:24 Min. : 0.000 Min. : 0.000
|
|
|
2228 |
35-44:15 40-79 :23 10-19 :24 1st Qu.: 0.000 1st Qu.: 1.000
|
|
|
2229 |
45-54:16 80-119 :21 20-29 :20 Median : 1.000 Median : 4.000
|
|
|
2230 |
55-64:16 120+ :21 30+ :20 Mean : 2.273 Mean : 8.807
|
|
|
2231 |
65-74:15 3rd Qu.: 4.000 3rd Qu.:10.000
|
|
|
2232 |
75+ :11 Max. :17.000 Max. :60.000
|
| 56186 |
murdoch |
2233 |
> ## effects of alcohol, tobacco and interaction, age-adjusted
|
|
|
2234 |
> model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp,
|
|
|
2235 |
+ data = esoph, family = binomial())
|
|
|
2236 |
> anova(model1)
|
|
|
2237 |
Analysis of Deviance Table
|
|
|
2238 |
|
|
|
2239 |
Model: binomial, link: logit
|
|
|
2240 |
|
|
|
2241 |
Response: cbind(ncases, ncontrols)
|
|
|
2242 |
|
|
|
2243 |
Terms added sequentially (first to last)
|
|
|
2244 |
|
|
|
2245 |
|
| 84104 |
plummer |
2246 |
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
|
|
|
2247 |
NULL 87 367.95
|
|
|
2248 |
agegp 5 121.045 82 246.91 < 2.2e-16 ***
|
|
|
2249 |
tobgp 3 36.639 79 210.27 5.485e-08 ***
|
|
|
2250 |
alcgp 3 127.933 76 82.34 < 2.2e-16 ***
|
|
|
2251 |
tobgp:alcgp 9 5.451 67 76.89 0.7934
|
|
|
2252 |
---
|
|
|
2253 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2254 |
> ## Try a linear effect of alcohol and tobacco
|
|
|
2255 |
> model2 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp)
|
|
|
2256 |
+ + unclass(alcgp),
|
|
|
2257 |
+ data = esoph, family = binomial())
|
|
|
2258 |
> summary(model2)
|
|
|
2259 |
|
|
|
2260 |
Call:
|
|
|
2261 |
glm(formula = cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) +
|
|
|
2262 |
unclass(alcgp), family = binomial(), data = esoph)
|
|
|
2263 |
|
|
|
2264 |
Coefficients:
|
|
|
2265 |
Estimate Std. Error z value Pr(>|z|)
|
| 79510 |
hornik |
2266 |
(Intercept) -4.92067 0.37344 -13.177 < 2e-16 ***
|
|
|
2267 |
agegp.L 3.81892 0.67862 5.627 1.83e-08 ***
|
|
|
2268 |
agegp.Q -1.49473 0.60671 -2.464 0.0138 *
|
|
|
2269 |
agegp.C 0.07923 0.46318 0.171 0.8642
|
|
|
2270 |
agegp^4 0.12136 0.32203 0.377 0.7063
|
|
|
2271 |
agegp^5 -0.24856 0.21153 -1.175 0.2400
|
|
|
2272 |
unclass(tobgp) 0.43955 0.09623 4.568 4.93e-06 ***
|
|
|
2273 |
unclass(alcgp) 1.06766 0.10493 10.175 < 2e-16 ***
|
| 56186 |
murdoch |
2274 |
---
|
| 61435 |
ripley |
2275 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2276 |
|
|
|
2277 |
(Dispersion parameter for binomial family taken to be 1)
|
|
|
2278 |
|
| 79510 |
hornik |
2279 |
Null deviance: 367.953 on 87 degrees of freedom
|
|
|
2280 |
Residual deviance: 91.121 on 80 degrees of freedom
|
|
|
2281 |
AIC: 222.18
|
| 56186 |
murdoch |
2282 |
|
|
|
2283 |
Number of Fisher Scoring iterations: 6
|
|
|
2284 |
|
|
|
2285 |
> ## Re-arrange data for a mosaic plot
|
|
|
2286 |
> ttt <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
|
| 57020 |
murdoch |
2287 |
> o <- with(esoph, order(tobgp, alcgp, agegp))
|
|
|
2288 |
> ttt[ttt == 1] <- esoph$ncases[o]
|
| 56186 |
murdoch |
2289 |
> tt1 <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
|
| 57020 |
murdoch |
2290 |
> tt1[tt1 == 1] <- esoph$ncontrols[o]
|
| 56186 |
murdoch |
2291 |
> tt <- array(c(ttt, tt1), c(dim(ttt),2),
|
|
|
2292 |
+ c(dimnames(ttt), list(c("Cancer", "control"))))
|
|
|
2293 |
> mosaicplot(tt, main = "esoph data set", color = TRUE)
|
|
|
2294 |
>
|
|
|
2295 |
>
|
|
|
2296 |
>
|
|
|
2297 |
> cleanEx()
|
|
|
2298 |
> nameEx("euro")
|
|
|
2299 |
> ### * euro
|
|
|
2300 |
>
|
|
|
2301 |
> flush(stderr()); flush(stdout())
|
|
|
2302 |
>
|
|
|
2303 |
> ### Name: euro
|
|
|
2304 |
> ### Title: Conversion Rates of Euro Currencies
|
|
|
2305 |
> ### Aliases: euro euro.cross
|
|
|
2306 |
> ### Keywords: datasets
|
|
|
2307 |
>
|
|
|
2308 |
> ### ** Examples
|
|
|
2309 |
>
|
|
|
2310 |
> cbind(euro)
|
|
|
2311 |
euro
|
|
|
2312 |
ATS 13.760300
|
|
|
2313 |
BEF 40.339900
|
|
|
2314 |
DEM 1.955830
|
|
|
2315 |
ESP 166.386000
|
|
|
2316 |
FIM 5.945730
|
|
|
2317 |
FRF 6.559570
|
|
|
2318 |
IEP 0.787564
|
|
|
2319 |
ITL 1936.270000
|
|
|
2320 |
LUF 40.339900
|
|
|
2321 |
NLG 2.203710
|
|
|
2322 |
PTE 200.482000
|
|
|
2323 |
>
|
|
|
2324 |
> ## These relations hold:
|
| 61169 |
ripley |
2325 |
> euro == signif(euro, 6) # [6 digit precision in Euro's definition]
|
| 56186 |
murdoch |
2326 |
ATS BEF DEM ESP FIM FRF IEP ITL LUF NLG PTE
|
|
|
2327 |
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
|
|
|
2328 |
> all(euro.cross == outer(1/euro, euro))
|
|
|
2329 |
[1] TRUE
|
|
|
2330 |
>
|
|
|
2331 |
> ## Convert 20 Euro to Belgian Franc
|
|
|
2332 |
> 20 * euro["BEF"]
|
|
|
2333 |
BEF
|
|
|
2334 |
806.798
|
|
|
2335 |
> ## Convert 20 Austrian Schilling to Euro
|
|
|
2336 |
> 20 / euro["ATS"]
|
|
|
2337 |
ATS
|
|
|
2338 |
1.453457
|
|
|
2339 |
> ## Convert 20 Spanish Pesetas to Italian Lira
|
|
|
2340 |
> 20 * euro.cross["ESP", "ITL"]
|
|
|
2341 |
[1] 232.7443
|
|
|
2342 |
>
|
|
|
2343 |
> require(graphics)
|
|
|
2344 |
> dotchart(euro,
|
|
|
2345 |
+ main = "euro data: 1 Euro in currency unit")
|
|
|
2346 |
> dotchart(1/euro,
|
|
|
2347 |
+ main = "euro data: 1 currency unit in Euros")
|
|
|
2348 |
> dotchart(log(euro, 10),
|
|
|
2349 |
+ main = "euro data: log10(1 Euro in currency unit)")
|
|
|
2350 |
>
|
|
|
2351 |
>
|
|
|
2352 |
>
|
|
|
2353 |
> cleanEx()
|
|
|
2354 |
> nameEx("faithful")
|
|
|
2355 |
> ### * faithful
|
|
|
2356 |
>
|
|
|
2357 |
> flush(stderr()); flush(stdout())
|
|
|
2358 |
>
|
|
|
2359 |
> ### Name: faithful
|
|
|
2360 |
> ### Title: Old Faithful Geyser Data
|
|
|
2361 |
> ### Aliases: faithful
|
|
|
2362 |
> ### Keywords: datasets
|
|
|
2363 |
>
|
|
|
2364 |
> ### ** Examples
|
|
|
2365 |
>
|
|
|
2366 |
> require(stats); require(graphics)
|
|
|
2367 |
> f.tit <- "faithful data: Eruptions of Old Faithful"
|
|
|
2368 |
>
|
|
|
2369 |
> ne60 <- round(e60 <- 60 * faithful$eruptions)
|
|
|
2370 |
> all.equal(e60, ne60) # relative diff. ~ 1/10000
|
|
|
2371 |
[1] "Mean relative difference: 9.515332e-05"
|
|
|
2372 |
> table(zapsmall(abs(e60 - ne60))) # 0, 0.02 or 0.04
|
|
|
2373 |
|
|
|
2374 |
|
|
|
2375 |
106 163 3
|
|
|
2376 |
> faithful$better.eruptions <- ne60 / 60
|
|
|
2377 |
> te <- table(ne60)
|
|
|
2378 |
> te[te >= 4] # (too) many multiples of 5 !
|
|
|
2379 |
ne60
|
|
|
2380 |
105 108 110 112 113 120 216 230 240 245 249 250 255 260 261 262 265 270 272 275
|
|
|
2381 |
6 4 7 8 4 4 4 5 6 5 4 4 4 5 4 4 4 8 5 4
|
|
|
2382 |
276 282 288
|
|
|
2383 |
4 6 6
|
| 61157 |
ripley |
2384 |
> plot(names(te), te, type = "h", main = f.tit, xlab = "Eruption time (sec)")
|
| 56186 |
murdoch |
2385 |
>
|
|
|
2386 |
> plot(faithful[, -3], main = f.tit,
|
|
|
2387 |
+ xlab = "Eruption time (min)",
|
|
|
2388 |
+ ylab = "Waiting time to next eruption (min)")
|
|
|
2389 |
> lines(lowess(faithful$eruptions, faithful$waiting, f = 2/3, iter = 3),
|
|
|
2390 |
+ col = "red")
|
|
|
2391 |
>
|
|
|
2392 |
>
|
|
|
2393 |
>
|
|
|
2394 |
> cleanEx()
|
|
|
2395 |
> nameEx("freeny")
|
|
|
2396 |
> ### * freeny
|
|
|
2397 |
>
|
|
|
2398 |
> flush(stderr()); flush(stdout())
|
|
|
2399 |
>
|
|
|
2400 |
> ### Name: freeny
|
|
|
2401 |
> ### Title: Freeny's Revenue Data
|
|
|
2402 |
> ### Aliases: freeny freeny.x freeny.y
|
|
|
2403 |
> ### Keywords: datasets
|
|
|
2404 |
>
|
|
|
2405 |
> ### ** Examples
|
|
|
2406 |
>
|
|
|
2407 |
> require(stats); require(graphics)
|
|
|
2408 |
> summary(freeny)
|
|
|
2409 |
y lag.quarterly.revenue price.index income.level
|
|
|
2410 |
Min. :8.791 Min. :8.791 Min. :4.278 Min. :5.821
|
|
|
2411 |
1st Qu.:9.045 1st Qu.:9.020 1st Qu.:4.392 1st Qu.:5.948
|
|
|
2412 |
Median :9.314 Median :9.284 Median :4.510 Median :6.061
|
|
|
2413 |
Mean :9.306 Mean :9.281 Mean :4.496 Mean :6.039
|
|
|
2414 |
3rd Qu.:9.591 3rd Qu.:9.561 3rd Qu.:4.605 3rd Qu.:6.139
|
|
|
2415 |
Max. :9.794 Max. :9.775 Max. :4.710 Max. :6.200
|
|
|
2416 |
market.potential
|
|
|
2417 |
Min. :12.97
|
|
|
2418 |
1st Qu.:13.01
|
|
|
2419 |
Median :13.07
|
|
|
2420 |
Mean :13.07
|
|
|
2421 |
3rd Qu.:13.12
|
|
|
2422 |
Max. :13.17
|
|
|
2423 |
> pairs(freeny, main = "freeny data")
|
|
|
2424 |
> # gives warning: freeny$y has class "ts"
|
|
|
2425 |
>
|
|
|
2426 |
> summary(fm1 <- lm(y ~ ., data = freeny))
|
|
|
2427 |
|
|
|
2428 |
Call:
|
|
|
2429 |
lm(formula = y ~ ., data = freeny)
|
|
|
2430 |
|
|
|
2431 |
Residuals:
|
|
|
2432 |
Min 1Q Median 3Q Max
|
|
|
2433 |
-0.0259426 -0.0101033 0.0003824 0.0103236 0.0267124
|
|
|
2434 |
|
|
|
2435 |
Coefficients:
|
|
|
2436 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
2437 |
(Intercept) -10.4726 6.0217 -1.739 0.0911 .
|
|
|
2438 |
lag.quarterly.revenue 0.1239 0.1424 0.870 0.3904
|
|
|
2439 |
price.index -0.7542 0.1607 -4.693 4.28e-05 ***
|
|
|
2440 |
income.level 0.7675 0.1339 5.730 1.93e-06 ***
|
|
|
2441 |
market.potential 1.3306 0.5093 2.613 0.0133 *
|
|
|
2442 |
---
|
| 61435 |
ripley |
2443 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2444 |
|
|
|
2445 |
Residual standard error: 0.01473 on 34 degrees of freedom
|
| 61435 |
ripley |
2446 |
Multiple R-squared: 0.9981, Adjusted R-squared: 0.9978
|
|
|
2447 |
F-statistic: 4354 on 4 and 34 DF, p-value: < 2.2e-16
|
| 56186 |
murdoch |
2448 |
|
|
|
2449 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
2450 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
2451 |
> plot(fm1)
|
|
|
2452 |
> par(opar)
|
|
|
2453 |
>
|
|
|
2454 |
>
|
|
|
2455 |
>
|
|
|
2456 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
2457 |
> cleanEx()
|
| 87723 |
maechler |
2458 |
> nameEx("gait")
|
|
|
2459 |
> ### * gait
|
|
|
2460 |
>
|
|
|
2461 |
> flush(stderr()); flush(stdout())
|
|
|
2462 |
>
|
|
|
2463 |
> ### Name: gait
|
|
|
2464 |
> ### Title: Hip and Knee Angle while Walking
|
|
|
2465 |
> ### Aliases: gait
|
|
|
2466 |
> ### Keywords: datasets
|
|
|
2467 |
>
|
|
|
2468 |
> ### ** Examples
|
|
|
2469 |
>
|
|
|
2470 |
> plot(gait[, 1, ], type = "b",
|
|
|
2471 |
+ xlim = range(gait[,,1]), ylim = range(gait[,,2]),
|
|
|
2472 |
+ xlab = "Hip Angle", ylab = "Knee Angle", main = "'gait' data : Boy 1")
|
|
|
2473 |
> mtext("all other boys", col = "thistle"); grid()
|
|
|
2474 |
> matlines(gait[, -1, 1], gait[, -1, 2], type = "l", lty = 1, col = adjustcolor("thistle", 1/3))
|
|
|
2475 |
>
|
| 87725 |
maechler |
2476 |
> ## The data array, two matrices :
|
|
|
2477 |
> op <- options(width = 128) # on a wide console
|
|
|
2478 |
> aperm(gait, c(2:1, 3))
|
|
|
2479 |
, , Variable = Hip Angle
|
|
|
2480 |
|
|
|
2481 |
Time
|
|
|
2482 |
Subject 0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
|
|
|
2483 |
boy1 37 36 33 29 23 18 15 12 9 6 4 6 13 22 31 38 43 44 40 35
|
|
|
2484 |
boy2 47 46 42 34 27 21 16 12 4 -1 -2 2 12 24 34 42 47 48 45 43
|
|
|
2485 |
boy3 46 44 39 34 33 27 21 15 7 1 -3 -1 10 22 34 43 51 53 50 49
|
|
|
2486 |
boy4 37 36 27 20 15 15 10 6 3 0 3 4 10 14 23 31 33 35 34 32
|
|
|
2487 |
boy5 20 18 11 8 7 5 1 0 -2 -4 -4 1 10 18 27 36 34 33 32 28
|
|
|
2488 |
boy6 57 48 44 35 31 27 21 18 14 10 8 8 11 21 32 40 49 55 56 55
|
|
|
2489 |
boy7 46 38 33 25 18 15 14 10 8 4 0 3 9 19 31 37 46 45 46 41
|
|
|
2490 |
boy8 46 46 43 40 36 30 24 17 12 7 1 3 12 24 35 44 49 49 47 43
|
|
|
2491 |
boy9 46 42 37 34 31 25 19 17 13 10 5 7 13 21 29 39 47 51 52 47
|
|
|
2492 |
boy10 35 34 29 28 19 15 13 6 0 -1 2 11 20 30 39 45 48 48 43 39
|
|
|
2493 |
boy11 38 37 33 29 26 20 14 9 3 -2 -1 5 15 24 36 44 47 45 42 40
|
|
|
2494 |
boy12 35 31 29 26 22 19 15 11 6 3 0 4 16 27 35 41 44 45 40 39
|
|
|
2495 |
boy13 34 31 27 23 19 15 13 10 6 1 0 3 8 16 27 34 40 41 39 36
|
|
|
2496 |
boy14 43 41 36 31 26 20 16 13 10 8 4 12 19 27 37 41 44 44 41 37
|
|
|
2497 |
boy15 43 37 35 28 26 21 20 15 13 7 3 8 16 25 35 43 49 50 45 46
|
|
|
2498 |
boy16 40 41 36 32 27 20 17 10 5 -2 -4 -2 7 24 38 49 57 59 54 46
|
|
|
2499 |
boy17 51 49 45 39 31 23 16 9 4 0 -5 1 10 21 33 44 51 55 56 51
|
|
|
2500 |
boy18 52 46 41 35 31 24 18 12 6 3 2 5 12 22 32 42 46 48 49 46
|
|
|
2501 |
boy19 36 33 28 22 18 13 11 7 1 -2 -3 -3 2 13 22 32 39 41 38 34
|
|
|
2502 |
boy20 35 37 33 27 22 14 9 6 4 -1 1 2 11 17 34 35 40 43 43 42
|
|
|
2503 |
boy21 46 38 30 23 17 13 7 3 0 -3 -3 1 11 20 36 48 55 57 56 50
|
|
|
2504 |
boy22 43 41 37 30 24 16 12 8 5 2 1 8 16 24 33 42 48 48 48 46
|
|
|
2505 |
boy23 55 51 47 41 35 30 26 23 20 11 8 10 19 35 41 52 57 61 63 58
|
|
|
2506 |
boy24 39 38 31 27 21 14 9 9 6 4 2 3 13 30 37 47 53 53 49 44
|
|
|
2507 |
boy25 37 34 30 27 26 19 15 10 4 0 0 4 12 22 31 40 43 43 38 38
|
|
|
2508 |
boy26 36 33 28 22 18 13 11 7 1 -2 -3 -3 2 13 22 32 39 41 38 34
|
|
|
2509 |
boy27 36 33 30 28 21 15 8 1 -5 -11 -12 -7 4 16 26 37 44 47 44 37
|
|
|
2510 |
boy28 42 40 40 34 23 15 11 7 5 6 8 12 22 33 43 51 57 58 54 46
|
|
|
2511 |
boy29 38 34 30 23 17 12 8 4 1 0 -4 -4 1 10 22 32 38 41 41 40
|
|
|
2512 |
boy30 46 47 44 37 29 23 19 14 8 3 -2 1 12 26 39 48 52 48 43 42
|
|
|
2513 |
boy31 54 48 44 37 30 27 21 18 15 11 12 20 30 41 50 56 61 59 57 58
|
|
|
2514 |
boy32 52 44 44 33 28 27 23 24 19 15 15 16 24 32 43 52 58 59 57 52
|
|
|
2515 |
boy33 32 28 26 22 19 13 8 5 -1 -6 -5 0 12 22 30 36 39 36 30 29
|
|
|
2516 |
boy34 46 41 38 31 25 20 13 7 1 -4 -3 -2 8 22 34 45 53 57 55 43
|
|
|
2517 |
boy35 46 44 40 35 31 25 19 15 10 5 3 6 14 27 40 48 53 53 50 47
|
|
|
2518 |
boy36 48 42 42 35 30 23 19 14 9 4 -1 3 9 25 37 45 52 53 52 46
|
|
|
2519 |
boy37 44 41 38 32 24 18 10 6 3 0 -2 -1 6 20 36 44 49 46 38 35
|
|
|
2520 |
boy38 55 56 51 46 41 36 30 25 21 15 8 5 9 19 31 43 52 56 59 59
|
|
|
2521 |
boy39 48 50 47 42 37 29 22 14 8 5 8 15 24 36 51 59 63 64 61 55
|
|
|
2522 |
|
|
|
2523 |
, , Variable = Knee Angle
|
|
|
2524 |
|
|
|
2525 |
Time
|
|
|
2526 |
Subject 0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
|
|
|
2527 |
boy1 10 15 18 18 15 14 13 13 15 21 26 39 52 64 70 66 57 40 22 11
|
|
|
2528 |
boy2 16 25 28 25 18 12 10 9 7 9 17 32 51 67 70 66 55 39 23 16
|
|
|
2529 |
boy3 18 27 32 32 28 23 16 11 8 11 18 33 56 72 79 77 67 46 22 14
|
|
|
2530 |
boy4 5 14 16 17 10 8 6 5 5 9 20 40 57 66 67 61 43 18 3 0
|
|
|
2531 |
boy5 2 6 6 6 5 6 6 7 8 15 24 41 59 69 71 66 44 23 5 4
|
|
|
2532 |
boy6 15 17 23 23 20 19 15 13 9 10 14 25 40 58 65 66 61 45 21 9
|
|
|
2533 |
boy7 13 16 22 17 12 9 8 10 14 16 22 36 51 64 71 68 65 40 17 3
|
|
|
2534 |
boy8 14 17 18 19 19 15 11 8 7 8 12 26 47 68 79 76 61 35 13 5
|
|
|
2535 |
boy9 15 20 23 26 25 21 16 13 8 8 10 22 40 56 65 67 61 47 26 10
|
|
|
2536 |
boy10 9 22 25 21 10 9 9 9 11 19 33 55 69 75 77 70 55 36 17 16
|
|
|
2537 |
boy11 13 24 27 23 18 13 12 15 16 19 27 43 60 70 75 69 50 25 13 6
|
|
|
2538 |
boy12 7 8 11 12 8 6 4 5 5 10 19 36 58 70 72 66 51 30 8 5
|
|
|
2539 |
boy13 9 14 16 15 15 12 13 14 14 15 20 31 45 60 71 70 62 45 26 12
|
|
|
2540 |
boy14 15 20 22 22 21 19 19 19 19 20 25 41 58 70 73 68 53 32 16 10
|
|
|
2541 |
boy15 6 11 20 18 13 9 9 8 11 12 20 38 56 67 71 65 48 25 5 8
|
|
|
2542 |
boy16 11 19 30 28 25 17 15 10 7 7 13 25 46 68 79 77 67 48 23 13
|
|
|
2543 |
boy17 24 32 35 33 29 24 18 14 12 14 20 36 55 69 75 72 60 43 26 16
|
|
|
2544 |
boy18 16 20 21 20 18 14 13 11 9 11 19 34 53 66 73 72 61 43 22 11
|
|
|
2545 |
boy19 16 20 22 21 20 20 20 19 16 16 21 31 49 69 74 71 60 41 23 15
|
|
|
2546 |
boy20 7 13 14 17 14 8 7 7 7 8 17 30 49 60 71 65 58 39 20 10
|
|
|
2547 |
boy21 21 24 25 21 16 9 4 4 5 9 19 36 56 70 80 80 68 48 27 16
|
|
|
2548 |
boy22 11 14 14 11 8 5 6 7 8 10 18 36 53 65 72 73 64 46 28 13
|
|
|
2549 |
boy23 12 17 20 20 18 12 11 9 11 10 16 30 49 68 69 69 63 44 28 8
|
|
|
2550 |
boy24 8 12 14 13 12 9 7 10 9 13 21 39 59 76 76 71 55 26 8 12
|
|
|
2551 |
boy25 11 20 22 21 21 17 16 17 15 19 28 42 59 70 73 66 49 24 6 9
|
|
|
2552 |
boy26 16 20 22 21 20 20 20 19 16 16 21 31 49 69 74 71 60 41 23 15
|
|
|
2553 |
boy27 19 26 28 28 24 18 12 8 7 8 16 33 52 68 78 80 74 57 35 18
|
|
|
2554 |
boy28 13 23 30 28 19 10 6 5 10 19 33 50 68 79 81 75 61 42 25 16
|
|
|
2555 |
boy29 11 15 19 20 18 17 15 11 9 9 12 23 40 58 71 70 56 36 19 11
|
|
|
2556 |
boy30 17 25 30 30 27 22 17 11 7 7 13 31 52 69 75 68 46 21 6 8
|
|
|
2557 |
boy31 20 20 22 16 10 10 10 12 13 17 28 47 64 72 72 63 47 23 9 19
|
|
|
2558 |
boy32 18 18 25 23 18 19 21 25 24 25 27 39 59 72 77 69 60 38 17 14
|
|
|
2559 |
boy33 9 12 16 15 14 11 8 6 5 9 18 39 56 68 71 67 54 30 9 6
|
|
|
2560 |
boy34 8 10 17 16 12 10 8 4 2 4 12 27 48 66 76 77 71 54 26 3
|
|
|
2561 |
boy35 9 18 19 19 19 15 12 11 10 12 20 36 56 72 80 75 59 34 10 9
|
|
|
2562 |
boy36 13 18 27 26 25 18 12 10 9 13 16 33 53 73 77 75 64 44 23 16
|
|
|
2563 |
boy37 19 23 26 25 21 18 14 13 14 17 24 37 56 73 82 76 59 31 7 7
|
|
|
2564 |
boy38 16 23 28 28 25 21 19 17 16 13 13 19 37 57 69 71 62 45 28 20
|
|
|
2565 |
boy39 14 25 32 34 30 20 12 7 6 11 24 45 67 79 82 76 65 46 25 15
|
|
|
2566 |
|
|
|
2567 |
> options(op)
|
| 87723 |
maechler |
2568 |
>
|
|
|
2569 |
>
|
| 87725 |
maechler |
2570 |
>
|
| 87723 |
maechler |
2571 |
> cleanEx()
|
| 56186 |
murdoch |
2572 |
> nameEx("infert")
|
|
|
2573 |
> ### * infert
|
|
|
2574 |
>
|
|
|
2575 |
> flush(stderr()); flush(stdout())
|
|
|
2576 |
>
|
|
|
2577 |
> ### Name: infert
|
|
|
2578 |
> ### Title: Infertility after Spontaneous and Induced Abortion
|
|
|
2579 |
> ### Aliases: infert
|
|
|
2580 |
> ### Keywords: datasets
|
|
|
2581 |
>
|
|
|
2582 |
> ### ** Examples
|
|
|
2583 |
>
|
|
|
2584 |
> require(stats)
|
| 61157 |
ripley |
2585 |
> model1 <- glm(case ~ spontaneous+induced, data = infert, family = binomial())
|
| 56186 |
murdoch |
2586 |
> summary(model1)
|
|
|
2587 |
|
|
|
2588 |
Call:
|
|
|
2589 |
glm(formula = case ~ spontaneous + induced, family = binomial(),
|
|
|
2590 |
data = infert)
|
|
|
2591 |
|
|
|
2592 |
Coefficients:
|
|
|
2593 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
2594 |
(Intercept) -1.7079 0.2677 -6.380 1.78e-10 ***
|
|
|
2595 |
spontaneous 1.1972 0.2116 5.657 1.54e-08 ***
|
|
|
2596 |
induced 0.4181 0.2056 2.033 0.042 *
|
|
|
2597 |
---
|
| 61435 |
ripley |
2598 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2599 |
|
|
|
2600 |
(Dispersion parameter for binomial family taken to be 1)
|
|
|
2601 |
|
|
|
2602 |
Null deviance: 316.17 on 247 degrees of freedom
|
|
|
2603 |
Residual deviance: 279.61 on 245 degrees of freedom
|
|
|
2604 |
AIC: 285.61
|
|
|
2605 |
|
|
|
2606 |
Number of Fisher Scoring iterations: 4
|
|
|
2607 |
|
|
|
2608 |
> ## adjusted for other potential confounders:
|
|
|
2609 |
> summary(model2 <- glm(case ~ age+parity+education+spontaneous+induced,
|
| 61157 |
ripley |
2610 |
+ data = infert, family = binomial()))
|
| 56186 |
murdoch |
2611 |
|
|
|
2612 |
Call:
|
|
|
2613 |
glm(formula = case ~ age + parity + education + spontaneous +
|
|
|
2614 |
induced, family = binomial(), data = infert)
|
|
|
2615 |
|
|
|
2616 |
Coefficients:
|
|
|
2617 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
2618 |
(Intercept) -1.14924 1.41220 -0.814 0.4158
|
|
|
2619 |
age 0.03958 0.03120 1.269 0.2046
|
|
|
2620 |
parity -0.82828 0.19649 -4.215 2.49e-05 ***
|
|
|
2621 |
education6-11yrs -1.04424 0.79255 -1.318 0.1876
|
|
|
2622 |
education12+ yrs -1.40321 0.83416 -1.682 0.0925 .
|
|
|
2623 |
spontaneous 2.04591 0.31016 6.596 4.21e-11 ***
|
|
|
2624 |
induced 1.28876 0.30146 4.275 1.91e-05 ***
|
|
|
2625 |
---
|
| 61435 |
ripley |
2626 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2627 |
|
|
|
2628 |
(Dispersion parameter for binomial family taken to be 1)
|
|
|
2629 |
|
|
|
2630 |
Null deviance: 316.17 on 247 degrees of freedom
|
|
|
2631 |
Residual deviance: 257.80 on 241 degrees of freedom
|
|
|
2632 |
AIC: 271.8
|
|
|
2633 |
|
|
|
2634 |
Number of Fisher Scoring iterations: 4
|
|
|
2635 |
|
|
|
2636 |
> ## Really should be analysed by conditional logistic regression
|
|
|
2637 |
> ## which is in the survival package
|
|
|
2638 |
>
|
|
|
2639 |
>
|
|
|
2640 |
>
|
|
|
2641 |
> cleanEx()
|
|
|
2642 |
> nameEx("iris")
|
|
|
2643 |
> ### * iris
|
|
|
2644 |
>
|
|
|
2645 |
> flush(stderr()); flush(stdout())
|
|
|
2646 |
>
|
|
|
2647 |
> ### Name: iris
|
|
|
2648 |
> ### Title: Edgar Anderson's Iris Data
|
|
|
2649 |
> ### Aliases: iris iris3
|
|
|
2650 |
> ### Keywords: datasets
|
|
|
2651 |
>
|
|
|
2652 |
> ### ** Examples
|
|
|
2653 |
>
|
| 85134 |
smeyer |
2654 |
> summary(iris)
|
|
|
2655 |
Sepal.Length Sepal.Width Petal.Length Petal.Width
|
|
|
2656 |
Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
|
|
|
2657 |
1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
|
|
|
2658 |
Median :5.800 Median :3.000 Median :4.350 Median :1.300
|
|
|
2659 |
Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
|
|
|
2660 |
3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
|
|
|
2661 |
Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
|
|
|
2662 |
Species
|
|
|
2663 |
setosa :50
|
|
|
2664 |
versicolor:50
|
|
|
2665 |
virginica :50
|
|
|
2666 |
|
|
|
2667 |
|
|
|
2668 |
|
|
|
2669 |
>
|
|
|
2670 |
> ## Fisher's (1936) research question: whether (compound measurements of)
|
|
|
2671 |
> ## Iris versicolor "differs twice as much from I. setosa as from I. virginica"
|
|
|
2672 |
> pairs(iris[1:4], col = iris$Species)
|
|
|
2673 |
> legend(0.5, 1, levels(iris$Species), fill = 1:3, bty = "n",
|
|
|
2674 |
+ horiz = TRUE, xjust = 0.5, yjust = 0, xpd = TRUE)
|
|
|
2675 |
>
|
|
|
2676 |
> ## equivalence of legacy array (iris3) and data.frame (iris) representation
|
| 56186 |
murdoch |
2677 |
> dni3 <- dimnames(iris3)
|
| 61157 |
ripley |
2678 |
> ii <- data.frame(matrix(aperm(iris3, c(1,3,2)), ncol = 4,
|
| 56186 |
murdoch |
2679 |
+ dimnames = list(NULL, sub(" L.",".Length",
|
|
|
2680 |
+ sub(" W.",".Width", dni3[[2]])))),
|
| 61157 |
ripley |
2681 |
+ Species = gl(3, 50, labels = sub("S", "s", sub("V", "v", dni3[[3]]))))
|
| 85134 |
smeyer |
2682 |
> stopifnot(all.equal(ii, iris))
|
| 56186 |
murdoch |
2683 |
>
|
|
|
2684 |
>
|
|
|
2685 |
>
|
|
|
2686 |
> cleanEx()
|
|
|
2687 |
> nameEx("islands")
|
|
|
2688 |
> ### * islands
|
|
|
2689 |
>
|
|
|
2690 |
> flush(stderr()); flush(stdout())
|
|
|
2691 |
>
|
|
|
2692 |
> ### Name: islands
|
|
|
2693 |
> ### Title: Areas of the World's Major Landmasses
|
|
|
2694 |
> ### Aliases: islands
|
|
|
2695 |
> ### Keywords: datasets
|
|
|
2696 |
>
|
|
|
2697 |
> ### ** Examples
|
|
|
2698 |
>
|
|
|
2699 |
> require(graphics)
|
|
|
2700 |
> dotchart(log(islands, 10),
|
|
|
2701 |
+ main = "islands data: log10(area) (log10(sq. miles))")
|
|
|
2702 |
> dotchart(log(islands[order(islands)], 10),
|
|
|
2703 |
+ main = "islands data: log10(area) (log10(sq. miles))")
|
|
|
2704 |
>
|
|
|
2705 |
>
|
|
|
2706 |
>
|
|
|
2707 |
> cleanEx()
|
|
|
2708 |
> nameEx("longley")
|
|
|
2709 |
> ### * longley
|
|
|
2710 |
>
|
|
|
2711 |
> flush(stderr()); flush(stdout())
|
|
|
2712 |
>
|
|
|
2713 |
> ### Name: longley
|
|
|
2714 |
> ### Title: Longley's Economic Regression Data
|
|
|
2715 |
> ### Aliases: longley
|
|
|
2716 |
> ### Keywords: datasets
|
|
|
2717 |
>
|
|
|
2718 |
> ### ** Examples
|
|
|
2719 |
>
|
|
|
2720 |
> require(stats); require(graphics)
|
| 85836 |
ripley |
2721 |
> ## give the data set in the form it was used in S-PLUS:
|
| 56186 |
murdoch |
2722 |
> longley.x <- data.matrix(longley[, 1:6])
|
|
|
2723 |
> longley.y <- longley[, "Employed"]
|
|
|
2724 |
> pairs(longley, main = "longley data")
|
|
|
2725 |
> summary(fm1 <- lm(Employed ~ ., data = longley))
|
|
|
2726 |
|
|
|
2727 |
Call:
|
|
|
2728 |
lm(formula = Employed ~ ., data = longley)
|
|
|
2729 |
|
|
|
2730 |
Residuals:
|
|
|
2731 |
Min 1Q Median 3Q Max
|
|
|
2732 |
-0.41011 -0.15767 -0.02816 0.10155 0.45539
|
|
|
2733 |
|
|
|
2734 |
Coefficients:
|
|
|
2735 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
2736 |
(Intercept) -3.482e+03 8.904e+02 -3.911 0.003560 **
|
|
|
2737 |
GNP.deflator 1.506e-02 8.492e-02 0.177 0.863141
|
|
|
2738 |
GNP -3.582e-02 3.349e-02 -1.070 0.312681
|
|
|
2739 |
Unemployed -2.020e-02 4.884e-03 -4.136 0.002535 **
|
|
|
2740 |
Armed.Forces -1.033e-02 2.143e-03 -4.822 0.000944 ***
|
|
|
2741 |
Population -5.110e-02 2.261e-01 -0.226 0.826212
|
|
|
2742 |
Year 1.829e+00 4.555e-01 4.016 0.003037 **
|
|
|
2743 |
---
|
| 61435 |
ripley |
2744 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2745 |
|
|
|
2746 |
Residual standard error: 0.3049 on 9 degrees of freedom
|
| 61435 |
ripley |
2747 |
Multiple R-squared: 0.9955, Adjusted R-squared: 0.9925
|
|
|
2748 |
F-statistic: 330.3 on 6 and 9 DF, p-value: 4.984e-10
|
| 56186 |
murdoch |
2749 |
|
|
|
2750 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
2751 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
2752 |
> plot(fm1)
|
|
|
2753 |
> par(opar)
|
|
|
2754 |
>
|
|
|
2755 |
>
|
|
|
2756 |
>
|
|
|
2757 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
2758 |
> cleanEx()
|
|
|
2759 |
> nameEx("morley")
|
|
|
2760 |
> ### * morley
|
|
|
2761 |
>
|
|
|
2762 |
> flush(stderr()); flush(stdout())
|
|
|
2763 |
>
|
|
|
2764 |
> ### Name: morley
|
| 59017 |
ripley |
2765 |
> ### Title: Michelson Speed of Light Data
|
| 56186 |
murdoch |
2766 |
> ### Aliases: morley
|
|
|
2767 |
> ### Keywords: datasets
|
|
|
2768 |
>
|
|
|
2769 |
> ### ** Examples
|
|
|
2770 |
>
|
|
|
2771 |
> require(stats); require(graphics)
|
| 58992 |
maechler |
2772 |
> michelson <- transform(morley,
|
|
|
2773 |
+ Expt = factor(Expt), Run = factor(Run))
|
| 61169 |
ripley |
2774 |
> xtabs(~ Expt + Run, data = michelson) # 5 x 20 balanced (two-way)
|
| 56186 |
murdoch |
2775 |
Run
|
|
|
2776 |
Expt 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
|
|
|
2777 |
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2778 |
2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2779 |
3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2780 |
4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2781 |
5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
| 58992 |
maechler |
2782 |
> plot(Speed ~ Expt, data = michelson,
|
| 56186 |
murdoch |
2783 |
+ main = "Speed of Light Data", xlab = "Experiment No.")
|
| 58992 |
maechler |
2784 |
> fm <- aov(Speed ~ Run + Expt, data = michelson)
|
| 56186 |
murdoch |
2785 |
> summary(fm)
|
| 57044 |
ripley |
2786 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
2787 |
Run 19 113344 5965 1.105 0.36321
|
|
|
2788 |
Expt 4 94514 23629 4.378 0.00307 **
|
|
|
2789 |
Residuals 76 410166 5397
|
| 56186 |
murdoch |
2790 |
---
|
| 61435 |
ripley |
2791 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2792 |
> fm0 <- update(fm, . ~ . - Run)
|
|
|
2793 |
> anova(fm0, fm)
|
|
|
2794 |
Analysis of Variance Table
|
|
|
2795 |
|
|
|
2796 |
Model 1: Speed ~ Expt
|
|
|
2797 |
Model 2: Speed ~ Run + Expt
|
|
|
2798 |
Res.Df RSS Df Sum of Sq F Pr(>F)
|
|
|
2799 |
1 95 523510
|
|
|
2800 |
2 76 410166 19 113344 1.1053 0.3632
|
|
|
2801 |
>
|
|
|
2802 |
>
|
|
|
2803 |
>
|
|
|
2804 |
> cleanEx()
|
|
|
2805 |
> nameEx("mtcars")
|
|
|
2806 |
> ### * mtcars
|
|
|
2807 |
>
|
|
|
2808 |
> flush(stderr()); flush(stdout())
|
|
|
2809 |
>
|
|
|
2810 |
> ### Name: mtcars
|
|
|
2811 |
> ### Title: Motor Trend Car Road Tests
|
|
|
2812 |
> ### Aliases: mtcars
|
|
|
2813 |
> ### Keywords: datasets
|
|
|
2814 |
>
|
|
|
2815 |
> ### ** Examples
|
|
|
2816 |
>
|
|
|
2817 |
> require(graphics)
|
| 74582 |
maechler |
2818 |
> pairs(mtcars, main = "mtcars data", gap = 1/4)
|
| 56186 |
murdoch |
2819 |
> coplot(mpg ~ disp | as.factor(cyl), data = mtcars,
|
|
|
2820 |
+ panel = panel.smooth, rows = 1)
|
| 74582 |
maechler |
2821 |
> ## possibly more meaningful, e.g., for summary() or bivariate plots:
|
|
|
2822 |
> mtcars2 <- within(mtcars, {
|
|
|
2823 |
+ vs <- factor(vs, labels = c("V", "S"))
|
| 74608 |
maechler |
2824 |
+ am <- factor(am, labels = c("automatic", "manual"))
|
| 74582 |
maechler |
2825 |
+ cyl <- ordered(cyl)
|
|
|
2826 |
+ gear <- ordered(gear)
|
|
|
2827 |
+ carb <- ordered(carb)
|
|
|
2828 |
+ })
|
|
|
2829 |
> summary(mtcars2)
|
|
|
2830 |
mpg cyl disp hp drat
|
|
|
2831 |
Min. :10.40 4:11 Min. : 71.1 Min. : 52.0 Min. :2.760
|
|
|
2832 |
1st Qu.:15.43 6: 7 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080
|
|
|
2833 |
Median :19.20 8:14 Median :196.3 Median :123.0 Median :3.695
|
|
|
2834 |
Mean :20.09 Mean :230.7 Mean :146.7 Mean :3.597
|
|
|
2835 |
3rd Qu.:22.80 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920
|
|
|
2836 |
Max. :33.90 Max. :472.0 Max. :335.0 Max. :4.930
|
|
|
2837 |
wt qsec vs am gear carb
|
| 74608 |
maechler |
2838 |
Min. :1.513 Min. :14.50 V:18 automatic:19 3:15 1: 7
|
|
|
2839 |
1st Qu.:2.581 1st Qu.:16.89 S:14 manual :13 4:12 2:10
|
| 74582 |
maechler |
2840 |
Median :3.325 Median :17.71 5: 5 3: 3
|
|
|
2841 |
Mean :3.217 Mean :17.85 4:10
|
|
|
2842 |
3rd Qu.:3.610 3rd Qu.:18.90 6: 1
|
|
|
2843 |
Max. :5.424 Max. :22.90 8: 1
|
| 56186 |
murdoch |
2844 |
>
|
|
|
2845 |
>
|
|
|
2846 |
>
|
|
|
2847 |
> cleanEx()
|
|
|
2848 |
> nameEx("nhtemp")
|
|
|
2849 |
> ### * nhtemp
|
|
|
2850 |
>
|
|
|
2851 |
> flush(stderr()); flush(stdout())
|
|
|
2852 |
>
|
|
|
2853 |
> ### Name: nhtemp
|
|
|
2854 |
> ### Title: Average Yearly Temperatures in New Haven
|
|
|
2855 |
> ### Aliases: nhtemp
|
|
|
2856 |
> ### Keywords: datasets
|
|
|
2857 |
>
|
|
|
2858 |
> ### ** Examples
|
|
|
2859 |
>
|
|
|
2860 |
> require(stats); require(graphics)
|
|
|
2861 |
> plot(nhtemp, main = "nhtemp data",
|
|
|
2862 |
+ ylab = "Mean annual temperature in New Haven, CT (deg. F)")
|
|
|
2863 |
>
|
|
|
2864 |
>
|
|
|
2865 |
>
|
|
|
2866 |
> cleanEx()
|
|
|
2867 |
> nameEx("nottem")
|
|
|
2868 |
> ### * nottem
|
|
|
2869 |
>
|
|
|
2870 |
> flush(stderr()); flush(stdout())
|
|
|
2871 |
>
|
|
|
2872 |
> ### Name: nottem
|
|
|
2873 |
> ### Title: Average Monthly Temperatures at Nottingham, 1920-1939
|
|
|
2874 |
> ### Aliases: nottem
|
|
|
2875 |
> ### Keywords: datasets
|
|
|
2876 |
>
|
|
|
2877 |
> ### ** Examples
|
|
|
2878 |
>
|
| 86027 |
smeyer |
2879 |
> require(stats); require(graphics)
|
|
|
2880 |
> nott <- window(nottem, end = c(1936,12))
|
|
|
2881 |
> fit <- arima(nott, order = c(1,0,0), list(order = c(2,1,0), period = 12))
|
|
|
2882 |
> nott.fore <- predict(fit, n.ahead = 36)
|
|
|
2883 |
> ts.plot(nott, nott.fore$pred, nott.fore$pred+2*nott.fore$se,
|
|
|
2884 |
+ nott.fore$pred-2*nott.fore$se, gpars = list(col = c(1,1,4,4)))
|
| 56186 |
murdoch |
2885 |
>
|
|
|
2886 |
>
|
| 86027 |
smeyer |
2887 |
>
|
| 56186 |
murdoch |
2888 |
> cleanEx()
|
| 63381 |
ripley |
2889 |
> nameEx("npk")
|
|
|
2890 |
> ### * npk
|
|
|
2891 |
>
|
|
|
2892 |
> flush(stderr()); flush(stdout())
|
|
|
2893 |
>
|
|
|
2894 |
> ### Name: npk
|
|
|
2895 |
> ### Title: Classical N, P, K Factorial Experiment
|
|
|
2896 |
> ### Aliases: npk
|
|
|
2897 |
> ### Keywords: datasets
|
|
|
2898 |
>
|
|
|
2899 |
> ### ** Examples
|
|
|
2900 |
>
|
|
|
2901 |
>
|
|
|
2902 |
> cleanEx()
|
| 56186 |
murdoch |
2903 |
> nameEx("occupationalStatus")
|
|
|
2904 |
> ### * occupationalStatus
|
|
|
2905 |
>
|
|
|
2906 |
> flush(stderr()); flush(stdout())
|
|
|
2907 |
>
|
|
|
2908 |
> ### Name: occupationalStatus
|
|
|
2909 |
> ### Title: Occupational Status of Fathers and their Sons
|
|
|
2910 |
> ### Aliases: occupationalStatus
|
|
|
2911 |
> ### Keywords: datasets
|
|
|
2912 |
>
|
|
|
2913 |
> ### ** Examples
|
|
|
2914 |
>
|
|
|
2915 |
> require(stats); require(graphics)
|
|
|
2916 |
>
|
|
|
2917 |
> plot(occupationalStatus)
|
|
|
2918 |
>
|
|
|
2919 |
> ## Fit a uniform association model separating diagonal effects
|
|
|
2920 |
> Diag <- as.factor(diag(1:8))
|
|
|
2921 |
> Rscore <- scale(as.numeric(row(occupationalStatus)), scale = FALSE)
|
|
|
2922 |
> Cscore <- scale(as.numeric(col(occupationalStatus)), scale = FALSE)
|
|
|
2923 |
> modUnif <- glm(Freq ~ origin + destination + Diag + Rscore:Cscore,
|
|
|
2924 |
+ family = poisson, data = occupationalStatus)
|
|
|
2925 |
>
|
|
|
2926 |
> summary(modUnif)
|
|
|
2927 |
|
|
|
2928 |
Call:
|
|
|
2929 |
glm(formula = Freq ~ origin + destination + Diag + Rscore:Cscore,
|
|
|
2930 |
family = poisson, data = occupationalStatus)
|
|
|
2931 |
|
|
|
2932 |
Coefficients:
|
|
|
2933 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
2934 |
(Intercept) 0.568592 0.183358 3.101 0.001929 **
|
|
|
2935 |
origin2 0.431314 0.149415 2.887 0.003893 **
|
|
|
2936 |
origin3 1.461862 0.131141 11.147 < 2e-16 ***
|
|
|
2937 |
origin4 1.788731 0.126588 14.130 < 2e-16 ***
|
|
|
2938 |
origin5 0.441011 0.144844 3.045 0.002329 **
|
|
|
2939 |
origin6 2.491735 0.121219 20.556 < 2e-16 ***
|
|
|
2940 |
origin7 1.127536 0.129032 8.738 < 2e-16 ***
|
|
|
2941 |
origin8 0.796445 0.131863 6.040 1.54e-09 ***
|
|
|
2942 |
destination2 0.873202 0.166902 5.232 1.68e-07 ***
|
|
|
2943 |
destination3 1.813992 0.153861 11.790 < 2e-16 ***
|
|
|
2944 |
destination4 2.082515 0.150887 13.802 < 2e-16 ***
|
|
|
2945 |
destination5 1.366383 0.155590 8.782 < 2e-16 ***
|
|
|
2946 |
destination6 2.816369 0.146054 19.283 < 2e-16 ***
|
|
|
2947 |
destination7 1.903918 0.147810 12.881 < 2e-16 ***
|
|
|
2948 |
destination8 1.398585 0.151658 9.222 < 2e-16 ***
|
|
|
2949 |
Diag1 1.665495 0.237383 7.016 2.28e-12 ***
|
|
|
2950 |
Diag2 0.959681 0.212122 4.524 6.06e-06 ***
|
|
|
2951 |
Diag3 0.021750 0.156530 0.139 0.889490
|
|
|
2952 |
Diag4 0.226399 0.124364 1.820 0.068689 .
|
|
|
2953 |
Diag5 0.808646 0.229754 3.520 0.000432 ***
|
|
|
2954 |
Diag6 0.132277 0.077191 1.714 0.086597 .
|
|
|
2955 |
Diag7 0.506709 0.115936 4.371 1.24e-05 ***
|
|
|
2956 |
Diag8 0.221880 0.134803 1.646 0.099771 .
|
|
|
2957 |
Rscore:Cscore 0.136974 0.007489 18.289 < 2e-16 ***
|
|
|
2958 |
---
|
| 61435 |
ripley |
2959 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2960 |
|
|
|
2961 |
(Dispersion parameter for poisson family taken to be 1)
|
|
|
2962 |
|
|
|
2963 |
Null deviance: 4679.004 on 63 degrees of freedom
|
|
|
2964 |
Residual deviance: 58.436 on 40 degrees of freedom
|
|
|
2965 |
AIC: 428.78
|
|
|
2966 |
|
|
|
2967 |
Number of Fisher Scoring iterations: 4
|
|
|
2968 |
|
|
|
2969 |
> plot(modUnif) # 4 plots, with warning about h_ii ~= 1
|
| 60841 |
ripley |
2970 |
Warning: not plotting observations with leverage one:
|
| 56186 |
murdoch |
2971 |
1, 10, 19, 28, 37, 46, 55, 64
|
|
|
2972 |
>
|
|
|
2973 |
>
|
|
|
2974 |
>
|
|
|
2975 |
> cleanEx()
|
| 87776 |
maechler |
2976 |
> nameEx("penguins")
|
|
|
2977 |
> ### * penguins
|
|
|
2978 |
>
|
|
|
2979 |
> flush(stderr()); flush(stdout())
|
|
|
2980 |
>
|
|
|
2981 |
> ### Name: penguins
|
|
|
2982 |
> ### Title: Measurements of Penguins near Palmer Station, Antarctica
|
|
|
2983 |
> ### Aliases: penguins penguins_raw
|
|
|
2984 |
> ### Keywords: datasets
|
|
|
2985 |
>
|
|
|
2986 |
> ### ** Examples
|
|
|
2987 |
>
|
|
|
2988 |
> ## view summaries
|
|
|
2989 |
> summary(penguins)
|
|
|
2990 |
species island bill_len bill_dep
|
|
|
2991 |
Adelie :152 Biscoe :168 Min. :32.10 Min. :13.10
|
|
|
2992 |
Chinstrap: 68 Dream :124 1st Qu.:39.23 1st Qu.:15.60
|
|
|
2993 |
Gentoo :124 Torgersen: 52 Median :44.45 Median :17.30
|
|
|
2994 |
Mean :43.92 Mean :17.15
|
|
|
2995 |
3rd Qu.:48.50 3rd Qu.:18.70
|
|
|
2996 |
Max. :59.60 Max. :21.50
|
| 89306 |
smeyer |
2997 |
NAs :2 NAs :2
|
| 87776 |
maechler |
2998 |
flipper_len body_mass sex year
|
|
|
2999 |
Min. :172.0 Min. :2700 female:165 Min. :2007
|
|
|
3000 |
1st Qu.:190.0 1st Qu.:3550 male :168 1st Qu.:2007
|
| 89306 |
smeyer |
3001 |
Median :197.0 Median :4050 NAs : 11 Median :2008
|
| 87776 |
maechler |
3002 |
Mean :200.9 Mean :4202 Mean :2008
|
|
|
3003 |
3rd Qu.:213.0 3rd Qu.:4750 3rd Qu.:2009
|
|
|
3004 |
Max. :231.0 Max. :6300 Max. :2009
|
| 89306 |
smeyer |
3005 |
NAs :2 NAs :2
|
| 89316 |
smeyer |
3006 |
> summary(penguins_raw) # categorical variables are (still) character vectors
|
|
|
3007 |
studyName Sample Number Species Region
|
|
|
3008 |
Length :344 Min. : 1.00 Length :344 Length :344
|
|
|
3009 |
N.unique : 3 1st Qu.: 29.00 N.unique : 3 N.unique : 1
|
|
|
3010 |
N.blank : 0 Median : 58.00 N.blank : 0 N.blank : 0
|
|
|
3011 |
Min.nchar: 7 Mean : 63.15 Min.nchar: 33 Min.nchar: 6
|
|
|
3012 |
Max.nchar: 7 3rd Qu.: 95.25 Max.nchar: 41 Max.nchar: 6
|
|
|
3013 |
Max. :152.00
|
|
|
3014 |
|
|
|
3015 |
Island Stage Individual ID Clutch Completion
|
|
|
3016 |
Length :344 Length :344 Length :344 Length :344
|
|
|
3017 |
N.unique : 3 N.unique : 1 N.unique :190 N.unique : 2
|
|
|
3018 |
N.blank : 0 N.blank : 0 N.blank : 0 N.blank : 0
|
|
|
3019 |
Min.nchar: 5 Min.nchar: 18 Min.nchar: 4 Min.nchar: 2
|
|
|
3020 |
Max.nchar: 9 Max.nchar: 18 Max.nchar: 6 Max.nchar: 3
|
|
|
3021 |
|
|
|
3022 |
|
| 87776 |
maechler |
3023 |
Date Egg Culmen Length (mm) Culmen Depth (mm) Flipper Length (mm)
|
|
|
3024 |
Min. :2007-11-09 Min. :32.10 Min. :13.10 Min. :172.0
|
|
|
3025 |
1st Qu.:2007-11-28 1st Qu.:39.23 1st Qu.:15.60 1st Qu.:190.0
|
|
|
3026 |
Median :2008-11-09 Median :44.45 Median :17.30 Median :197.0
|
|
|
3027 |
Mean :2008-11-27 Mean :43.92 Mean :17.15 Mean :200.9
|
|
|
3028 |
3rd Qu.:2009-11-16 3rd Qu.:48.50 3rd Qu.:18.70 3rd Qu.:213.0
|
|
|
3029 |
Max. :2009-12-01 Max. :59.60 Max. :21.50 Max. :231.0
|
| 89306 |
smeyer |
3030 |
NAs :2 NAs :2 NAs :2
|
| 89316 |
smeyer |
3031 |
Body Mass (g) Sex Delta 15 N (o/oo) Delta 13 C (o/oo)
|
|
|
3032 |
Min. :2700 Length :344 Min. : 7.632 Min. :-27.02
|
|
|
3033 |
1st Qu.:3550 N.unique : 2 1st Qu.: 8.300 1st Qu.:-26.32
|
|
|
3034 |
Median :4050 N.blank : 0 Median : 8.652 Median :-25.83
|
|
|
3035 |
Mean :4202 Min.nchar: 4 Mean : 8.733 Mean :-25.69
|
|
|
3036 |
3rd Qu.:4750 Max.nchar: 6 3rd Qu.: 9.172 3rd Qu.:-25.06
|
|
|
3037 |
Max. :6300 NAs : 11 Max. :10.025 Max. :-23.79
|
|
|
3038 |
NAs :2 NAs :14 NAs :13
|
|
|
3039 |
Comments
|
|
|
3040 |
Length :344
|
|
|
3041 |
N.unique : 10
|
|
|
3042 |
N.blank : 0
|
|
|
3043 |
Min.nchar: 18
|
|
|
3044 |
Max.nchar: 68
|
|
|
3045 |
NAs :290
|
|
|
3046 |
|
|
|
3047 |
> ## summarize character vectors as factors:
|
|
|
3048 |
> summary(penguins_raw, character.method = "factor")
|
| 87776 |
maechler |
3049 |
studyName Sample Number Species
|
|
|
3050 |
PAL0708:110 Min. : 1.00 Adelie Penguin (Pygoscelis adeliae) :152
|
|
|
3051 |
PAL0809:114 1st Qu.: 29.00 Chinstrap penguin (Pygoscelis antarctica): 68
|
|
|
3052 |
PAL0910:120 Median : 58.00 Gentoo penguin (Pygoscelis papua) :124
|
|
|
3053 |
Mean : 63.15
|
|
|
3054 |
3rd Qu.: 95.25
|
|
|
3055 |
Max. :152.00
|
|
|
3056 |
|
|
|
3057 |
Region Island Stage Individual ID
|
|
|
3058 |
Anvers:344 Biscoe :168 Adult, 1 Egg Stage:344 N13A1 : 3
|
|
|
3059 |
Dream :124 N13A2 : 3
|
|
|
3060 |
Torgersen: 52 N18A1 : 3
|
|
|
3061 |
N18A2 : 3
|
|
|
3062 |
N21A1 : 3
|
|
|
3063 |
N21A2 : 3
|
|
|
3064 |
(Other):326
|
|
|
3065 |
Clutch Completion Date Egg Culmen Length (mm) Culmen Depth (mm)
|
|
|
3066 |
No : 36 Min. :2007-11-09 Min. :32.10 Min. :13.10
|
|
|
3067 |
Yes:308 1st Qu.:2007-11-28 1st Qu.:39.23 1st Qu.:15.60
|
|
|
3068 |
Median :2008-11-09 Median :44.45 Median :17.30
|
|
|
3069 |
Mean :2008-11-27 Mean :43.92 Mean :17.15
|
|
|
3070 |
3rd Qu.:2009-11-16 3rd Qu.:48.50 3rd Qu.:18.70
|
|
|
3071 |
Max. :2009-12-01 Max. :59.60 Max. :21.50
|
| 89306 |
smeyer |
3072 |
NAs :2 NAs :2
|
| 87776 |
maechler |
3073 |
Flipper Length (mm) Body Mass (g) Sex Delta 15 N (o/oo)
|
|
|
3074 |
Min. :172.0 Min. :2700 FEMALE:165 Min. : 7.632
|
|
|
3075 |
1st Qu.:190.0 1st Qu.:3550 MALE :168 1st Qu.: 8.300
|
| 89306 |
smeyer |
3076 |
Median :197.0 Median :4050 NAs : 11 Median : 8.652
|
| 87776 |
maechler |
3077 |
Mean :200.9 Mean :4202 Mean : 8.733
|
|
|
3078 |
3rd Qu.:213.0 3rd Qu.:4750 3rd Qu.: 9.172
|
|
|
3079 |
Max. :231.0 Max. :6300 Max. :10.025
|
| 89306 |
smeyer |
3080 |
NAs :2 NAs :2 NAs :14
|
| 87776 |
maechler |
3081 |
Delta 13 C (o/oo) Comments
|
|
|
3082 |
Min. :-27.02 Nest never observed with full clutch.: 34
|
|
|
3083 |
1st Qu.:-26.32 Not enough blood for isotopes. : 7
|
|
|
3084 |
Median :-25.83 Sexing primers did not amplify. : 4
|
|
|
3085 |
Mean :-25.69 No blood sample obtained for sexing. : 2
|
|
|
3086 |
3rd Qu.:-25.06 No blood sample obtained. : 2
|
|
|
3087 |
Max. :-23.79 (Other) : 5
|
| 89306 |
smeyer |
3088 |
NAs :13 NAs :290
|
| 87776 |
maechler |
3089 |
>
|
|
|
3090 |
> ## visualise distribution across factors
|
|
|
3091 |
> plot(island ~ species, data = penguins)
|
|
|
3092 |
> plot(sex ~ interaction(island, species, sep = "\n"), data = penguins)
|
|
|
3093 |
>
|
|
|
3094 |
> ## bill depth vs. length by species (color) and sex (symbol):
|
|
|
3095 |
> ## positive correlations for all species, males tend to have bigger bills
|
|
|
3096 |
> sym <- c(1, 16)
|
|
|
3097 |
> pal <- c("darkorange","purple","cyan4")
|
|
|
3098 |
> plot(bill_dep ~ bill_len, data = penguins, pch = sym[sex], col = pal[species])
|
|
|
3099 |
>
|
|
|
3100 |
> ## simplified sex dimorphism analysis for Adelie species:
|
|
|
3101 |
> ## proportion of males increases with several size measurements
|
|
|
3102 |
> adelie <- subset(penguins, species == "Adelie")
|
|
|
3103 |
> plot(sex ~ bill_len, data = adelie)
|
|
|
3104 |
> plot(sex ~ bill_dep, data = adelie)
|
|
|
3105 |
> plot(sex ~ body_mass, data = adelie)
|
|
|
3106 |
> m <- glm(sex ~ bill_len + bill_dep + body_mass, data = adelie, family = binomial)
|
|
|
3107 |
> summary(m)
|
|
|
3108 |
|
|
|
3109 |
Call:
|
|
|
3110 |
glm(formula = sex ~ bill_len + bill_dep + body_mass, family = binomial,
|
|
|
3111 |
data = adelie)
|
|
|
3112 |
|
|
|
3113 |
Coefficients:
|
|
|
3114 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
3115 |
(Intercept) -85.088438 17.811711 -4.777 1.78e-06 ***
|
|
|
3116 |
bill_len 0.840401 0.237656 3.536 0.000406 ***
|
|
|
3117 |
bill_dep 1.305989 0.423911 3.081 0.002064 **
|
|
|
3118 |
body_mass 0.007790 0.001969 3.957 7.59e-05 ***
|
|
|
3119 |
---
|
|
|
3120 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
|
|
3121 |
|
|
|
3122 |
(Dispersion parameter for binomial family taken to be 1)
|
|
|
3123 |
|
|
|
3124 |
Null deviance: 202.399 on 145 degrees of freedom
|
|
|
3125 |
Residual deviance: 54.024 on 142 degrees of freedom
|
|
|
3126 |
(6 observations deleted due to missingness)
|
|
|
3127 |
AIC: 62.024
|
|
|
3128 |
|
|
|
3129 |
Number of Fisher Scoring iterations: 8
|
|
|
3130 |
|
|
|
3131 |
>
|
|
|
3132 |
> ## Produce the long variable names as from {palmerpenguins} pkg:
|
|
|
3133 |
> long_nms <- sub("len", "length_mm",
|
|
|
3134 |
+ sub("dep","depth_mm",
|
|
|
3135 |
+ sub("mass", "mass_g", colnames(penguins))))
|
|
|
3136 |
> ## compare long and short names:
|
|
|
3137 |
> noquote(rbind(long_nms, nms = colnames(penguins)))
|
|
|
3138 |
[,1] [,2] [,3] [,4] [,5]
|
|
|
3139 |
long_nms species island bill_length_mm bill_depth_mm flipper_length_mm
|
|
|
3140 |
nms species island bill_len bill_dep flipper_len
|
|
|
3141 |
[,6] [,7] [,8]
|
|
|
3142 |
long_nms body_mass_g sex year
|
|
|
3143 |
nms body_mass sex year
|
|
|
3144 |
>
|
|
|
3145 |
> ## Not run:
|
|
|
3146 |
> ##D # << keeping shorter 'penguins' names in this example:
|
|
|
3147 |
> ##D colnames(penguins) <- long_nms
|
|
|
3148 |
> ## End(Not run)
|
|
|
3149 |
>
|
|
|
3150 |
>
|
|
|
3151 |
>
|
|
|
3152 |
> cleanEx()
|
| 56186 |
murdoch |
3153 |
> nameEx("precip")
|
|
|
3154 |
> ### * precip
|
|
|
3155 |
>
|
|
|
3156 |
> flush(stderr()); flush(stdout())
|
|
|
3157 |
>
|
|
|
3158 |
> ### Name: precip
|
| 88211 |
maechler |
3159 |
> ### Title: Annual Precipitation in Selected US Cities
|
| 56186 |
murdoch |
3160 |
> ### Aliases: precip
|
|
|
3161 |
> ### Keywords: datasets
|
|
|
3162 |
>
|
|
|
3163 |
> ### ** Examples
|
|
|
3164 |
>
|
| 88388 |
smeyer |
3165 |
> require(stats); require(graphics)
|
| 56186 |
murdoch |
3166 |
> dotchart(precip[order(precip)], main = "precip data")
|
|
|
3167 |
> title(sub = "Average annual precipitation (in.)")
|
|
|
3168 |
>
|
| 88388 |
smeyer |
3169 |
> ## Recreate pre-4.6.0 version of the dataset with typos in 2 names
|
|
|
3170 |
> precip.2 <- setNames(precip, local({
|
|
|
3171 |
+ nms <- names(precip)
|
|
|
3172 |
+ nms[nms == "Pittsburgh"] <- "Pittsburg"
|
|
|
3173 |
+ nms[nms == "Bismarck"] <- "Bismark"
|
|
|
3174 |
+ nms
|
|
|
3175 |
+ }))
|
|
|
3176 |
> ## Recreate original, pre-3.3.3 version of the dataset with one more typo
|
|
|
3177 |
> precip.O <- setNames(precip.2,
|
|
|
3178 |
+ replace(names(precip.2), names(precip.2) == "Cincinnati", "Cincinati"))
|
| 56186 |
murdoch |
3179 |
>
|
| 88388 |
smeyer |
3180 |
> stopifnot(
|
|
|
3181 |
+ all(precip == precip.2), all(precip == precip.O), # just name changes
|
|
|
3182 |
+ setequal(match(c("Bismarck", "Cincinnati", "Pittsburgh"), names(precip)),
|
|
|
3183 |
+ c(45, 46, 52)),
|
|
|
3184 |
+ identical(names(precip)[-c(45, 46, 52)], names(precip.O)[-c(45, 46, 52)])
|
|
|
3185 |
+ )
|
| 56186 |
murdoch |
3186 |
>
|
| 71678 |
ripley |
3187 |
>
|
| 88388 |
smeyer |
3188 |
>
|
| 56186 |
murdoch |
3189 |
> cleanEx()
|
|
|
3190 |
> nameEx("presidents")
|
|
|
3191 |
> ### * presidents
|
|
|
3192 |
>
|
|
|
3193 |
> flush(stderr()); flush(stdout())
|
|
|
3194 |
>
|
|
|
3195 |
> ### Name: presidents
|
|
|
3196 |
> ### Title: Quarterly Approval Ratings of US Presidents
|
|
|
3197 |
> ### Aliases: presidents
|
|
|
3198 |
> ### Keywords: datasets
|
|
|
3199 |
>
|
|
|
3200 |
> ### ** Examples
|
|
|
3201 |
>
|
|
|
3202 |
> require(stats); require(graphics)
|
|
|
3203 |
> plot(presidents, las = 1, ylab = "Approval rating (%)",
|
|
|
3204 |
+ main = "presidents data")
|
|
|
3205 |
>
|
|
|
3206 |
>
|
|
|
3207 |
>
|
|
|
3208 |
> cleanEx()
|
|
|
3209 |
> nameEx("pressure")
|
|
|
3210 |
> ### * pressure
|
|
|
3211 |
>
|
|
|
3212 |
> flush(stderr()); flush(stdout())
|
|
|
3213 |
>
|
|
|
3214 |
> ### Name: pressure
|
|
|
3215 |
> ### Title: Vapor Pressure of Mercury as a Function of Temperature
|
|
|
3216 |
> ### Aliases: pressure
|
|
|
3217 |
> ### Keywords: datasets
|
|
|
3218 |
>
|
|
|
3219 |
> ### ** Examples
|
|
|
3220 |
>
|
|
|
3221 |
> require(graphics)
|
|
|
3222 |
> plot(pressure, xlab = "Temperature (deg C)",
|
|
|
3223 |
+ ylab = "Pressure (mm of Hg)",
|
|
|
3224 |
+ main = "pressure data: Vapor Pressure of Mercury")
|
|
|
3225 |
> plot(pressure, xlab = "Temperature (deg C)", log = "y",
|
|
|
3226 |
+ ylab = "Pressure (mm of Hg)",
|
|
|
3227 |
+ main = "pressure data: Vapor Pressure of Mercury")
|
|
|
3228 |
>
|
|
|
3229 |
>
|
|
|
3230 |
>
|
|
|
3231 |
> cleanEx()
|
|
|
3232 |
> nameEx("quakes")
|
|
|
3233 |
> ### * quakes
|
|
|
3234 |
>
|
|
|
3235 |
> flush(stderr()); flush(stdout())
|
|
|
3236 |
>
|
|
|
3237 |
> ### Name: quakes
|
|
|
3238 |
> ### Title: Locations of Earthquakes off Fiji
|
|
|
3239 |
> ### Aliases: quakes
|
|
|
3240 |
> ### Keywords: datasets
|
|
|
3241 |
>
|
|
|
3242 |
> ### ** Examples
|
|
|
3243 |
>
|
|
|
3244 |
> require(graphics)
|
| 61157 |
ripley |
3245 |
> pairs(quakes, main = "Fiji Earthquakes, N = 1000", cex.main = 1.2, pch = ".")
|
| 56186 |
murdoch |
3246 |
>
|
|
|
3247 |
>
|
|
|
3248 |
>
|
|
|
3249 |
> cleanEx()
|
|
|
3250 |
> nameEx("randu")
|
|
|
3251 |
> ### * randu
|
|
|
3252 |
>
|
|
|
3253 |
> flush(stderr()); flush(stdout())
|
|
|
3254 |
>
|
|
|
3255 |
> ### Name: randu
|
|
|
3256 |
> ### Title: Random Numbers from Congruential Generator RANDU
|
|
|
3257 |
> ### Aliases: randu
|
|
|
3258 |
> ### Keywords: datasets
|
|
|
3259 |
>
|
|
|
3260 |
> ### ** Examples
|
|
|
3261 |
>
|
| 86023 |
smeyer |
3262 |
> ## We could re-generate the dataset by the following R code
|
|
|
3263 |
> seed <- as.double(1)
|
|
|
3264 |
> RANDU <- function() {
|
|
|
3265 |
+ seed <<- ((2^16 + 3) * seed) %% (2^31)
|
|
|
3266 |
+ seed/(2^31)
|
|
|
3267 |
+ }
|
|
|
3268 |
> myrandu <- matrix(NA_real_, 400, 3, dimnames = list(NULL, c("x","y","z")))
|
|
|
3269 |
> for(i in 1:400) {
|
|
|
3270 |
+ U <- c(RANDU(), RANDU(), RANDU(), RANDU(), RANDU())
|
|
|
3271 |
+ myrandu[i,] <- round(U[1:3], 6)
|
|
|
3272 |
+ }
|
|
|
3273 |
> stopifnot(all.equal(randu, as.data.frame(myrandu), tolerance = 1e-5))
|
| 56186 |
murdoch |
3274 |
>
|
|
|
3275 |
>
|
|
|
3276 |
>
|
|
|
3277 |
> cleanEx()
|
|
|
3278 |
> nameEx("sleep")
|
|
|
3279 |
> ### * sleep
|
|
|
3280 |
>
|
|
|
3281 |
> flush(stderr()); flush(stdout())
|
|
|
3282 |
>
|
|
|
3283 |
> ### Name: sleep
|
|
|
3284 |
> ### Title: Student's Sleep Data
|
|
|
3285 |
> ### Aliases: sleep
|
|
|
3286 |
> ### Keywords: datasets
|
|
|
3287 |
>
|
|
|
3288 |
> ### ** Examples
|
|
|
3289 |
>
|
|
|
3290 |
> require(stats)
|
|
|
3291 |
> ## Student's paired t-test
|
|
|
3292 |
> with(sleep,
|
|
|
3293 |
+ t.test(extra[group == 1],
|
|
|
3294 |
+ extra[group == 2], paired = TRUE))
|
|
|
3295 |
|
|
|
3296 |
Paired t-test
|
|
|
3297 |
|
| 61430 |
ripley |
3298 |
data: extra[group == 1] and extra[group == 2]
|
| 56186 |
murdoch |
3299 |
t = -4.0621, df = 9, p-value = 0.002833
|
| 80819 |
hornik |
3300 |
alternative hypothesis: true mean difference is not equal to 0
|
| 56186 |
murdoch |
3301 |
95 percent confidence interval:
|
| 61430 |
ripley |
3302 |
-2.4598858 -0.7001142
|
| 56186 |
murdoch |
3303 |
sample estimates:
|
| 80819 |
hornik |
3304 |
mean difference
|
|
|
3305 |
-1.58
|
| 56186 |
murdoch |
3306 |
|
|
|
3307 |
>
|
| 59811 |
maechler |
3308 |
> ## The sleep *prolongations*
|
|
|
3309 |
> sleep1 <- with(sleep, extra[group == 2] - extra[group == 1])
|
|
|
3310 |
> summary(sleep1)
|
|
|
3311 |
Min. 1st Qu. Median Mean 3rd Qu. Max.
|
|
|
3312 |
0.00 1.05 1.30 1.58 1.70 4.60
|
|
|
3313 |
> stripchart(sleep1, method = "stack", xlab = "hours",
|
|
|
3314 |
+ main = "Sleep prolongation (n = 10)")
|
| 61157 |
ripley |
3315 |
> boxplot(sleep1, horizontal = TRUE, add = TRUE,
|
|
|
3316 |
+ at = .6, pars = list(boxwex = 0.5, staplewex = 0.25))
|
| 56186 |
murdoch |
3317 |
>
|
|
|
3318 |
>
|
| 59811 |
maechler |
3319 |
>
|
| 56186 |
murdoch |
3320 |
> cleanEx()
|
|
|
3321 |
> nameEx("stackloss")
|
|
|
3322 |
> ### * stackloss
|
|
|
3323 |
>
|
|
|
3324 |
> flush(stderr()); flush(stdout())
|
|
|
3325 |
>
|
|
|
3326 |
> ### Name: stackloss
|
|
|
3327 |
> ### Title: Brownlee's Stack Loss Plant Data
|
|
|
3328 |
> ### Aliases: stackloss stack.loss stack.x
|
|
|
3329 |
> ### Keywords: datasets
|
|
|
3330 |
>
|
|
|
3331 |
> ### ** Examples
|
|
|
3332 |
>
|
|
|
3333 |
> require(stats)
|
|
|
3334 |
> summary(lm.stack <- lm(stack.loss ~ stack.x))
|
|
|
3335 |
|
|
|
3336 |
Call:
|
|
|
3337 |
lm(formula = stack.loss ~ stack.x)
|
|
|
3338 |
|
|
|
3339 |
Residuals:
|
|
|
3340 |
Min 1Q Median 3Q Max
|
|
|
3341 |
-7.2377 -1.7117 -0.4551 2.3614 5.6978
|
|
|
3342 |
|
|
|
3343 |
Coefficients:
|
|
|
3344 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3345 |
(Intercept) -39.9197 11.8960 -3.356 0.00375 **
|
|
|
3346 |
stack.xAir.Flow 0.7156 0.1349 5.307 5.8e-05 ***
|
|
|
3347 |
stack.xWater.Temp 1.2953 0.3680 3.520 0.00263 **
|
|
|
3348 |
stack.xAcid.Conc. -0.1521 0.1563 -0.973 0.34405
|
|
|
3349 |
---
|
| 61435 |
ripley |
3350 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3351 |
|
|
|
3352 |
Residual standard error: 3.243 on 17 degrees of freedom
|
| 61435 |
ripley |
3353 |
Multiple R-squared: 0.9136, Adjusted R-squared: 0.8983
|
|
|
3354 |
F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09
|
| 56186 |
murdoch |
3355 |
|
|
|
3356 |
>
|
|
|
3357 |
>
|
|
|
3358 |
>
|
|
|
3359 |
> cleanEx()
|
| 82568 |
maechler |
3360 |
> nameEx("state")
|
|
|
3361 |
> ### * state
|
|
|
3362 |
>
|
|
|
3363 |
> flush(stderr()); flush(stdout())
|
|
|
3364 |
>
|
|
|
3365 |
> ### Name: state
|
|
|
3366 |
> ### Title: US State Facts and Figures
|
|
|
3367 |
> ### Aliases: state state.abb state.area state.center state.division
|
|
|
3368 |
> ### state.name state.region state.x77
|
|
|
3369 |
> ### Keywords: datasets
|
|
|
3370 |
>
|
|
|
3371 |
> ### ** Examples
|
|
|
3372 |
>
|
|
|
3373 |
> (dst <- dxy <- data.frame(state.center, row.names=state.abb))
|
|
|
3374 |
x y
|
|
|
3375 |
AL -86.7509 32.5901
|
|
|
3376 |
AK -127.2500 49.2500
|
|
|
3377 |
AZ -111.6250 34.2192
|
|
|
3378 |
AR -92.2992 34.7336
|
|
|
3379 |
CA -119.7730 36.5341
|
|
|
3380 |
CO -105.5130 38.6777
|
|
|
3381 |
CT -72.3573 41.5928
|
|
|
3382 |
DE -74.9841 38.6777
|
|
|
3383 |
FL -81.6850 27.8744
|
|
|
3384 |
GA -83.3736 32.3329
|
|
|
3385 |
HI -126.2500 31.7500
|
|
|
3386 |
ID -113.9300 43.5648
|
|
|
3387 |
IL -89.3776 40.0495
|
|
|
3388 |
IN -86.0808 40.0495
|
|
|
3389 |
IA -93.3714 41.9358
|
|
|
3390 |
KS -98.1156 38.4204
|
|
|
3391 |
KY -84.7674 37.3915
|
|
|
3392 |
LA -92.2724 30.6181
|
|
|
3393 |
ME -68.9801 45.6226
|
|
|
3394 |
MD -76.6459 39.2778
|
|
|
3395 |
MA -71.5800 42.3645
|
|
|
3396 |
MI -84.6870 43.1361
|
|
|
3397 |
MN -94.6043 46.3943
|
|
|
3398 |
MS -89.8065 32.6758
|
|
|
3399 |
MO -92.5137 38.3347
|
|
|
3400 |
MT -109.3200 46.8230
|
|
|
3401 |
NE -99.5898 41.3356
|
|
|
3402 |
NV -116.8510 39.1063
|
|
|
3403 |
NH -71.3924 43.3934
|
|
|
3404 |
NJ -74.2336 39.9637
|
|
|
3405 |
NM -105.9420 34.4764
|
|
|
3406 |
NY -75.1449 43.1361
|
|
|
3407 |
NC -78.4686 35.4195
|
|
|
3408 |
ND -100.0990 47.2517
|
|
|
3409 |
OH -82.5963 40.2210
|
|
|
3410 |
OK -97.1239 35.5053
|
|
|
3411 |
OR -120.0680 43.9078
|
|
|
3412 |
PA -77.4500 40.9069
|
|
|
3413 |
RI -71.1244 41.5928
|
|
|
3414 |
SC -80.5056 33.6190
|
|
|
3415 |
SD -99.7238 44.3365
|
|
|
3416 |
TN -86.4560 35.6767
|
|
|
3417 |
TX -98.7857 31.3897
|
|
|
3418 |
UT -111.3300 39.1063
|
|
|
3419 |
VT -72.5450 44.2508
|
|
|
3420 |
VA -78.2005 37.5630
|
|
|
3421 |
WA -119.7460 47.4231
|
|
|
3422 |
WV -80.6665 38.4204
|
|
|
3423 |
WI -89.9941 44.5937
|
|
|
3424 |
WY -107.2560 43.0504
|
|
|
3425 |
> ## Alaska and Hawaii are placed just off the West Coast (for compact map drawing):
|
|
|
3426 |
> dst[c("AK", "HI"),]
|
|
|
3427 |
x y
|
|
|
3428 |
AK -127.25 49.25
|
|
|
3429 |
HI -126.25 31.75
|
|
|
3430 |
> ## state.center2 := version of state.center with "correct" coordinates for AK & HI:
|
|
|
3431 |
> ## From https://pubs.usgs.gov/gip/Elevations-Distances/elvadist.html#Geographic%20Centers
|
|
|
3432 |
> ## Alaska 63°50' N., 152°00' W., 60 miles northwest of Mount McKinley
|
|
|
3433 |
> ## Hawaii 20°15' N., 156°20' W., off Maui Island
|
|
|
3434 |
> dxy["AK",] <- c(-152. , 63.83) # or c(-152.11, 65.17)
|
|
|
3435 |
> dxy["HI",] <- c(-156.33, 20.25) # or c(-156.69, 20.89)
|
|
|
3436 |
> state.center2 <- as.list(dxy)
|
|
|
3437 |
>
|
|
|
3438 |
> plot(dxy, asp=1.2, pch=3, col=2)
|
|
|
3439 |
> text(state.center2, state.abb, cex=1/2, pos=4, offset=1/4)
|
|
|
3440 |
> i <- c("AK","HI")
|
|
|
3441 |
> do.call(arrows, c(setNames(c(dst[i,], dxy[i,]), c("x0","y0", "x1","y1")),
|
|
|
3442 |
+ col=adjustcolor(4, .7), length=1/8))
|
|
|
3443 |
> points(dst[i,], col=2)
|
|
|
3444 |
> if(FALSE) { # if(require("maps")) {
|
|
|
3445 |
+ map("state", interior = FALSE, add = TRUE)
|
|
|
3446 |
+ map("state", boundary = FALSE, lty = 2, add = TRUE)
|
|
|
3447 |
+ }
|
|
|
3448 |
>
|
|
|
3449 |
>
|
|
|
3450 |
>
|
|
|
3451 |
> cleanEx()
|
| 56186 |
murdoch |
3452 |
> nameEx("sunspot.month")
|
|
|
3453 |
> ### * sunspot.month
|
|
|
3454 |
>
|
|
|
3455 |
> flush(stderr()); flush(stdout())
|
|
|
3456 |
>
|
|
|
3457 |
> ### Name: sunspot.month
|
| 65249 |
maechler |
3458 |
> ### Title: Monthly Sunspot Data, from 1749 to "Present"
|
| 87340 |
maechler |
3459 |
> ### Aliases: sunspot.m2014 sunspot.month
|
| 56186 |
murdoch |
3460 |
> ### Keywords: datasets
|
|
|
3461 |
>
|
|
|
3462 |
> ### ** Examples
|
|
|
3463 |
>
|
|
|
3464 |
> require(stats); require(graphics)
|
| 61435 |
ripley |
3465 |
> ## Compare the monthly series
|
| 65249 |
maechler |
3466 |
> plot (sunspot.month,
|
| 87311 |
ripley |
3467 |
+ main="sunspot.month & sunspots [package 'datasets']", col=2)
|
| 87340 |
maechler |
3468 |
> lines(sunspots) # -> clear recalibration (to *larger* values)
|
|
|
3469 |
> at. <- seq(1750, 2030, by=10)
|
|
|
3470 |
> atyr <- at.[at. %% 50 != 0] #%
|
|
|
3471 |
> Axis(time(sunspot.month), at = atyr, side = 1,
|
|
|
3472 |
+ tck= -1/100, padj = -1.5, cex.axis = 3/4)
|
| 56186 |
murdoch |
3473 |
>
|
|
|
3474 |
> ## Now look at the difference :
|
|
|
3475 |
> all(tsp(sunspots) [c(1,3)] ==
|
|
|
3476 |
+ tsp(sunspot.month)[c(1,3)]) ## Start & Periodicity are the same
|
|
|
3477 |
[1] TRUE
|
|
|
3478 |
> n1 <- length(sunspots)
|
| 87340 |
maechler |
3479 |
> table(eq <- sunspots == sunspot.m2014[1:n1]) #> 143 are different !
|
| 56186 |
murdoch |
3480 |
|
|
|
3481 |
FALSE TRUE
|
| 65249 |
maechler |
3482 |
143 2677
|
| 61435 |
ripley |
3483 |
> i <- which(!eq)
|
| 56186 |
murdoch |
3484 |
> rug(time(eq)[i])
|
| 87340 |
maechler |
3485 |
> s1 <- sunspots[i] ; s2 <- sunspot.m2014[i]
|
| 65249 |
maechler |
3486 |
> cbind(i = i, time = time(sunspots)[i], sunspots = s1, ss.month = s2,
|
| 56186 |
murdoch |
3487 |
+ perc.diff = round(100*2*abs(s1-s2)/(s1+s2), 1))
|
| 65249 |
maechler |
3488 |
i time sunspots ss.month perc.diff
|
|
|
3489 |
[1,] 55 1753.500 22.2 22.0 0.9
|
|
|
3490 |
[2,] 838 1818.750 31.7 31.6 0.3
|
|
|
3491 |
[3,] 841 1819.000 32.5 32.8 0.9
|
|
|
3492 |
[4,] 862 1820.750 9.0 8.9 1.1
|
|
|
3493 |
[5,] 864 1820.917 9.7 9.1 6.4
|
|
|
3494 |
[6,] 866 1821.083 4.3 4.2 2.4
|
|
|
3495 |
[7,] 876 1821.917 0.0 0.2 200.0
|
|
|
3496 |
[8,] 901 1824.000 21.6 21.7 0.5
|
|
|
3497 |
[9,] 917 1825.333 15.4 15.5 0.6
|
|
|
3498 |
[10,] 920 1825.583 25.4 25.7 1.2
|
|
|
3499 |
[11,] 943 1827.500 42.9 42.3 1.4
|
|
|
3500 |
[12,] 946 1827.750 57.2 56.1 1.9
|
|
|
3501 |
[13,] 955 1828.500 54.3 54.2 0.2
|
|
|
3502 |
[14,] 960 1828.917 46.6 46.9 0.6
|
|
|
3503 |
[15,] 965 1829.333 67.5 67.4 0.1
|
|
|
3504 |
[16,] 968 1829.583 78.3 77.6 0.9
|
|
|
3505 |
[17,] 976 1830.250 107.1 106.3 0.7
|
|
|
3506 |
[18,] 988 1831.250 54.6 54.5 0.2
|
|
|
3507 |
[19,] 992 1831.583 54.9 55.0 0.2
|
|
|
3508 |
[20,] 994 1831.750 46.2 46.3 0.2
|
|
|
3509 |
[21,] 998 1832.083 55.5 55.6 0.2
|
|
|
3510 |
[22,] 1003 1832.500 13.9 14.0 0.7
|
|
|
3511 |
[23,] 1047 1836.167 98.1 98.2 0.1
|
|
|
3512 |
[24,] 1061 1837.333 111.3 111.7 0.4
|
|
|
3513 |
[25,] 1081 1839.000 107.6 105.6 1.9
|
|
|
3514 |
[26,] 1087 1839.500 84.7 84.8 0.1
|
|
|
3515 |
[27,] 1090 1839.750 90.8 90.9 0.1
|
|
|
3516 |
[28,] 1092 1839.917 63.6 63.7 0.2
|
|
|
3517 |
[29,] 1095 1840.167 55.5 67.8 20.0
|
|
|
3518 |
[30,] 1102 1840.750 49.8 55.0 9.9
|
|
|
3519 |
[31,] 1105 1841.000 24.0 24.1 0.4
|
|
|
3520 |
[32,] 1108 1841.250 42.6 40.2 5.8
|
|
|
3521 |
[33,] 1109 1841.333 67.4 67.5 0.1
|
|
|
3522 |
[34,] 1113 1841.667 35.1 36.5 3.9
|
|
|
3523 |
[35,] 1124 1842.583 26.5 26.6 0.4
|
|
|
3524 |
[36,] 1125 1842.667 18.5 18.4 0.5
|
|
|
3525 |
[37,] 1132 1843.250 8.8 9.5 7.7
|
|
|
3526 |
[38,] 1145 1844.333 12.0 11.6 3.4
|
|
|
3527 |
[39,] 1149 1844.667 6.9 7.0 1.4
|
|
|
3528 |
[40,] 1156 1845.250 56.9 57.0 0.2
|
|
|
3529 |
[41,] 1168 1846.250 69.2 69.3 0.1
|
|
|
3530 |
[42,] 1185 1847.667 161.2 160.9 0.2
|
|
|
3531 |
[43,] 1191 1848.167 108.9 108.6 0.3
|
|
|
3532 |
[44,] 1194 1848.417 123.8 129.0 4.1
|
|
|
3533 |
[45,] 1196 1848.583 132.5 132.6 0.1
|
|
|
3534 |
[46,] 1200 1848.917 159.9 159.5 0.3
|
|
|
3535 |
[47,] 1201 1849.000 156.7 157.0 0.2
|
|
|
3536 |
[48,] 1202 1849.083 131.7 131.8 0.1
|
|
|
3537 |
[49,] 1203 1849.167 96.5 96.2 0.3
|
|
|
3538 |
[50,] 1206 1849.417 81.2 81.1 0.1
|
|
|
3539 |
[51,] 1208 1849.583 61.3 67.7 9.9
|
|
|
3540 |
[52,] 1211 1849.833 99.7 99.0 0.7
|
|
|
3541 |
[53,] 1224 1850.917 60.0 61.0 1.7
|
|
|
3542 |
[54,] 1235 1851.833 50.9 51.0 0.2
|
|
|
3543 |
[55,] 1238 1852.083 67.5 66.4 1.6
|
|
|
3544 |
[56,] 1243 1852.500 42.0 42.1 0.2
|
|
|
3545 |
[57,] 1256 1853.583 50.4 50.5 0.2
|
|
|
3546 |
[58,] 1258 1853.750 42.3 42.4 0.2
|
|
|
3547 |
[59,] 1264 1854.250 26.4 26.5 0.4
|
|
|
3548 |
[60,] 1270 1854.750 12.7 12.6 0.8
|
|
|
3549 |
[61,] 1272 1854.917 21.4 21.6 0.9
|
|
|
3550 |
[62,] 1282 1855.750 9.7 9.6 1.0
|
|
|
3551 |
[63,] 1283 1855.833 4.3 4.2 2.4
|
|
|
3552 |
[64,] 1290 1856.417 5.0 5.2 3.9
|
|
|
3553 |
[65,] 1301 1857.333 29.2 28.5 2.4
|
|
|
3554 |
[66,] 1333 1860.000 81.5 82.4 1.1
|
|
|
3555 |
[67,] 1334 1860.083 88.0 88.3 0.3
|
|
|
3556 |
[68,] 1346 1861.083 77.8 77.7 0.1
|
|
|
3557 |
[69,] 1350 1861.417 87.8 88.1 0.3
|
|
|
3558 |
[70,] 1366 1862.750 42.0 41.9 0.2
|
|
|
3559 |
[71,] 1407 1866.167 24.6 24.5 0.4
|
|
|
3560 |
[72,] 1424 1867.583 4.9 4.8 2.1
|
|
|
3561 |
[73,] 1427 1867.833 9.3 9.6 3.2
|
|
|
3562 |
[74,] 1429 1868.000 15.6 15.5 0.6
|
|
|
3563 |
[75,] 1430 1868.083 15.8 15.7 0.6
|
|
|
3564 |
[76,] 1435 1868.500 28.6 29.0 1.4
|
|
|
3565 |
[77,] 1437 1868.667 43.8 47.2 7.5
|
|
|
3566 |
[78,] 1438 1868.750 61.7 61.6 0.2
|
|
|
3567 |
[79,] 1442 1869.083 59.3 59.9 1.0
|
|
|
3568 |
[80,] 1445 1869.333 104.0 103.9 0.1
|
|
|
3569 |
[81,] 1450 1869.750 59.4 59.3 0.2
|
|
|
3570 |
[82,] 1451 1869.833 77.4 78.1 0.9
|
|
|
3571 |
[83,] 1452 1869.917 104.3 104.4 0.1
|
|
|
3572 |
[84,] 1455 1870.167 159.4 157.5 1.2
|
|
|
3573 |
[85,] 1472 1871.583 110.0 110.1 0.1
|
|
|
3574 |
[86,] 1476 1871.917 90.3 90.4 0.1
|
|
|
3575 |
[87,] 1486 1872.750 103.5 102.6 0.9
|
|
|
3576 |
[88,] 1497 1873.667 47.5 47.1 0.8
|
|
|
3577 |
[89,] 1498 1873.750 47.4 47.1 0.6
|
|
|
3578 |
[90,] 1514 1875.083 22.2 21.5 3.2
|
|
|
3579 |
[91,] 1527 1876.167 31.2 30.6 1.9
|
|
|
3580 |
[92,] 1539 1877.167 11.7 11.9 1.7
|
|
|
3581 |
[93,] 1541 1877.333 21.2 21.6 1.9
|
|
|
3582 |
[94,] 1542 1877.417 13.4 14.2 5.8
|
|
|
3583 |
[95,] 1543 1877.500 5.9 6.0 1.7
|
|
|
3584 |
[96,] 1545 1877.667 16.4 16.9 3.0
|
|
|
3585 |
[97,] 1547 1877.833 14.5 14.2 2.1
|
|
|
3586 |
[98,] 1548 1877.917 2.3 2.2 4.4
|
|
|
3587 |
[99,] 1550 1878.083 6.0 6.6 9.5
|
|
|
3588 |
[100,] 1553 1878.333 5.8 5.9 1.7
|
|
|
3589 |
[101,] 1561 1879.000 0.8 1.0 22.2
|
|
|
3590 |
[102,] 1571 1879.833 12.9 13.1 1.5
|
|
|
3591 |
[103,] 1572 1879.917 7.2 7.3 1.4
|
|
|
3592 |
[104,] 1574 1880.083 27.5 27.2 1.1
|
|
|
3593 |
[105,] 1575 1880.167 19.5 19.3 1.0
|
|
|
3594 |
[106,] 1576 1880.250 19.3 19.5 1.0
|
|
|
3595 |
[107,] 1588 1881.250 51.7 51.6 0.2
|
|
|
3596 |
[108,] 1592 1881.583 58.0 58.4 0.7
|
|
|
3597 |
[109,] 1594 1881.750 64.0 64.4 0.6
|
|
|
3598 |
[110,] 1598 1882.083 69.3 69.5 0.3
|
|
|
3599 |
[111,] 1599 1882.167 67.5 66.8 1.0
|
|
|
3600 |
[112,] 1613 1883.333 32.1 31.5 1.9
|
|
|
3601 |
[113,] 1614 1883.417 76.5 76.3 0.3
|
|
|
3602 |
[114,] 1623 1884.167 86.8 87.5 0.8
|
|
|
3603 |
[115,] 1643 1885.833 33.3 30.9 7.5
|
|
|
3604 |
[116,] 1656 1886.917 12.4 13.0 4.7
|
|
|
3605 |
[117,] 1663 1887.500 23.3 23.4 0.4
|
|
|
3606 |
[118,] 1683 1889.167 7.0 6.7 4.4
|
|
|
3607 |
[119,] 1687 1889.500 9.7 9.4 3.1
|
|
|
3608 |
[120,] 1712 1891.583 33.2 33.0 0.6
|
|
|
3609 |
[121,] 1716 1891.917 32.3 32.5 0.6
|
|
|
3610 |
[122,] 1723 1892.500 76.8 76.5 0.4
|
|
|
3611 |
[123,] 1734 1893.417 88.2 89.9 1.9
|
|
|
3612 |
[124,] 1735 1893.500 88.8 88.6 0.2
|
|
|
3613 |
[125,] 1738 1893.750 79.7 80.0 0.4
|
|
|
3614 |
[126,] 1774 1896.750 28.4 28.7 1.1
|
|
|
3615 |
[127,] 1837 1902.000 5.2 5.5 5.6
|
|
|
3616 |
[128,] 2126 1926.083 70.0 69.9 0.1
|
|
|
3617 |
[129,] 2151 1928.167 85.4 85.5 0.1
|
|
|
3618 |
[130,] 2153 1928.333 76.9 77.0 0.1
|
|
|
3619 |
[131,] 2162 1929.083 64.1 62.8 2.0
|
|
|
3620 |
[132,] 2174 1930.083 49.2 49.9 1.4
|
|
|
3621 |
[133,] 2233 1935.000 18.9 18.6 1.6
|
|
|
3622 |
[134,] 2315 1941.833 38.3 38.4 0.3
|
|
|
3623 |
[135,] 2329 1943.000 12.4 12.5 0.8
|
|
|
3624 |
[136,] 2378 1947.083 113.4 133.4 16.2
|
|
|
3625 |
[137,] 2427 1951.167 59.9 55.9 6.9
|
|
|
3626 |
[138,] 2498 1957.083 130.2 130.3 0.1
|
|
|
3627 |
[139,] 2552 1961.583 55.9 55.8 0.2
|
|
|
3628 |
[140,] 2556 1961.917 40.0 39.9 0.3
|
|
|
3629 |
[141,] 2594 1965.083 14.2 14.3 0.7
|
|
|
3630 |
[142,] 2790 1981.417 90.0 90.9 1.0
|
|
|
3631 |
[143,] 2819 1983.833 33.3 33.4 0.3
|
| 56186 |
murdoch |
3632 |
>
|
| 87340 |
maechler |
3633 |
> ## How to recreate the "old" sunspot.month (R <= 3.0.3) =: sunspot.month.0
|
| 65253 |
maechler |
3634 |
> .sunspot.diff <- cbind(
|
|
|
3635 |
+ i = c(1202L, 1256L, 1258L, 1301L, 1407L, 1429L, 1452L, 1455L,
|
|
|
3636 |
+ 1663L, 2151L, 2329L, 2498L, 2594L, 2694L, 2819L),
|
|
|
3637 |
+ res10 = c(1L, 1L, 1L, -1L, -1L, -1L, 1L, -1L,
|
|
|
3638 |
+ 1L, 1L, 1L, 1L, 1L, 20L, 1L))
|
| 87340 |
maechler |
3639 |
> ssm0 <- sunspot.m2014[1:2988]
|
| 65253 |
maechler |
3640 |
> with(as.data.frame(.sunspot.diff), ssm0[i] <<- ssm0[i] - res10/10)
|
|
|
3641 |
> sunspot.month.0 <- ts(ssm0, start = 1749, frequency = 12)
|
| 87311 |
ripley |
3642 |
> stopifnot(length(sunspot.month.0) == 2988)
|
| 56186 |
murdoch |
3643 |
>
|
|
|
3644 |
>
|
| 65253 |
maechler |
3645 |
>
|
| 56186 |
murdoch |
3646 |
> cleanEx()
|
| 65249 |
maechler |
3647 |
> nameEx("sunspot.year")
|
|
|
3648 |
> ### * sunspot.year
|
|
|
3649 |
>
|
|
|
3650 |
> flush(stderr()); flush(stdout())
|
|
|
3651 |
>
|
|
|
3652 |
> ### Name: sunspot.year
|
|
|
3653 |
> ### Title: Yearly Sunspot Data, 1700-1988
|
|
|
3654 |
> ### Aliases: sunspot.year
|
|
|
3655 |
> ### Keywords: datasets
|
|
|
3656 |
>
|
|
|
3657 |
> ### ** Examples
|
|
|
3658 |
>
|
| 65261 |
ripley |
3659 |
> utils::str(sm <- sunspots)# the monthly version we keep unchanged
|
|
|
3660 |
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 ...
|
| 65249 |
maechler |
3661 |
> utils::str(sy <- sunspot.year)
|
|
|
3662 |
Time-Series [1:289] from 1700 to 1988: 5 11 16 23 36 58 29 20 10 8 ...
|
|
|
3663 |
> ## The common time interval
|
|
|
3664 |
> (t1 <- c(max(start(sm), start(sy)), 1)) # Jan 1749
|
|
|
3665 |
[1] 1749 1
|
| 87311 |
ripley |
3666 |
> (t2 <- c(min( end(sm)[1],end(sy)[1]), 12)) # Dec 1983 (will not be updated!)
|
| 65261 |
ripley |
3667 |
[1] 1983 12
|
| 65249 |
maechler |
3668 |
> s.m <- window(sm, start=t1, end=t2)
|
| 87311 |
ripley |
3669 |
> s.y <- window(sy, start=t1, end=t2[1])
|
| 65249 |
maechler |
3670 |
> stopifnot(length(s.y) * 12 == length(s.m),
|
|
|
3671 |
+ ## The yearly series *is* close to the averages of the monthly one:
|
| 80081 |
hornik |
3672 |
+ all.equal(s.y, aggregate(s.m, FUN = mean), tolerance = 0.0020))
|
| 65249 |
maechler |
3673 |
> ## NOTE: Strangely, correctly weighting the number of days per month
|
|
|
3674 |
> ## (using 28.25 for February) is *not* closer than the simple mean:
|
|
|
3675 |
> ndays <- c(31, 28.25, rep(c(31,30, 31,30, 31), 2))
|
| 65261 |
ripley |
3676 |
> all.equal(s.y, aggregate(s.m, FUN = mean)) # 0.0013
|
|
|
3677 |
[1] "Mean relative difference: 0.001312539"
|
|
|
3678 |
> all.equal(s.y, aggregate(s.m, FUN = weighted.mean, w = ndays)) # 0.0017
|
|
|
3679 |
[1] "Mean relative difference: 0.001692215"
|
| 65249 |
maechler |
3680 |
>
|
|
|
3681 |
>
|
|
|
3682 |
>
|
|
|
3683 |
> cleanEx()
|
| 56186 |
murdoch |
3684 |
> nameEx("sunspots")
|
|
|
3685 |
> ### * sunspots
|
|
|
3686 |
>
|
|
|
3687 |
> flush(stderr()); flush(stdout())
|
|
|
3688 |
>
|
|
|
3689 |
> ### Name: sunspots
|
|
|
3690 |
> ### Title: Monthly Sunspot Numbers, 1749-1983
|
|
|
3691 |
> ### Aliases: sunspots
|
|
|
3692 |
> ### Keywords: datasets
|
|
|
3693 |
>
|
|
|
3694 |
> ### ** Examples
|
|
|
3695 |
>
|
|
|
3696 |
> require(graphics)
|
|
|
3697 |
> plot(sunspots, main = "sunspots data", xlab = "Year",
|
|
|
3698 |
+ ylab = "Monthly sunspot numbers")
|
|
|
3699 |
>
|
|
|
3700 |
>
|
|
|
3701 |
>
|
|
|
3702 |
> cleanEx()
|
|
|
3703 |
> nameEx("swiss")
|
|
|
3704 |
> ### * swiss
|
|
|
3705 |
>
|
|
|
3706 |
> flush(stderr()); flush(stdout())
|
|
|
3707 |
>
|
|
|
3708 |
> ### Name: swiss
|
|
|
3709 |
> ### Title: Swiss Fertility and Socioeconomic Indicators (1888) Data
|
|
|
3710 |
> ### Aliases: swiss
|
|
|
3711 |
> ### Keywords: datasets
|
|
|
3712 |
>
|
|
|
3713 |
> ### ** Examples
|
|
|
3714 |
>
|
|
|
3715 |
> require(stats); require(graphics)
|
|
|
3716 |
> pairs(swiss, panel = panel.smooth, main = "swiss data",
|
|
|
3717 |
+ col = 3 + (swiss$Catholic > 50))
|
|
|
3718 |
> summary(lm(Fertility ~ . , data = swiss))
|
|
|
3719 |
|
|
|
3720 |
Call:
|
|
|
3721 |
lm(formula = Fertility ~ ., data = swiss)
|
|
|
3722 |
|
|
|
3723 |
Residuals:
|
|
|
3724 |
Min 1Q Median 3Q Max
|
|
|
3725 |
-15.2743 -5.2617 0.5032 4.1198 15.3213
|
|
|
3726 |
|
|
|
3727 |
Coefficients:
|
|
|
3728 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3729 |
(Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
|
|
|
3730 |
Agriculture -0.17211 0.07030 -2.448 0.01873 *
|
|
|
3731 |
Examination -0.25801 0.25388 -1.016 0.31546
|
|
|
3732 |
Education -0.87094 0.18303 -4.758 2.43e-05 ***
|
|
|
3733 |
Catholic 0.10412 0.03526 2.953 0.00519 **
|
|
|
3734 |
Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
|
|
|
3735 |
---
|
| 61435 |
ripley |
3736 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3737 |
|
|
|
3738 |
Residual standard error: 7.165 on 41 degrees of freedom
|
| 61435 |
ripley |
3739 |
Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
|
|
|
3740 |
F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10
|
| 56186 |
murdoch |
3741 |
|
|
|
3742 |
>
|
|
|
3743 |
>
|
|
|
3744 |
>
|
|
|
3745 |
> cleanEx()
|
|
|
3746 |
> nameEx("trees")
|
|
|
3747 |
> ### * trees
|
|
|
3748 |
>
|
|
|
3749 |
> flush(stderr()); flush(stdout())
|
|
|
3750 |
>
|
|
|
3751 |
> ### Name: trees
|
| 75009 |
ripley |
3752 |
> ### Title: Diameter, Height and Volume for Black Cherry Trees
|
| 56186 |
murdoch |
3753 |
> ### Aliases: trees
|
|
|
3754 |
> ### Keywords: datasets
|
|
|
3755 |
>
|
|
|
3756 |
> ### ** Examples
|
|
|
3757 |
>
|
|
|
3758 |
> require(stats); require(graphics)
|
|
|
3759 |
> pairs(trees, panel = panel.smooth, main = "trees data")
|
|
|
3760 |
> plot(Volume ~ Girth, data = trees, log = "xy")
|
|
|
3761 |
> coplot(log(Volume) ~ log(Girth) | Height, data = trees,
|
|
|
3762 |
+ panel = panel.smooth)
|
|
|
3763 |
> summary(fm1 <- lm(log(Volume) ~ log(Girth), data = trees))
|
|
|
3764 |
|
|
|
3765 |
Call:
|
|
|
3766 |
lm(formula = log(Volume) ~ log(Girth), data = trees)
|
|
|
3767 |
|
|
|
3768 |
Residuals:
|
|
|
3769 |
Min 1Q Median 3Q Max
|
|
|
3770 |
-0.205999 -0.068702 0.001011 0.072585 0.247963
|
|
|
3771 |
|
|
|
3772 |
Coefficients:
|
|
|
3773 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3774 |
(Intercept) -2.35332 0.23066 -10.20 4.18e-11 ***
|
|
|
3775 |
log(Girth) 2.19997 0.08983 24.49 < 2e-16 ***
|
|
|
3776 |
---
|
| 61435 |
ripley |
3777 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3778 |
|
|
|
3779 |
Residual standard error: 0.115 on 29 degrees of freedom
|
| 61435 |
ripley |
3780 |
Multiple R-squared: 0.9539, Adjusted R-squared: 0.9523
|
|
|
3781 |
F-statistic: 599.7 on 1 and 29 DF, p-value: < 2.2e-16
|
| 56186 |
murdoch |
3782 |
|
|
|
3783 |
> summary(fm2 <- update(fm1, ~ . + log(Height), data = trees))
|
|
|
3784 |
|
|
|
3785 |
Call:
|
|
|
3786 |
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
|
|
|
3787 |
|
|
|
3788 |
Residuals:
|
|
|
3789 |
Min 1Q Median 3Q Max
|
|
|
3790 |
-0.168561 -0.048488 0.002431 0.063637 0.129223
|
|
|
3791 |
|
|
|
3792 |
Coefficients:
|
|
|
3793 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3794 |
(Intercept) -6.63162 0.79979 -8.292 5.06e-09 ***
|
|
|
3795 |
log(Girth) 1.98265 0.07501 26.432 < 2e-16 ***
|
|
|
3796 |
log(Height) 1.11712 0.20444 5.464 7.81e-06 ***
|
|
|
3797 |
---
|
| 61435 |
ripley |
3798 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3799 |
|
|
|
3800 |
Residual standard error: 0.08139 on 28 degrees of freedom
|
| 61435 |
ripley |
3801 |
Multiple R-squared: 0.9777, Adjusted R-squared: 0.9761
|
|
|
3802 |
F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16
|
| 56186 |
murdoch |
3803 |
|
|
|
3804 |
> step(fm2)
|
|
|
3805 |
Start: AIC=-152.69
|
|
|
3806 |
log(Volume) ~ log(Girth) + log(Height)
|
|
|
3807 |
|
| 57809 |
ripley |
3808 |
Df Sum of Sq RSS AIC
|
|
|
3809 |
<none> 0.1855 -152.685
|
|
|
3810 |
- log(Height) 1 0.1978 0.3832 -132.185
|
|
|
3811 |
- log(Girth) 1 4.6275 4.8130 -53.743
|
| 56186 |
murdoch |
3812 |
|
|
|
3813 |
Call:
|
|
|
3814 |
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
|
|
|
3815 |
|
|
|
3816 |
Coefficients:
|
|
|
3817 |
(Intercept) log(Girth) log(Height)
|
|
|
3818 |
-6.632 1.983 1.117
|
|
|
3819 |
|
|
|
3820 |
> ## i.e., Volume ~= c * Height * Girth^2 seems reasonable
|
|
|
3821 |
>
|
|
|
3822 |
>
|
|
|
3823 |
>
|
|
|
3824 |
> cleanEx()
|
|
|
3825 |
> nameEx("uspop")
|
|
|
3826 |
> ### * uspop
|
|
|
3827 |
>
|
|
|
3828 |
> flush(stderr()); flush(stdout())
|
|
|
3829 |
>
|
|
|
3830 |
> ### Name: uspop
|
|
|
3831 |
> ### Title: Populations Recorded by the US Census
|
|
|
3832 |
> ### Aliases: uspop
|
|
|
3833 |
> ### Keywords: datasets
|
|
|
3834 |
>
|
|
|
3835 |
> ### ** Examples
|
|
|
3836 |
>
|
|
|
3837 |
> require(graphics)
|
|
|
3838 |
> plot(uspop, log = "y", main = "uspop data", xlab = "Year",
|
|
|
3839 |
+ ylab = "U.S. Population (millions)")
|
|
|
3840 |
>
|
|
|
3841 |
>
|
|
|
3842 |
>
|
|
|
3843 |
> cleanEx()
|
|
|
3844 |
> nameEx("volcano")
|
|
|
3845 |
> ### * volcano
|
|
|
3846 |
>
|
|
|
3847 |
> flush(stderr()); flush(stdout())
|
|
|
3848 |
>
|
|
|
3849 |
> ### Name: volcano
|
|
|
3850 |
> ### Title: Topographic Information on Auckland's Maunga Whau Volcano
|
|
|
3851 |
> ### Aliases: volcano
|
|
|
3852 |
> ### Keywords: datasets
|
|
|
3853 |
>
|
|
|
3854 |
> ### ** Examples
|
|
|
3855 |
>
|
|
|
3856 |
> require(grDevices); require(graphics)
|
|
|
3857 |
> filled.contour(volcano, color.palette = terrain.colors, asp = 1)
|
|
|
3858 |
> title(main = "volcano data: filled contour map")
|
|
|
3859 |
>
|
|
|
3860 |
>
|
|
|
3861 |
>
|
|
|
3862 |
> cleanEx()
|
|
|
3863 |
> nameEx("warpbreaks")
|
|
|
3864 |
> ### * warpbreaks
|
|
|
3865 |
>
|
|
|
3866 |
> flush(stderr()); flush(stdout())
|
|
|
3867 |
>
|
|
|
3868 |
> ### Name: warpbreaks
|
|
|
3869 |
> ### Title: The Number of Breaks in Yarn during Weaving
|
|
|
3870 |
> ### Aliases: warpbreaks
|
|
|
3871 |
> ### Keywords: datasets
|
|
|
3872 |
>
|
|
|
3873 |
> ### ** Examples
|
|
|
3874 |
>
|
|
|
3875 |
> require(stats); require(graphics)
|
|
|
3876 |
> summary(warpbreaks)
|
|
|
3877 |
breaks wool tension
|
|
|
3878 |
Min. :10.00 A:27 L:18
|
|
|
3879 |
1st Qu.:18.25 B:27 M:18
|
|
|
3880 |
Median :26.00 H:18
|
|
|
3881 |
Mean :28.15
|
|
|
3882 |
3rd Qu.:34.00
|
|
|
3883 |
Max. :70.00
|
| 61169 |
ripley |
3884 |
> opar <- par(mfrow = c(1, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
3885 |
> plot(breaks ~ tension, data = warpbreaks, col = "lightgray",
|
|
|
3886 |
+ varwidth = TRUE, subset = wool == "A", main = "Wool A")
|
|
|
3887 |
> plot(breaks ~ tension, data = warpbreaks, col = "lightgray",
|
|
|
3888 |
+ varwidth = TRUE, subset = wool == "B", main = "Wool B")
|
|
|
3889 |
> mtext("warpbreaks data", side = 3, outer = TRUE)
|
|
|
3890 |
> par(opar)
|
|
|
3891 |
> summary(fm1 <- lm(breaks ~ wool*tension, data = warpbreaks))
|
|
|
3892 |
|
|
|
3893 |
Call:
|
|
|
3894 |
lm(formula = breaks ~ wool * tension, data = warpbreaks)
|
|
|
3895 |
|
|
|
3896 |
Residuals:
|
|
|
3897 |
Min 1Q Median 3Q Max
|
|
|
3898 |
-19.5556 -6.8889 -0.6667 7.1944 25.4444
|
|
|
3899 |
|
|
|
3900 |
Coefficients:
|
|
|
3901 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3902 |
(Intercept) 44.556 3.647 12.218 2.43e-16 ***
|
|
|
3903 |
woolB -16.333 5.157 -3.167 0.002677 **
|
|
|
3904 |
tensionM -20.556 5.157 -3.986 0.000228 ***
|
|
|
3905 |
tensionH -20.000 5.157 -3.878 0.000320 ***
|
|
|
3906 |
woolB:tensionM 21.111 7.294 2.895 0.005698 **
|
|
|
3907 |
woolB:tensionH 10.556 7.294 1.447 0.154327
|
|
|
3908 |
---
|
| 61435 |
ripley |
3909 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3910 |
|
|
|
3911 |
Residual standard error: 10.94 on 48 degrees of freedom
|
| 61435 |
ripley |
3912 |
Multiple R-squared: 0.3778, Adjusted R-squared: 0.3129
|
|
|
3913 |
F-statistic: 5.828 on 5 and 48 DF, p-value: 0.0002772
|
| 56186 |
murdoch |
3914 |
|
|
|
3915 |
> anova(fm1)
|
|
|
3916 |
Analysis of Variance Table
|
|
|
3917 |
|
|
|
3918 |
Response: breaks
|
|
|
3919 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
3920 |
wool 1 450.7 450.67 3.7653 0.0582130 .
|
|
|
3921 |
tension 2 2034.3 1017.13 8.4980 0.0006926 ***
|
|
|
3922 |
wool:tension 2 1002.8 501.39 4.1891 0.0210442 *
|
|
|
3923 |
Residuals 48 5745.1 119.69
|
|
|
3924 |
---
|
| 61435 |
ripley |
3925 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3926 |
>
|
|
|
3927 |
>
|
|
|
3928 |
>
|
|
|
3929 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
3930 |
> cleanEx()
|
|
|
3931 |
> nameEx("women")
|
|
|
3932 |
> ### * women
|
|
|
3933 |
>
|
|
|
3934 |
> flush(stderr()); flush(stdout())
|
|
|
3935 |
>
|
|
|
3936 |
> ### Name: women
|
|
|
3937 |
> ### Title: Average Heights and Weights for American Women
|
|
|
3938 |
> ### Aliases: women
|
|
|
3939 |
> ### Keywords: datasets
|
|
|
3940 |
>
|
|
|
3941 |
> ### ** Examples
|
|
|
3942 |
>
|
|
|
3943 |
> require(graphics)
|
|
|
3944 |
> plot(women, xlab = "Height (in)", ylab = "Weight (lb)",
|
|
|
3945 |
+ main = "women data: American women aged 30-39")
|
|
|
3946 |
>
|
|
|
3947 |
>
|
|
|
3948 |
>
|
|
|
3949 |
> cleanEx()
|
|
|
3950 |
> nameEx("zCO2")
|
|
|
3951 |
> ### * zCO2
|
|
|
3952 |
>
|
|
|
3953 |
> flush(stderr()); flush(stdout())
|
|
|
3954 |
>
|
|
|
3955 |
> ### Name: CO2
|
|
|
3956 |
> ### Title: Carbon Dioxide Uptake in Grass Plants
|
|
|
3957 |
> ### Aliases: CO2
|
|
|
3958 |
> ### Keywords: datasets
|
|
|
3959 |
>
|
|
|
3960 |
> ### ** Examples
|
|
|
3961 |
>
|
|
|
3962 |
> require(stats); require(graphics)
|
|
|
3963 |
> ## Don't show:
|
|
|
3964 |
> options(show.nls.convergence=FALSE)
|
| 66943 |
maechler |
3965 |
> ## End(Don't show)
|
| 56186 |
murdoch |
3966 |
> coplot(uptake ~ conc | Plant, data = CO2, show.given = FALSE, type = "b")
|
|
|
3967 |
> ## fit the data for the first plant
|
|
|
3968 |
> fm1 <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0),
|
| 61169 |
ripley |
3969 |
+ data = CO2, subset = Plant == "Qn1")
|
| 56186 |
murdoch |
3970 |
> summary(fm1)
|
|
|
3971 |
|
|
|
3972 |
Formula: uptake ~ SSasymp(conc, Asym, lrc, c0)
|
|
|
3973 |
|
|
|
3974 |
Parameters:
|
|
|
3975 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3976 |
Asym 38.1398 0.9164 41.620 1.99e-06 ***
|
|
|
3977 |
lrc -34.2766 18.9661 -1.807 0.145
|
|
|
3978 |
c0 -4.3806 0.2042 -21.457 2.79e-05 ***
|
|
|
3979 |
---
|
| 61435 |
ripley |
3980 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3981 |
|
| 62439 |
ripley |
3982 |
Residual standard error: 1.663 on 4 degrees of freedom
|
| 56186 |
murdoch |
3983 |
|
|
|
3984 |
> ## fit each plant separately
|
|
|
3985 |
> fmlist <- list()
|
|
|
3986 |
> for (pp in levels(CO2$Plant)) {
|
|
|
3987 |
+ fmlist[[pp]] <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0),
|
|
|
3988 |
+ data = CO2, subset = Plant == pp)
|
|
|
3989 |
+ }
|
|
|
3990 |
> ## check the coefficients by plant
|
| 61157 |
ripley |
3991 |
> print(sapply(fmlist, coef), digits = 3)
|
| 56186 |
murdoch |
3992 |
Qn1 Qn2 Qn3 Qc1 Qc3 Qc2 Mn3 Mn2 Mn1 Mc2 Mc3
|
|
|
3993 |
Asym 38.14 42.87 44.23 36.43 40.68 39.82 28.48 32.13 34.08 13.56 18.54
|
|
|
3994 |
lrc -34.28 -29.66 -37.63 -9.90 -11.54 -51.53 -17.37 -29.04 -8.81 -1.98 -136.11
|
|
|
3995 |
c0 -4.38 -4.67 -4.49 -4.86 -4.95 -4.46 -4.59 -4.47 -5.06 -4.56 -3.47
|
|
|
3996 |
Mc1
|
|
|
3997 |
Asym 21.79
|
|
|
3998 |
lrc 2.45
|
|
|
3999 |
c0 -5.14
|
|
|
4000 |
>
|
|
|
4001 |
>
|
|
|
4002 |
>
|
|
|
4003 |
> ### * <FOOTER>
|
|
|
4004 |
> ###
|
| 73819 |
hornik |
4005 |
> cleanEx()
|
| 62439 |
ripley |
4006 |
> options(digits = 7L)
|
| 61787 |
ripley |
4007 |
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
|
| 89316 |
smeyer |
4008 |
Time elapsed: 1.734 0.076 1.809 0 0
|
| 56186 |
murdoch |
4009 |
> grDevices::dev.off()
|
|
|
4010 |
null device
|
|
|
4011 |
1
|
|
|
4012 |
> ###
|
|
|
4013 |
> ### Local variables: ***
|
|
|
4014 |
> ### mode: outline-minor ***
|
|
|
4015 |
> ### outline-regexp: "\\(> \\)?### [*]+" ***
|
|
|
4016 |
> ### End: ***
|
|
|
4017 |
> quit('no')
|
The R Project SVN
R