| 56186 |
murdoch |
1 |
|
| 63202 |
ripley |
2 |
R Under development (unstable) (2013-07-06 r63201) -- "Unsuffered Consequences"
|
| 61787 |
ripley |
3 |
Copyright (C) 2013 The R Foundation for Statistical Computing
|
| 57821 |
maechler |
4 |
Platform: x86_64-unknown-linux-gnu (64-bit)
|
| 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 |
Natural language support but running in an English locale
|
|
|
11 |
|
|
|
12 |
R is a collaborative project with many contributors.
|
|
|
13 |
Type 'contributors()' for more information and
|
|
|
14 |
'citation()' on how to cite R or R packages in publications.
|
|
|
15 |
|
|
|
16 |
Type 'demo()' for some demos, 'help()' for on-line help, or
|
|
|
17 |
'help.start()' for an HTML browser interface to help.
|
|
|
18 |
Type 'q()' to quit R.
|
|
|
19 |
|
|
|
20 |
> pkgname <- "datasets"
|
|
|
21 |
> source(file.path(R.home("share"), "R", "examples-header.R"))
|
|
|
22 |
> options(warn = 1)
|
|
|
23 |
> library('datasets')
|
|
|
24 |
>
|
| 61787 |
ripley |
25 |
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
|
| 56186 |
murdoch |
26 |
> cleanEx()
|
|
|
27 |
> nameEx("AirPassengers")
|
|
|
28 |
> ### * AirPassengers
|
|
|
29 |
>
|
|
|
30 |
> flush(stderr()); flush(stdout())
|
|
|
31 |
>
|
|
|
32 |
> ### Name: AirPassengers
|
|
|
33 |
> ### Title: Monthly Airline Passenger Numbers 1949-1960
|
|
|
34 |
> ### Aliases: AirPassengers
|
|
|
35 |
> ### Keywords: datasets
|
|
|
36 |
>
|
|
|
37 |
> ### ** Examples
|
|
|
38 |
>
|
|
|
39 |
> ## Not run:
|
|
|
40 |
> ##D ## These are quite slow and so not run by example(AirPassengers)
|
|
|
41 |
> ##D
|
|
|
42 |
> ##D ## The classic 'airline model', by full ML
|
|
|
43 |
> ##D (fit <- arima(log10(AirPassengers), c(0, 1, 1),
|
| 61169 |
ripley |
44 |
> ##D seasonal = list(order = c(0, 1, 1), period = 12)))
|
| 56186 |
murdoch |
45 |
> ##D update(fit, method = "CSS")
|
| 61157 |
ripley |
46 |
> ##D update(fit, x = window(log10(AirPassengers), start = 1954))
|
| 56186 |
murdoch |
47 |
> ##D pred <- predict(fit, n.ahead = 24)
|
|
|
48 |
> ##D tl <- pred$pred - 1.96 * pred$se
|
|
|
49 |
> ##D tu <- pred$pred + 1.96 * pred$se
|
| 61169 |
ripley |
50 |
> ##D ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty = c(1, 2, 2))
|
| 56186 |
murdoch |
51 |
> ##D
|
|
|
52 |
> ##D ## full ML fit is the same if the series is reversed, CSS fit is not
|
|
|
53 |
> ##D ap0 <- rev(log10(AirPassengers))
|
|
|
54 |
> ##D attributes(ap0) <- attributes(AirPassengers)
|
| 61169 |
ripley |
55 |
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
|
|
|
56 |
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),
|
| 56186 |
murdoch |
57 |
> ##D method = "CSS")
|
|
|
58 |
> ##D
|
|
|
59 |
> ##D ## Structural Time Series
|
|
|
60 |
> ##D ap <- log10(AirPassengers) - 2
|
| 61157 |
ripley |
61 |
> ##D (fit <- StructTS(ap, type = "BSM"))
|
| 61169 |
ripley |
62 |
> ##D par(mfrow = c(1, 2))
|
| 56186 |
murdoch |
63 |
> ##D plot(cbind(ap, fitted(fit)), plot.type = "single")
|
|
|
64 |
> ##D plot(cbind(ap, tsSmooth(fit)), plot.type = "single")
|
|
|
65 |
> ## End(Not run)
|
|
|
66 |
>
|
|
|
67 |
>
|
|
|
68 |
> cleanEx()
|
|
|
69 |
> nameEx("BOD")
|
|
|
70 |
> ### * BOD
|
|
|
71 |
>
|
|
|
72 |
> flush(stderr()); flush(stdout())
|
|
|
73 |
>
|
|
|
74 |
> ### Name: BOD
|
|
|
75 |
> ### Title: Biochemical Oxygen Demand
|
|
|
76 |
> ### Aliases: BOD
|
|
|
77 |
> ### Keywords: datasets
|
|
|
78 |
>
|
|
|
79 |
> ### ** Examples
|
|
|
80 |
>
|
|
|
81 |
> ## Don't show:
|
|
|
82 |
> options(show.nls.convergence=FALSE)
|
|
|
83 |
> old <- options(digits = 5)
|
|
|
84 |
> ## End Don't show
|
|
|
85 |
> require(stats)
|
|
|
86 |
> # simplest form of fitting a first-order model to these data
|
|
|
87 |
> fm1 <- nls(demand ~ A*(1-exp(-exp(lrc)*Time)), data = BOD,
|
|
|
88 |
+ start = c(A = 20, lrc = log(.35)))
|
|
|
89 |
> coef(fm1)
|
|
|
90 |
A lrc
|
|
|
91 |
19.14258 -0.63282
|
|
|
92 |
> fm1
|
|
|
93 |
Nonlinear regression model
|
| 61435 |
ripley |
94 |
model: demand ~ A * (1 - exp(-exp(lrc) * Time))
|
|
|
95 |
data: BOD
|
| 56186 |
murdoch |
96 |
A lrc
|
|
|
97 |
19.143 -0.633
|
|
|
98 |
residual sum-of-squares: 26
|
|
|
99 |
> # using the plinear algorithm
|
|
|
100 |
> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,
|
|
|
101 |
+ start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)
|
| 61435 |
ripley |
102 |
32.946 : -1.0498 22.1260
|
|
|
103 |
25.992 : -0.62572 19.10319
|
|
|
104 |
25.99 : -0.6327 19.1419
|
|
|
105 |
25.99 : -0.63282 19.14256
|
| 56186 |
murdoch |
106 |
> # using a self-starting model
|
|
|
107 |
> fm3 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD)
|
|
|
108 |
> summary(fm3)
|
|
|
109 |
|
|
|
110 |
Formula: demand ~ SSasympOrig(Time, A, lrc)
|
|
|
111 |
|
|
|
112 |
Parameters:
|
|
|
113 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
114 |
A 19.143 2.496 7.67 0.0016 **
|
|
|
115 |
lrc -0.633 0.382 -1.65 0.1733
|
|
|
116 |
---
|
| 61435 |
ripley |
117 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
118 |
|
| 62439 |
ripley |
119 |
Residual standard error: 2.55 on 4 degrees of freedom
|
| 56186 |
murdoch |
120 |
|
|
|
121 |
> ## Don't show:
|
|
|
122 |
> options(old)
|
|
|
123 |
> ## End Don't show
|
|
|
124 |
>
|
|
|
125 |
>
|
|
|
126 |
>
|
|
|
127 |
> cleanEx()
|
|
|
128 |
> nameEx("ChickWeight")
|
|
|
129 |
> ### * ChickWeight
|
|
|
130 |
>
|
|
|
131 |
> flush(stderr()); flush(stdout())
|
|
|
132 |
>
|
|
|
133 |
> ### Name: ChickWeight
|
|
|
134 |
> ### Title: Weight versus age of chicks on different diets
|
|
|
135 |
> ### Aliases: ChickWeight
|
|
|
136 |
> ### Keywords: datasets
|
|
|
137 |
>
|
|
|
138 |
> ### ** Examples
|
|
|
139 |
>
|
|
|
140 |
>
|
|
|
141 |
> cleanEx()
|
|
|
142 |
> nameEx("DNase")
|
|
|
143 |
> ### * DNase
|
|
|
144 |
>
|
|
|
145 |
> flush(stderr()); flush(stdout())
|
|
|
146 |
>
|
|
|
147 |
> ### Name: DNase
|
|
|
148 |
> ### Title: Elisa assay of DNase
|
|
|
149 |
> ### Aliases: DNase
|
|
|
150 |
> ### Keywords: datasets
|
|
|
151 |
>
|
|
|
152 |
> ### ** Examples
|
|
|
153 |
>
|
|
|
154 |
> require(stats); require(graphics)
|
|
|
155 |
> ## Don't show:
|
|
|
156 |
> options(show.nls.convergence=FALSE)
|
|
|
157 |
> ## End Don't show
|
|
|
158 |
> coplot(density ~ conc | Run, data = DNase,
|
|
|
159 |
+ show.given = FALSE, type = "b")
|
|
|
160 |
> coplot(density ~ log(conc) | Run, data = DNase,
|
|
|
161 |
+ show.given = FALSE, type = "b")
|
|
|
162 |
> ## fit a representative run
|
|
|
163 |
> fm1 <- nls(density ~ SSlogis( log(conc), Asym, xmid, scal ),
|
|
|
164 |
+ data = DNase, subset = Run == 1)
|
|
|
165 |
> ## compare with a four-parameter logistic
|
|
|
166 |
> fm2 <- nls(density ~ SSfpl( log(conc), A, B, xmid, scal ),
|
|
|
167 |
+ data = DNase, subset = Run == 1)
|
|
|
168 |
> summary(fm2)
|
|
|
169 |
|
|
|
170 |
Formula: density ~ SSfpl(log(conc), A, B, xmid, scal)
|
|
|
171 |
|
|
|
172 |
Parameters:
|
|
|
173 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
174 |
A -0.007897 0.017200 -0.459 0.654
|
|
|
175 |
B 2.377239 0.109516 21.707 5.35e-11 ***
|
|
|
176 |
xmid 1.507403 0.102080 14.767 4.65e-09 ***
|
|
|
177 |
scal 1.062579 0.056996 18.643 3.16e-10 ***
|
|
|
178 |
---
|
| 61435 |
ripley |
179 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
180 |
|
| 62439 |
ripley |
181 |
Residual standard error: 0.01981 on 12 degrees of freedom
|
| 56186 |
murdoch |
182 |
|
|
|
183 |
> anova(fm1, fm2)
|
|
|
184 |
Analysis of Variance Table
|
|
|
185 |
|
|
|
186 |
Model 1: density ~ SSlogis(log(conc), Asym, xmid, scal)
|
|
|
187 |
Model 2: density ~ SSfpl(log(conc), A, B, xmid, scal)
|
|
|
188 |
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
|
|
|
189 |
1 13 0.0047896
|
|
|
190 |
2 12 0.0047073 1 8.2314e-05 0.2098 0.6551
|
|
|
191 |
>
|
|
|
192 |
>
|
|
|
193 |
>
|
|
|
194 |
> cleanEx()
|
|
|
195 |
> nameEx("Formaldehyde")
|
|
|
196 |
> ### * Formaldehyde
|
|
|
197 |
>
|
|
|
198 |
> flush(stderr()); flush(stdout())
|
|
|
199 |
>
|
|
|
200 |
> ### Name: Formaldehyde
|
|
|
201 |
> ### Title: Determination of Formaldehyde
|
|
|
202 |
> ### Aliases: Formaldehyde
|
|
|
203 |
> ### Keywords: datasets
|
|
|
204 |
>
|
|
|
205 |
> ### ** Examples
|
|
|
206 |
>
|
|
|
207 |
> require(stats); require(graphics)
|
|
|
208 |
> plot(optden ~ carb, data = Formaldehyde,
|
|
|
209 |
+ xlab = "Carbohydrate (ml)", ylab = "Optical Density",
|
|
|
210 |
+ main = "Formaldehyde data", col = 4, las = 1)
|
|
|
211 |
> abline(fm1 <- lm(optden ~ carb, data = Formaldehyde))
|
|
|
212 |
> summary(fm1)
|
|
|
213 |
|
|
|
214 |
Call:
|
|
|
215 |
lm(formula = optden ~ carb, data = Formaldehyde)
|
|
|
216 |
|
|
|
217 |
Residuals:
|
|
|
218 |
1 2 3 4 5 6
|
|
|
219 |
-0.006714 0.001029 0.002771 0.007143 0.007514 -0.011743
|
|
|
220 |
|
|
|
221 |
Coefficients:
|
|
|
222 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
223 |
(Intercept) 0.005086 0.007834 0.649 0.552
|
|
|
224 |
carb 0.876286 0.013535 64.744 3.41e-07 ***
|
|
|
225 |
---
|
| 61435 |
ripley |
226 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
227 |
|
|
|
228 |
Residual standard error: 0.008649 on 4 degrees of freedom
|
| 61435 |
ripley |
229 |
Multiple R-squared: 0.999, Adjusted R-squared: 0.9988
|
|
|
230 |
F-statistic: 4192 on 1 and 4 DF, p-value: 3.409e-07
|
| 56186 |
murdoch |
231 |
|
| 61169 |
ripley |
232 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
233 |
> plot(fm1)
|
|
|
234 |
> par(opar)
|
|
|
235 |
>
|
|
|
236 |
>
|
|
|
237 |
>
|
|
|
238 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
239 |
> cleanEx()
|
|
|
240 |
> nameEx("HairEyeColor")
|
|
|
241 |
> ### * HairEyeColor
|
|
|
242 |
>
|
|
|
243 |
> flush(stderr()); flush(stdout())
|
|
|
244 |
>
|
|
|
245 |
> ### Name: HairEyeColor
|
|
|
246 |
> ### Title: Hair and Eye Color of Statistics Students
|
|
|
247 |
> ### Aliases: HairEyeColor
|
|
|
248 |
> ### Keywords: datasets
|
|
|
249 |
>
|
|
|
250 |
> ### ** Examples
|
|
|
251 |
>
|
|
|
252 |
> require(graphics)
|
|
|
253 |
> ## Full mosaic
|
|
|
254 |
> mosaicplot(HairEyeColor)
|
|
|
255 |
> ## Aggregate over sex (as in Snee's original data)
|
|
|
256 |
> x <- apply(HairEyeColor, c(1, 2), sum)
|
|
|
257 |
> x
|
|
|
258 |
Eye
|
|
|
259 |
Hair Brown Blue Hazel Green
|
|
|
260 |
Black 68 20 15 5
|
|
|
261 |
Brown 119 84 54 29
|
|
|
262 |
Red 26 17 14 14
|
|
|
263 |
Blond 7 94 10 16
|
|
|
264 |
> mosaicplot(x, main = "Relation between hair and eye color")
|
|
|
265 |
>
|
|
|
266 |
>
|
|
|
267 |
>
|
|
|
268 |
> cleanEx()
|
|
|
269 |
> nameEx("Harman23.cor")
|
|
|
270 |
> ### * Harman23.cor
|
|
|
271 |
>
|
|
|
272 |
> flush(stderr()); flush(stdout())
|
|
|
273 |
>
|
|
|
274 |
> ### Name: Harman23.cor
|
|
|
275 |
> ### Title: Harman Example 2.3
|
|
|
276 |
> ### Aliases: Harman23.cor
|
|
|
277 |
> ### Keywords: datasets
|
|
|
278 |
>
|
|
|
279 |
> ### ** Examples
|
|
|
280 |
>
|
|
|
281 |
> require(stats)
|
|
|
282 |
> (Harman23.FA <- factanal(factors = 1, covmat = Harman23.cor))
|
|
|
283 |
|
|
|
284 |
Call:
|
|
|
285 |
factanal(factors = 1, covmat = Harman23.cor)
|
|
|
286 |
|
|
|
287 |
Uniquenesses:
|
|
|
288 |
height arm.span forearm lower.leg weight
|
|
|
289 |
0.158 0.135 0.190 0.187 0.760
|
|
|
290 |
bitro.diameter chest.girth chest.width
|
|
|
291 |
0.829 0.877 0.801
|
|
|
292 |
|
|
|
293 |
Loadings:
|
|
|
294 |
Factor1
|
|
|
295 |
height 0.918
|
|
|
296 |
arm.span 0.930
|
|
|
297 |
forearm 0.900
|
|
|
298 |
lower.leg 0.902
|
|
|
299 |
weight 0.490
|
|
|
300 |
bitro.diameter 0.413
|
|
|
301 |
chest.girth 0.351
|
|
|
302 |
chest.width 0.446
|
|
|
303 |
|
|
|
304 |
Factor1
|
|
|
305 |
SS loadings 4.064
|
|
|
306 |
Proportion Var 0.508
|
|
|
307 |
|
|
|
308 |
Test of the hypothesis that 1 factor is sufficient.
|
|
|
309 |
The chi square statistic is 611.44 on 20 degrees of freedom.
|
|
|
310 |
The p-value is 1.12e-116
|
|
|
311 |
> for(factors in 2:4) print(update(Harman23.FA, factors = factors))
|
|
|
312 |
|
|
|
313 |
Call:
|
|
|
314 |
factanal(factors = factors, covmat = Harman23.cor)
|
|
|
315 |
|
|
|
316 |
Uniquenesses:
|
|
|
317 |
height arm.span forearm lower.leg weight
|
|
|
318 |
0.170 0.107 0.166 0.199 0.089
|
|
|
319 |
bitro.diameter chest.girth chest.width
|
|
|
320 |
0.364 0.416 0.537
|
|
|
321 |
|
|
|
322 |
Loadings:
|
|
|
323 |
Factor1 Factor2
|
|
|
324 |
height 0.865 0.287
|
|
|
325 |
arm.span 0.927 0.181
|
|
|
326 |
forearm 0.895 0.179
|
|
|
327 |
lower.leg 0.859 0.252
|
|
|
328 |
weight 0.233 0.925
|
|
|
329 |
bitro.diameter 0.194 0.774
|
|
|
330 |
chest.girth 0.134 0.752
|
|
|
331 |
chest.width 0.278 0.621
|
|
|
332 |
|
|
|
333 |
Factor1 Factor2
|
|
|
334 |
SS loadings 3.335 2.617
|
|
|
335 |
Proportion Var 0.417 0.327
|
|
|
336 |
Cumulative Var 0.417 0.744
|
|
|
337 |
|
|
|
338 |
Test of the hypothesis that 2 factors are sufficient.
|
|
|
339 |
The chi square statistic is 75.74 on 13 degrees of freedom.
|
|
|
340 |
The p-value is 6.94e-11
|
|
|
341 |
|
|
|
342 |
Call:
|
|
|
343 |
factanal(factors = factors, covmat = Harman23.cor)
|
|
|
344 |
|
|
|
345 |
Uniquenesses:
|
|
|
346 |
height arm.span forearm lower.leg weight
|
|
|
347 |
0.127 0.005 0.193 0.157 0.090
|
|
|
348 |
bitro.diameter chest.girth chest.width
|
|
|
349 |
0.359 0.411 0.490
|
|
|
350 |
|
|
|
351 |
Loadings:
|
|
|
352 |
Factor1 Factor2 Factor3
|
|
|
353 |
height 0.886 0.267 -0.130
|
|
|
354 |
arm.span 0.937 0.195 0.280
|
|
|
355 |
forearm 0.874 0.188
|
|
|
356 |
lower.leg 0.877 0.230 -0.145
|
|
|
357 |
weight 0.242 0.916 -0.106
|
|
|
358 |
bitro.diameter 0.193 0.777
|
|
|
359 |
chest.girth 0.137 0.755
|
|
|
360 |
chest.width 0.261 0.646 0.159
|
|
|
361 |
|
|
|
362 |
Factor1 Factor2 Factor3
|
|
|
363 |
SS loadings 3.379 2.628 0.162
|
|
|
364 |
Proportion Var 0.422 0.329 0.020
|
|
|
365 |
Cumulative Var 0.422 0.751 0.771
|
|
|
366 |
|
|
|
367 |
Test of the hypothesis that 3 factors are sufficient.
|
|
|
368 |
The chi square statistic is 22.81 on 7 degrees of freedom.
|
|
|
369 |
The p-value is 0.00184
|
|
|
370 |
|
|
|
371 |
Call:
|
|
|
372 |
factanal(factors = factors, covmat = Harman23.cor)
|
|
|
373 |
|
|
|
374 |
Uniquenesses:
|
|
|
375 |
height arm.span forearm lower.leg weight
|
|
|
376 |
0.137 0.005 0.191 0.116 0.138
|
|
|
377 |
bitro.diameter chest.girth chest.width
|
|
|
378 |
0.283 0.178 0.488
|
|
|
379 |
|
|
|
380 |
Loadings:
|
|
|
381 |
Factor1 Factor2 Factor3 Factor4
|
|
|
382 |
height 0.879 0.277 -0.115
|
|
|
383 |
arm.span 0.937 0.194 0.277
|
|
|
384 |
forearm 0.875 0.191
|
|
|
385 |
lower.leg 0.887 0.209 0.135 -0.188
|
|
|
386 |
weight 0.246 0.882 0.111 -0.109
|
|
|
387 |
bitro.diameter 0.187 0.822
|
|
|
388 |
chest.girth 0.117 0.729 0.526
|
|
|
389 |
chest.width 0.263 0.644 0.141
|
|
|
390 |
|
|
|
391 |
Factor1 Factor2 Factor3 Factor4
|
|
|
392 |
SS loadings 3.382 2.595 0.323 0.165
|
|
|
393 |
Proportion Var 0.423 0.324 0.040 0.021
|
|
|
394 |
Cumulative Var 0.423 0.747 0.787 0.808
|
|
|
395 |
|
|
|
396 |
Test of the hypothesis that 4 factors are sufficient.
|
|
|
397 |
The chi square statistic is 4.63 on 2 degrees of freedom.
|
|
|
398 |
The p-value is 0.0988
|
|
|
399 |
>
|
|
|
400 |
>
|
|
|
401 |
>
|
|
|
402 |
> cleanEx()
|
|
|
403 |
> nameEx("Harman74.cor")
|
|
|
404 |
> ### * Harman74.cor
|
|
|
405 |
>
|
|
|
406 |
> flush(stderr()); flush(stdout())
|
|
|
407 |
>
|
|
|
408 |
> ### Name: Harman74.cor
|
|
|
409 |
> ### Title: Harman Example 7.4
|
|
|
410 |
> ### Aliases: Harman74.cor
|
|
|
411 |
> ### Keywords: datasets
|
|
|
412 |
>
|
|
|
413 |
> ### ** Examples
|
|
|
414 |
>
|
|
|
415 |
> require(stats)
|
|
|
416 |
> (Harman74.FA <- factanal(factors = 1, covmat = Harman74.cor))
|
|
|
417 |
|
|
|
418 |
Call:
|
|
|
419 |
factanal(factors = 1, covmat = Harman74.cor)
|
|
|
420 |
|
|
|
421 |
Uniquenesses:
|
|
|
422 |
VisualPerception Cubes PaperFormBoard
|
|
|
423 |
0.677 0.866 0.830
|
|
|
424 |
Flags GeneralInformation PargraphComprehension
|
|
|
425 |
0.768 0.487 0.491
|
|
|
426 |
SentenceCompletion WordClassification WordMeaning
|
|
|
427 |
0.500 0.514 0.474
|
|
|
428 |
Addition Code CountingDots
|
|
|
429 |
0.818 0.731 0.824
|
|
|
430 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
431 |
0.681 0.833 0.863
|
|
|
432 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
433 |
0.775 0.812 0.778
|
|
|
434 |
FigureWord Deduction NumericalPuzzles
|
|
|
435 |
0.816 0.612 0.676
|
|
|
436 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
437 |
0.619 0.524 0.593
|
|
|
438 |
|
|
|
439 |
Loadings:
|
|
|
440 |
Factor1
|
|
|
441 |
VisualPerception 0.569
|
|
|
442 |
Cubes 0.366
|
|
|
443 |
PaperFormBoard 0.412
|
|
|
444 |
Flags 0.482
|
|
|
445 |
GeneralInformation 0.716
|
|
|
446 |
PargraphComprehension 0.713
|
|
|
447 |
SentenceCompletion 0.707
|
|
|
448 |
WordClassification 0.697
|
|
|
449 |
WordMeaning 0.725
|
|
|
450 |
Addition 0.426
|
|
|
451 |
Code 0.519
|
|
|
452 |
CountingDots 0.419
|
|
|
453 |
StraightCurvedCapitals 0.565
|
|
|
454 |
WordRecognition 0.408
|
|
|
455 |
NumberRecognition 0.370
|
|
|
456 |
FigureRecognition 0.474
|
|
|
457 |
ObjectNumber 0.434
|
|
|
458 |
NumberFigure 0.471
|
|
|
459 |
FigureWord 0.429
|
|
|
460 |
Deduction 0.623
|
|
|
461 |
NumericalPuzzles 0.569
|
|
|
462 |
ProblemReasoning 0.617
|
|
|
463 |
SeriesCompletion 0.690
|
|
|
464 |
ArithmeticProblems 0.638
|
|
|
465 |
|
|
|
466 |
Factor1
|
|
|
467 |
SS loadings 7.438
|
|
|
468 |
Proportion Var 0.310
|
|
|
469 |
|
|
|
470 |
Test of the hypothesis that 1 factor is sufficient.
|
|
|
471 |
The chi square statistic is 622.91 on 252 degrees of freedom.
|
|
|
472 |
The p-value is 2.28e-33
|
|
|
473 |
> for(factors in 2:5) print(update(Harman74.FA, factors = factors))
|
|
|
474 |
|
|
|
475 |
Call:
|
|
|
476 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
477 |
|
|
|
478 |
Uniquenesses:
|
|
|
479 |
VisualPerception Cubes PaperFormBoard
|
|
|
480 |
0.650 0.864 0.844
|
|
|
481 |
Flags GeneralInformation PargraphComprehension
|
|
|
482 |
0.778 0.375 0.316
|
|
|
483 |
SentenceCompletion WordClassification WordMeaning
|
|
|
484 |
0.319 0.503 0.258
|
|
|
485 |
Addition Code CountingDots
|
|
|
486 |
0.670 0.608 0.581
|
|
|
487 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
488 |
0.567 0.832 0.850
|
|
|
489 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
490 |
0.743 0.770 0.625
|
|
|
491 |
FigureWord Deduction NumericalPuzzles
|
|
|
492 |
0.792 0.629 0.579
|
|
|
493 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
494 |
0.634 0.539 0.553
|
|
|
495 |
|
|
|
496 |
Loadings:
|
|
|
497 |
Factor1 Factor2
|
|
|
498 |
VisualPerception 0.506 0.306
|
|
|
499 |
Cubes 0.304 0.209
|
|
|
500 |
PaperFormBoard 0.297 0.260
|
|
|
501 |
Flags 0.327 0.339
|
|
|
502 |
GeneralInformation 0.240 0.753
|
|
|
503 |
PargraphComprehension 0.171 0.809
|
|
|
504 |
SentenceCompletion 0.163 0.809
|
|
|
505 |
WordClassification 0.344 0.615
|
|
|
506 |
WordMeaning 0.148 0.849
|
|
|
507 |
Addition 0.563 0.115
|
|
|
508 |
Code 0.591 0.207
|
|
|
509 |
CountingDots 0.647
|
|
|
510 |
StraightCurvedCapitals 0.612 0.241
|
|
|
511 |
WordRecognition 0.315 0.263
|
|
|
512 |
NumberRecognition 0.328 0.205
|
|
|
513 |
FigureRecognition 0.457 0.218
|
|
|
514 |
ObjectNumber 0.431 0.209
|
|
|
515 |
NumberFigure 0.601 0.116
|
|
|
516 |
FigureWord 0.399 0.222
|
|
|
517 |
Deduction 0.379 0.477
|
|
|
518 |
NumericalPuzzles 0.604 0.237
|
|
|
519 |
ProblemReasoning 0.390 0.462
|
|
|
520 |
SeriesCompletion 0.486 0.474
|
|
|
521 |
ArithmeticProblems 0.544 0.389
|
|
|
522 |
|
|
|
523 |
Factor1 Factor2
|
|
|
524 |
SS loadings 4.573 4.548
|
|
|
525 |
Proportion Var 0.191 0.190
|
|
|
526 |
Cumulative Var 0.191 0.380
|
|
|
527 |
|
|
|
528 |
Test of the hypothesis that 2 factors are sufficient.
|
|
|
529 |
The chi square statistic is 420.24 on 229 degrees of freedom.
|
|
|
530 |
The p-value is 2.01e-13
|
|
|
531 |
|
|
|
532 |
Call:
|
|
|
533 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
534 |
|
|
|
535 |
Uniquenesses:
|
|
|
536 |
VisualPerception Cubes PaperFormBoard
|
|
|
537 |
0.500 0.793 0.662
|
|
|
538 |
Flags GeneralInformation PargraphComprehension
|
|
|
539 |
0.694 0.352 0.316
|
|
|
540 |
SentenceCompletion WordClassification WordMeaning
|
|
|
541 |
0.300 0.502 0.256
|
|
|
542 |
Addition Code CountingDots
|
|
|
543 |
0.200 0.586 0.494
|
|
|
544 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
545 |
0.569 0.838 0.848
|
|
|
546 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
547 |
0.643 0.780 0.635
|
|
|
548 |
FigureWord Deduction NumericalPuzzles
|
|
|
549 |
0.788 0.590 0.580
|
|
|
550 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
551 |
0.597 0.498 0.500
|
|
|
552 |
|
|
|
553 |
Loadings:
|
|
|
554 |
Factor1 Factor2 Factor3
|
|
|
555 |
VisualPerception 0.176 0.656 0.198
|
|
|
556 |
Cubes 0.122 0.428
|
|
|
557 |
PaperFormBoard 0.145 0.563
|
|
|
558 |
Flags 0.239 0.487 0.107
|
|
|
559 |
GeneralInformation 0.745 0.191 0.237
|
|
|
560 |
PargraphComprehension 0.780 0.249 0.118
|
|
|
561 |
SentenceCompletion 0.802 0.175 0.160
|
|
|
562 |
WordClassification 0.571 0.327 0.256
|
|
|
563 |
WordMeaning 0.821 0.248
|
|
|
564 |
Addition 0.162 -0.118 0.871
|
|
|
565 |
Code 0.198 0.219 0.572
|
|
|
566 |
CountingDots 0.179 0.688
|
|
|
567 |
StraightCurvedCapitals 0.190 0.381 0.499
|
|
|
568 |
WordRecognition 0.231 0.253 0.210
|
|
|
569 |
NumberRecognition 0.158 0.299 0.195
|
|
|
570 |
FigureRecognition 0.108 0.557 0.186
|
|
|
571 |
ObjectNumber 0.178 0.267 0.342
|
|
|
572 |
NumberFigure 0.427 0.424
|
|
|
573 |
FigureWord 0.167 0.355 0.240
|
|
|
574 |
Deduction 0.392 0.472 0.181
|
|
|
575 |
NumericalPuzzles 0.178 0.406 0.473
|
|
|
576 |
ProblemReasoning 0.382 0.473 0.182
|
|
|
577 |
SeriesCompletion 0.379 0.528 0.283
|
|
|
578 |
ArithmeticProblems 0.377 0.226 0.554
|
|
|
579 |
|
|
|
580 |
Factor1 Factor2 Factor3
|
|
|
581 |
SS loadings 3.802 3.488 3.186
|
|
|
582 |
Proportion Var 0.158 0.145 0.133
|
|
|
583 |
Cumulative Var 0.158 0.304 0.436
|
|
|
584 |
|
|
|
585 |
Test of the hypothesis that 3 factors are sufficient.
|
|
|
586 |
The chi square statistic is 295.59 on 207 degrees of freedom.
|
|
|
587 |
The p-value is 5.12e-05
|
|
|
588 |
|
|
|
589 |
Call:
|
|
|
590 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
591 |
|
|
|
592 |
Uniquenesses:
|
|
|
593 |
VisualPerception Cubes PaperFormBoard
|
|
|
594 |
0.438 0.780 0.644
|
|
|
595 |
Flags GeneralInformation PargraphComprehension
|
|
|
596 |
0.651 0.352 0.312
|
|
|
597 |
SentenceCompletion WordClassification WordMeaning
|
|
|
598 |
0.283 0.485 0.257
|
|
|
599 |
Addition Code CountingDots
|
|
|
600 |
0.240 0.551 0.435
|
|
|
601 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
602 |
0.491 0.646 0.696
|
|
|
603 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
604 |
0.549 0.598 0.593
|
|
|
605 |
FigureWord Deduction NumericalPuzzles
|
|
|
606 |
0.762 0.592 0.583
|
|
|
607 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
608 |
0.601 0.497 0.500
|
|
|
609 |
|
|
|
610 |
Loadings:
|
|
|
611 |
Factor1 Factor2 Factor3 Factor4
|
|
|
612 |
VisualPerception 0.160 0.689 0.187 0.160
|
|
|
613 |
Cubes 0.117 0.436
|
|
|
614 |
PaperFormBoard 0.137 0.570 0.110
|
|
|
615 |
Flags 0.233 0.527
|
|
|
616 |
GeneralInformation 0.739 0.185 0.213 0.150
|
|
|
617 |
PargraphComprehension 0.767 0.205 0.233
|
|
|
618 |
SentenceCompletion 0.806 0.197 0.153
|
|
|
619 |
WordClassification 0.569 0.339 0.242 0.132
|
|
|
620 |
WordMeaning 0.806 0.201 0.227
|
|
|
621 |
Addition 0.167 -0.118 0.831 0.166
|
|
|
622 |
Code 0.180 0.120 0.512 0.374
|
|
|
623 |
CountingDots 0.210 0.716
|
|
|
624 |
StraightCurvedCapitals 0.188 0.438 0.525
|
|
|
625 |
WordRecognition 0.197 0.553
|
|
|
626 |
NumberRecognition 0.122 0.116 0.520
|
|
|
627 |
FigureRecognition 0.408 0.525
|
|
|
628 |
ObjectNumber 0.142 0.219 0.574
|
|
|
629 |
NumberFigure 0.293 0.336 0.456
|
|
|
630 |
FigureWord 0.148 0.239 0.161 0.365
|
|
|
631 |
Deduction 0.378 0.402 0.118 0.301
|
|
|
632 |
NumericalPuzzles 0.175 0.381 0.438 0.223
|
|
|
633 |
ProblemReasoning 0.366 0.399 0.123 0.301
|
|
|
634 |
SeriesCompletion 0.369 0.500 0.244 0.239
|
|
|
635 |
ArithmeticProblems 0.370 0.158 0.496 0.304
|
|
|
636 |
|
|
|
637 |
Factor1 Factor2 Factor3 Factor4
|
|
|
638 |
SS loadings 3.647 2.872 2.657 2.290
|
|
|
639 |
Proportion Var 0.152 0.120 0.111 0.095
|
|
|
640 |
Cumulative Var 0.152 0.272 0.382 0.478
|
|
|
641 |
|
|
|
642 |
Test of the hypothesis that 4 factors are sufficient.
|
|
|
643 |
The chi square statistic is 226.68 on 186 degrees of freedom.
|
|
|
644 |
The p-value is 0.0224
|
|
|
645 |
|
|
|
646 |
Call:
|
|
|
647 |
factanal(factors = factors, covmat = Harman74.cor)
|
|
|
648 |
|
|
|
649 |
Uniquenesses:
|
|
|
650 |
VisualPerception Cubes PaperFormBoard
|
|
|
651 |
0.450 0.781 0.639
|
|
|
652 |
Flags GeneralInformation PargraphComprehension
|
|
|
653 |
0.649 0.357 0.288
|
|
|
654 |
SentenceCompletion WordClassification WordMeaning
|
|
|
655 |
0.277 0.485 0.262
|
|
|
656 |
Addition Code CountingDots
|
|
|
657 |
0.215 0.386 0.444
|
|
|
658 |
StraightCurvedCapitals WordRecognition NumberRecognition
|
|
|
659 |
0.256 0.639 0.706
|
|
|
660 |
FigureRecognition ObjectNumber NumberFigure
|
|
|
661 |
0.550 0.614 0.596
|
|
|
662 |
FigureWord Deduction NumericalPuzzles
|
|
|
663 |
0.764 0.521 0.564
|
|
|
664 |
ProblemReasoning SeriesCompletion ArithmeticProblems
|
|
|
665 |
0.580 0.442 0.478
|
|
|
666 |
|
|
|
667 |
Loadings:
|
|
|
668 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
669 |
VisualPerception 0.161 0.658 0.136 0.182 0.199
|
|
|
670 |
Cubes 0.113 0.435 0.107
|
|
|
671 |
PaperFormBoard 0.135 0.562 0.107 0.116
|
|
|
672 |
Flags 0.231 0.533
|
|
|
673 |
GeneralInformation 0.736 0.188 0.192 0.162
|
|
|
674 |
PargraphComprehension 0.775 0.187 0.251 0.113
|
|
|
675 |
SentenceCompletion 0.809 0.208 0.136
|
|
|
676 |
WordClassification 0.568 0.348 0.223 0.131
|
|
|
677 |
WordMeaning 0.800 0.215 0.224
|
|
|
678 |
Addition 0.175 -0.100 0.844 0.176
|
|
|
679 |
Code 0.185 0.438 0.451 0.426
|
|
|
680 |
CountingDots 0.222 0.690 0.101 0.140
|
|
|
681 |
StraightCurvedCapitals 0.186 0.425 0.458 0.559
|
|
|
682 |
WordRecognition 0.197 0.557
|
|
|
683 |
NumberRecognition 0.121 0.130 0.508
|
|
|
684 |
FigureRecognition 0.400 0.529
|
|
|
685 |
ObjectNumber 0.145 0.208 0.562
|
|
|
686 |
NumberFigure 0.306 0.325 0.452
|
|
|
687 |
FigureWord 0.147 0.242 0.145 0.364
|
|
|
688 |
Deduction 0.370 0.452 0.139 0.287 -0.190
|
|
|
689 |
NumericalPuzzles 0.170 0.402 0.439 0.230
|
|
|
690 |
ProblemReasoning 0.358 0.423 0.126 0.302
|
|
|
691 |
SeriesCompletion 0.360 0.549 0.256 0.223 -0.107
|
|
|
692 |
ArithmeticProblems 0.371 0.185 0.502 0.307
|
|
|
693 |
|
|
|
694 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
695 |
SS loadings 3.632 2.964 2.456 2.345 0.663
|
|
|
696 |
Proportion Var 0.151 0.124 0.102 0.098 0.028
|
|
|
697 |
Cumulative Var 0.151 0.275 0.377 0.475 0.503
|
|
|
698 |
|
|
|
699 |
Test of the hypothesis that 5 factors are sufficient.
|
|
|
700 |
The chi square statistic is 186.82 on 166 degrees of freedom.
|
|
|
701 |
The p-value is 0.128
|
|
|
702 |
> Harman74.FA <- factanal(factors = 5, covmat = Harman74.cor,
|
| 61157 |
ripley |
703 |
+ rotation = "promax")
|
| 56186 |
murdoch |
704 |
> print(Harman74.FA$loadings, sort = TRUE)
|
|
|
705 |
|
|
|
706 |
Loadings:
|
|
|
707 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
708 |
VisualPerception 0.831 -0.127 0.230
|
|
|
709 |
Cubes 0.534
|
|
|
710 |
PaperFormBoard 0.736 -0.290 0.136
|
|
|
711 |
Flags 0.647 -0.104
|
|
|
712 |
SeriesCompletion 0.555 0.126 0.127
|
|
|
713 |
GeneralInformation 0.764
|
|
|
714 |
PargraphComprehension 0.845 -0.140 0.140
|
|
|
715 |
SentenceCompletion 0.872 -0.140
|
|
|
716 |
WordClassification 0.277 0.505 0.104
|
|
|
717 |
WordMeaning 0.846 -0.108
|
|
|
718 |
Addition -0.334 1.012
|
|
|
719 |
CountingDots 0.206 -0.200 0.722 0.185
|
|
|
720 |
ArithmeticProblems 0.197 0.500 0.139
|
|
|
721 |
WordRecognition -0.126 0.127 -0.103 0.657
|
|
|
722 |
NumberRecognition 0.568
|
|
|
723 |
FigureRecognition 0.399 -0.142 -0.207 0.562
|
|
|
724 |
ObjectNumber -0.108 0.107 0.613
|
|
|
725 |
StraightCurvedCapitals 0.542 0.247 0.618
|
|
|
726 |
Code 0.112 0.288 0.486 0.424
|
|
|
727 |
NumberFigure 0.255 -0.230 0.211 0.413
|
|
|
728 |
FigureWord 0.187 0.347
|
|
|
729 |
Deduction 0.404 0.169 0.117 -0.203
|
|
|
730 |
NumericalPuzzles 0.393 0.368
|
|
|
731 |
ProblemReasoning 0.381 0.188 0.169
|
|
|
732 |
|
|
|
733 |
Factor1 Factor2 Factor3 Factor4 Factor5
|
|
|
734 |
SS loadings 3.529 3.311 2.367 2.109 0.762
|
|
|
735 |
Proportion Var 0.147 0.138 0.099 0.088 0.032
|
|
|
736 |
Cumulative Var 0.147 0.285 0.384 0.471 0.503
|
|
|
737 |
>
|
|
|
738 |
>
|
|
|
739 |
>
|
|
|
740 |
> cleanEx()
|
|
|
741 |
> nameEx("InsectSprays")
|
|
|
742 |
> ### * InsectSprays
|
|
|
743 |
>
|
|
|
744 |
> flush(stderr()); flush(stdout())
|
|
|
745 |
>
|
|
|
746 |
> ### Name: InsectSprays
|
|
|
747 |
> ### Title: Effectiveness of Insect Sprays
|
|
|
748 |
> ### Aliases: InsectSprays
|
|
|
749 |
> ### Keywords: datasets
|
|
|
750 |
>
|
|
|
751 |
> ### ** Examples
|
|
|
752 |
>
|
|
|
753 |
> require(stats); require(graphics)
|
|
|
754 |
> boxplot(count ~ spray, data = InsectSprays,
|
|
|
755 |
+ xlab = "Type of spray", ylab = "Insect count",
|
|
|
756 |
+ main = "InsectSprays data", varwidth = TRUE, col = "lightgray")
|
|
|
757 |
> fm1 <- aov(count ~ spray, data = InsectSprays)
|
|
|
758 |
> summary(fm1)
|
| 57044 |
ripley |
759 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
760 |
spray 5 2669 533.8 34.7 <2e-16 ***
|
|
|
761 |
Residuals 66 1015 15.4
|
| 56186 |
murdoch |
762 |
---
|
| 61435 |
ripley |
763 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 61169 |
ripley |
764 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
765 |
> plot(fm1)
|
|
|
766 |
> fm2 <- aov(sqrt(count) ~ spray, data = InsectSprays)
|
|
|
767 |
> summary(fm2)
|
| 57044 |
ripley |
768 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
769 |
spray 5 88.44 17.688 44.8 <2e-16 ***
|
|
|
770 |
Residuals 66 26.06 0.395
|
| 56186 |
murdoch |
771 |
---
|
| 61435 |
ripley |
772 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
773 |
> plot(fm2)
|
|
|
774 |
> par(opar)
|
|
|
775 |
>
|
|
|
776 |
>
|
|
|
777 |
>
|
|
|
778 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
779 |
> cleanEx()
|
|
|
780 |
> nameEx("JohnsonJohnson")
|
|
|
781 |
> ### * JohnsonJohnson
|
|
|
782 |
>
|
|
|
783 |
> flush(stderr()); flush(stdout())
|
|
|
784 |
>
|
|
|
785 |
> ### Name: JohnsonJohnson
|
|
|
786 |
> ### Title: Quarterly Earnings per Johnson & Johnson Share
|
|
|
787 |
> ### Aliases: JohnsonJohnson
|
|
|
788 |
> ### Keywords: datasets
|
|
|
789 |
>
|
|
|
790 |
> ### ** Examples
|
|
|
791 |
>
|
|
|
792 |
>
|
|
|
793 |
> cleanEx()
|
|
|
794 |
> nameEx("LifeCycleSavings")
|
|
|
795 |
> ### * LifeCycleSavings
|
|
|
796 |
>
|
|
|
797 |
> flush(stderr()); flush(stdout())
|
|
|
798 |
>
|
|
|
799 |
> ### Name: LifeCycleSavings
|
|
|
800 |
> ### Title: Intercountry Life-Cycle Savings Data
|
|
|
801 |
> ### Aliases: LifeCycleSavings
|
|
|
802 |
> ### Keywords: datasets
|
|
|
803 |
>
|
|
|
804 |
> ### ** Examples
|
|
|
805 |
>
|
|
|
806 |
> require(stats); require(graphics)
|
|
|
807 |
> pairs(LifeCycleSavings, panel = panel.smooth,
|
|
|
808 |
+ main = "LifeCycleSavings data")
|
|
|
809 |
> fm1 <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
|
|
|
810 |
> summary(fm1)
|
|
|
811 |
|
|
|
812 |
Call:
|
|
|
813 |
lm(formula = sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
|
|
|
814 |
|
|
|
815 |
Residuals:
|
|
|
816 |
Min 1Q Median 3Q Max
|
|
|
817 |
-8.2422 -2.6857 -0.2488 2.4280 9.7509
|
|
|
818 |
|
|
|
819 |
Coefficients:
|
|
|
820 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
821 |
(Intercept) 28.5660865 7.3545161 3.884 0.000334 ***
|
|
|
822 |
pop15 -0.4611931 0.1446422 -3.189 0.002603 **
|
|
|
823 |
pop75 -1.6914977 1.0835989 -1.561 0.125530
|
|
|
824 |
dpi -0.0003369 0.0009311 -0.362 0.719173
|
|
|
825 |
ddpi 0.4096949 0.1961971 2.088 0.042471 *
|
|
|
826 |
---
|
| 61435 |
ripley |
827 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
828 |
|
|
|
829 |
Residual standard error: 3.803 on 45 degrees of freedom
|
| 61435 |
ripley |
830 |
Multiple R-squared: 0.3385, Adjusted R-squared: 0.2797
|
|
|
831 |
F-statistic: 5.756 on 4 and 45 DF, p-value: 0.0007904
|
| 56186 |
murdoch |
832 |
|
|
|
833 |
>
|
|
|
834 |
>
|
|
|
835 |
>
|
|
|
836 |
> cleanEx()
|
|
|
837 |
> nameEx("Loblolly")
|
|
|
838 |
> ### * Loblolly
|
|
|
839 |
>
|
|
|
840 |
> flush(stderr()); flush(stdout())
|
|
|
841 |
>
|
|
|
842 |
> ### Name: Loblolly
|
|
|
843 |
> ### Title: Growth of Loblolly pine trees
|
|
|
844 |
> ### Aliases: Loblolly
|
|
|
845 |
> ### Keywords: datasets
|
|
|
846 |
>
|
|
|
847 |
> ### ** Examples
|
|
|
848 |
>
|
|
|
849 |
> require(stats); require(graphics)
|
|
|
850 |
> plot(height ~ age, data = Loblolly, subset = Seed == 329,
|
|
|
851 |
+ xlab = "Tree age (yr)", las = 1,
|
|
|
852 |
+ ylab = "Tree height (ft)",
|
|
|
853 |
+ main = "Loblolly data and fitted curve (Seed 329 only)")
|
|
|
854 |
> fm1 <- nls(height ~ SSasymp(age, Asym, R0, lrc),
|
|
|
855 |
+ data = Loblolly, subset = Seed == 329)
|
|
|
856 |
> age <- seq(0, 30, length.out = 101)
|
|
|
857 |
> lines(age, predict(fm1, list(age = age)))
|
|
|
858 |
>
|
|
|
859 |
>
|
|
|
860 |
>
|
|
|
861 |
> cleanEx()
|
|
|
862 |
> nameEx("Nile")
|
|
|
863 |
> ### * Nile
|
|
|
864 |
>
|
|
|
865 |
> flush(stderr()); flush(stdout())
|
|
|
866 |
>
|
|
|
867 |
> ### Name: Nile
|
|
|
868 |
> ### Title: Flow of the River Nile
|
|
|
869 |
> ### Aliases: Nile
|
|
|
870 |
> ### Keywords: datasets
|
|
|
871 |
>
|
|
|
872 |
> ### ** Examples
|
|
|
873 |
>
|
|
|
874 |
> require(stats); require(graphics)
|
| 61169 |
ripley |
875 |
> par(mfrow = c(2, 2))
|
| 56186 |
murdoch |
876 |
> plot(Nile)
|
|
|
877 |
> acf(Nile)
|
|
|
878 |
> pacf(Nile)
|
|
|
879 |
> ar(Nile) # selects order 2
|
|
|
880 |
|
|
|
881 |
Call:
|
|
|
882 |
ar(x = Nile)
|
|
|
883 |
|
|
|
884 |
Coefficients:
|
|
|
885 |
1 2
|
|
|
886 |
0.4081 0.1812
|
|
|
887 |
|
| 61435 |
ripley |
888 |
Order selected 2 sigma^2 estimated as 21247
|
| 56186 |
murdoch |
889 |
> cpgram(ar(Nile)$resid)
|
| 61169 |
ripley |
890 |
> par(mfrow = c(1, 1))
|
| 56186 |
murdoch |
891 |
> arima(Nile, c(2, 0, 0))
|
|
|
892 |
|
|
|
893 |
Call:
|
|
|
894 |
arima(x = Nile, order = c(2, 0, 0))
|
|
|
895 |
|
|
|
896 |
Coefficients:
|
|
|
897 |
ar1 ar2 intercept
|
|
|
898 |
0.4096 0.1987 919.8397
|
|
|
899 |
s.e. 0.0974 0.0990 35.6410
|
|
|
900 |
|
|
|
901 |
sigma^2 estimated as 20291: log likelihood = -637.98, aic = 1283.96
|
|
|
902 |
>
|
|
|
903 |
> ## Now consider missing values, following Durbin & Koopman
|
|
|
904 |
> NileNA <- Nile
|
|
|
905 |
> NileNA[c(21:40, 61:80)] <- NA
|
|
|
906 |
> arima(NileNA, c(2, 0, 0))
|
|
|
907 |
|
|
|
908 |
Call:
|
|
|
909 |
arima(x = NileNA, order = c(2, 0, 0))
|
|
|
910 |
|
|
|
911 |
Coefficients:
|
|
|
912 |
ar1 ar2 intercept
|
|
|
913 |
0.3622 0.1678 918.3103
|
|
|
914 |
s.e. 0.1273 0.1323 39.5037
|
|
|
915 |
|
|
|
916 |
sigma^2 estimated as 23676: log likelihood = -387.7, aic = 783.41
|
|
|
917 |
> plot(NileNA)
|
|
|
918 |
> pred <-
|
| 61169 |
ripley |
919 |
+ predict(arima(window(NileNA, 1871, 1890), c(2, 0, 0)), n.ahead = 20)
|
| 56186 |
murdoch |
920 |
> lines(pred$pred, lty = 3, col = "red")
|
| 61157 |
ripley |
921 |
> lines(pred$pred + 2*pred$se, lty = 2, col = "blue")
|
|
|
922 |
> lines(pred$pred - 2*pred$se, lty = 2, col = "blue")
|
| 56186 |
murdoch |
923 |
> pred <-
|
| 61169 |
ripley |
924 |
+ predict(arima(window(NileNA, 1871, 1930), c(2, 0, 0)), n.ahead = 20)
|
| 56186 |
murdoch |
925 |
> lines(pred$pred, lty = 3, col = "red")
|
| 61157 |
ripley |
926 |
> lines(pred$pred + 2*pred$se, lty = 2, col = "blue")
|
|
|
927 |
> lines(pred$pred - 2*pred$se, lty = 2, col = "blue")
|
| 56186 |
murdoch |
928 |
>
|
|
|
929 |
> ## Structural time series models
|
|
|
930 |
> par(mfrow = c(3, 1))
|
|
|
931 |
> plot(Nile)
|
|
|
932 |
> ## local level model
|
|
|
933 |
> (fit <- StructTS(Nile, type = "level"))
|
|
|
934 |
|
|
|
935 |
Call:
|
|
|
936 |
StructTS(x = Nile, type = "level")
|
|
|
937 |
|
|
|
938 |
Variances:
|
|
|
939 |
level epsilon
|
|
|
940 |
1469 15099
|
|
|
941 |
> lines(fitted(fit), lty = 2) # contemporaneous smoothing
|
|
|
942 |
> lines(tsSmooth(fit), lty = 2, col = 4) # fixed-interval smoothing
|
|
|
943 |
> plot(residuals(fit)); abline(h = 0, lty = 3)
|
|
|
944 |
> ## local trend model
|
|
|
945 |
> (fit2 <- StructTS(Nile, type = "trend")) ## constant trend fitted
|
|
|
946 |
|
|
|
947 |
Call:
|
|
|
948 |
StructTS(x = Nile, type = "trend")
|
|
|
949 |
|
|
|
950 |
Variances:
|
|
|
951 |
level slope epsilon
|
|
|
952 |
1427 0 15047
|
|
|
953 |
> pred <- predict(fit, n.ahead = 30)
|
|
|
954 |
> ## with 50% confidence interval
|
|
|
955 |
> ts.plot(Nile, pred$pred,
|
|
|
956 |
+ pred$pred + 0.67*pred$se, pred$pred -0.67*pred$se)
|
|
|
957 |
>
|
|
|
958 |
> ## Now consider missing values
|
|
|
959 |
> plot(NileNA)
|
|
|
960 |
> (fit3 <- StructTS(NileNA, type = "level"))
|
|
|
961 |
|
|
|
962 |
Call:
|
|
|
963 |
StructTS(x = NileNA, type = "level")
|
|
|
964 |
|
|
|
965 |
Variances:
|
|
|
966 |
level epsilon
|
|
|
967 |
685.8 17899.8
|
|
|
968 |
> lines(fitted(fit3), lty = 2)
|
|
|
969 |
> lines(tsSmooth(fit3), lty = 3)
|
|
|
970 |
> plot(residuals(fit3)); abline(h = 0, lty = 3)
|
|
|
971 |
>
|
|
|
972 |
>
|
|
|
973 |
>
|
|
|
974 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
975 |
> cleanEx()
|
|
|
976 |
> nameEx("Orange")
|
|
|
977 |
> ### * Orange
|
|
|
978 |
>
|
|
|
979 |
> flush(stderr()); flush(stdout())
|
|
|
980 |
>
|
|
|
981 |
> ### Name: Orange
|
|
|
982 |
> ### Title: Growth of Orange Trees
|
|
|
983 |
> ### Aliases: Orange
|
|
|
984 |
> ### Keywords: datasets
|
|
|
985 |
>
|
|
|
986 |
> ### ** Examples
|
|
|
987 |
>
|
|
|
988 |
> require(stats); require(graphics)
|
|
|
989 |
> coplot(circumference ~ age | Tree, data = Orange, show.given = FALSE)
|
|
|
990 |
> fm1 <- nls(circumference ~ SSlogis(age, Asym, xmid, scal),
|
|
|
991 |
+ data = Orange, subset = Tree == 3)
|
|
|
992 |
> plot(circumference ~ age, data = Orange, subset = Tree == 3,
|
|
|
993 |
+ xlab = "Tree age (days since 1968/12/31)",
|
|
|
994 |
+ ylab = "Tree circumference (mm)", las = 1,
|
|
|
995 |
+ main = "Orange tree data and fitted model (Tree 3 only)")
|
|
|
996 |
> age <- seq(0, 1600, length.out = 101)
|
|
|
997 |
> lines(age, predict(fm1, list(age = age)))
|
|
|
998 |
>
|
|
|
999 |
>
|
|
|
1000 |
>
|
|
|
1001 |
> cleanEx()
|
|
|
1002 |
> nameEx("OrchardSprays")
|
|
|
1003 |
> ### * OrchardSprays
|
|
|
1004 |
>
|
|
|
1005 |
> flush(stderr()); flush(stdout())
|
|
|
1006 |
>
|
|
|
1007 |
> ### Name: OrchardSprays
|
|
|
1008 |
> ### Title: Potency of Orchard Sprays
|
|
|
1009 |
> ### Aliases: OrchardSprays
|
|
|
1010 |
> ### Keywords: datasets
|
|
|
1011 |
>
|
|
|
1012 |
> ### ** Examples
|
|
|
1013 |
>
|
|
|
1014 |
> require(graphics)
|
|
|
1015 |
> pairs(OrchardSprays, main = "OrchardSprays data")
|
|
|
1016 |
>
|
|
|
1017 |
>
|
|
|
1018 |
>
|
|
|
1019 |
> cleanEx()
|
|
|
1020 |
> nameEx("PlantGrowth")
|
|
|
1021 |
> ### * PlantGrowth
|
|
|
1022 |
>
|
|
|
1023 |
> flush(stderr()); flush(stdout())
|
|
|
1024 |
>
|
|
|
1025 |
> ### Name: PlantGrowth
|
|
|
1026 |
> ### Title: Results from an Experiment on Plant Growth
|
|
|
1027 |
> ### Aliases: PlantGrowth
|
|
|
1028 |
> ### Keywords: datasets
|
|
|
1029 |
>
|
|
|
1030 |
> ### ** Examples
|
|
|
1031 |
>
|
|
|
1032 |
> ## One factor ANOVA example from Dobson's book, cf. Table 7.4:
|
|
|
1033 |
> require(stats); require(graphics)
|
|
|
1034 |
> boxplot(weight ~ group, data = PlantGrowth, main = "PlantGrowth data",
|
|
|
1035 |
+ ylab = "Dried weight of plants", col = "lightgray",
|
|
|
1036 |
+ notch = TRUE, varwidth = TRUE)
|
|
|
1037 |
Warning in bxp(list(stats = c(4.17, 4.53, 5.155, 5.33, 6.11, 3.59, 4.17, :
|
|
|
1038 |
some notches went outside hinges ('box'): maybe set notch=FALSE
|
|
|
1039 |
> anova(lm(weight ~ group, data = PlantGrowth))
|
|
|
1040 |
Analysis of Variance Table
|
|
|
1041 |
|
|
|
1042 |
Response: weight
|
|
|
1043 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1044 |
group 2 3.7663 1.8832 4.8461 0.01591 *
|
|
|
1045 |
Residuals 27 10.4921 0.3886
|
|
|
1046 |
---
|
| 61435 |
ripley |
1047 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1048 |
>
|
|
|
1049 |
>
|
|
|
1050 |
>
|
|
|
1051 |
> cleanEx()
|
|
|
1052 |
> nameEx("Puromycin")
|
|
|
1053 |
> ### * Puromycin
|
|
|
1054 |
>
|
|
|
1055 |
> flush(stderr()); flush(stdout())
|
|
|
1056 |
>
|
|
|
1057 |
> ### Name: Puromycin
|
|
|
1058 |
> ### Title: Reaction Velocity of an Enzymatic Reaction
|
|
|
1059 |
> ### Aliases: Puromycin
|
|
|
1060 |
> ### Keywords: datasets
|
|
|
1061 |
>
|
|
|
1062 |
> ### ** Examples
|
|
|
1063 |
>
|
|
|
1064 |
> require(stats); require(graphics)
|
|
|
1065 |
> ## Don't show:
|
|
|
1066 |
> options(show.nls.convergence=FALSE)
|
|
|
1067 |
> ## End Don't show
|
|
|
1068 |
> plot(rate ~ conc, data = Puromycin, las = 1,
|
|
|
1069 |
+ xlab = "Substrate concentration (ppm)",
|
|
|
1070 |
+ ylab = "Reaction velocity (counts/min/min)",
|
|
|
1071 |
+ pch = as.integer(Puromycin$state),
|
|
|
1072 |
+ col = as.integer(Puromycin$state),
|
|
|
1073 |
+ main = "Puromycin data and fitted Michaelis-Menten curves")
|
|
|
1074 |
> ## simplest form of fitting the Michaelis-Menten model to these data
|
|
|
1075 |
> fm1 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
|
|
|
1076 |
+ subset = state == "treated",
|
|
|
1077 |
+ start = c(Vm = 200, K = 0.05))
|
|
|
1078 |
> fm2 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
|
|
|
1079 |
+ subset = state == "untreated",
|
|
|
1080 |
+ start = c(Vm = 160, K = 0.05))
|
|
|
1081 |
> summary(fm1)
|
|
|
1082 |
|
|
|
1083 |
Formula: rate ~ Vm * conc/(K + conc)
|
|
|
1084 |
|
|
|
1085 |
Parameters:
|
|
|
1086 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1087 |
Vm 2.127e+02 6.947e+00 30.615 3.24e-11 ***
|
|
|
1088 |
K 6.412e-02 8.281e-03 7.743 1.57e-05 ***
|
|
|
1089 |
---
|
| 61435 |
ripley |
1090 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1091 |
|
| 62439 |
ripley |
1092 |
Residual standard error: 10.93 on 10 degrees of freedom
|
| 56186 |
murdoch |
1093 |
|
|
|
1094 |
> summary(fm2)
|
|
|
1095 |
|
|
|
1096 |
Formula: rate ~ Vm * conc/(K + conc)
|
|
|
1097 |
|
|
|
1098 |
Parameters:
|
|
|
1099 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1100 |
Vm 1.603e+02 6.480e+00 24.734 1.38e-09 ***
|
|
|
1101 |
K 4.771e-02 7.782e-03 6.131 0.000173 ***
|
|
|
1102 |
---
|
| 61435 |
ripley |
1103 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1104 |
|
| 62439 |
ripley |
1105 |
Residual standard error: 9.773 on 9 degrees of freedom
|
| 56186 |
murdoch |
1106 |
|
|
|
1107 |
> ## add fitted lines to the plot
|
|
|
1108 |
> conc <- seq(0, 1.2, length.out = 101)
|
|
|
1109 |
> lines(conc, predict(fm1, list(conc = conc)), lty = 1, col = 1)
|
|
|
1110 |
> lines(conc, predict(fm2, list(conc = conc)), lty = 2, col = 2)
|
|
|
1111 |
> legend(0.8, 120, levels(Puromycin$state),
|
|
|
1112 |
+ col = 1:2, lty = 1:2, pch = 1:2)
|
|
|
1113 |
>
|
|
|
1114 |
> ## using partial linearity
|
|
|
1115 |
> fm3 <- nls(rate ~ conc/(K + conc), data = Puromycin,
|
|
|
1116 |
+ subset = state == "treated", start = c(K = 0.05),
|
|
|
1117 |
+ algorithm = "plinear")
|
|
|
1118 |
>
|
|
|
1119 |
>
|
|
|
1120 |
>
|
|
|
1121 |
> cleanEx()
|
|
|
1122 |
> nameEx("Theoph")
|
|
|
1123 |
> ### * Theoph
|
|
|
1124 |
>
|
|
|
1125 |
> flush(stderr()); flush(stdout())
|
|
|
1126 |
>
|
|
|
1127 |
> ### Name: Theoph
|
|
|
1128 |
> ### Title: Pharmacokinetics of Theophylline
|
|
|
1129 |
> ### Aliases: Theoph
|
|
|
1130 |
> ### Keywords: datasets
|
|
|
1131 |
>
|
|
|
1132 |
> ### ** Examples
|
|
|
1133 |
>
|
|
|
1134 |
> require(stats); require(graphics)
|
|
|
1135 |
> ## Don't show:
|
|
|
1136 |
> options(show.nls.convergence=FALSE)
|
|
|
1137 |
> ## End Don't show
|
|
|
1138 |
> coplot(conc ~ Time | Subject, data = Theoph, show.given = FALSE)
|
|
|
1139 |
> Theoph.4 <- subset(Theoph, Subject == 4)
|
|
|
1140 |
> fm1 <- nls(conc ~ SSfol(Dose, Time, lKe, lKa, lCl),
|
|
|
1141 |
+ data = Theoph.4)
|
|
|
1142 |
> summary(fm1)
|
|
|
1143 |
|
|
|
1144 |
Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl)
|
|
|
1145 |
|
|
|
1146 |
Parameters:
|
|
|
1147 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1148 |
lKe -2.4365 0.2257 -10.797 4.77e-06 ***
|
|
|
1149 |
lKa 0.1583 0.2297 0.689 0.51
|
|
|
1150 |
lCl -3.2861 0.1448 -22.695 1.51e-08 ***
|
|
|
1151 |
---
|
| 61435 |
ripley |
1152 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1153 |
|
| 62439 |
ripley |
1154 |
Residual standard error: 0.8465 on 8 degrees of freedom
|
| 56186 |
murdoch |
1155 |
|
|
|
1156 |
> plot(conc ~ Time, data = Theoph.4,
|
|
|
1157 |
+ xlab = "Time since drug administration (hr)",
|
|
|
1158 |
+ ylab = "Theophylline concentration (mg/L)",
|
|
|
1159 |
+ main = "Observed concentrations and fitted model",
|
|
|
1160 |
+ sub = "Theophylline data - Subject 4 only",
|
|
|
1161 |
+ las = 1, col = 4)
|
|
|
1162 |
> xvals <- seq(0, par("usr")[2], length.out = 55)
|
|
|
1163 |
> lines(xvals, predict(fm1, newdata = list(Time = xvals)),
|
|
|
1164 |
+ col = 4)
|
|
|
1165 |
>
|
|
|
1166 |
>
|
|
|
1167 |
>
|
|
|
1168 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1169 |
> cleanEx()
|
|
|
1170 |
> nameEx("Titanic")
|
|
|
1171 |
> ### * Titanic
|
|
|
1172 |
>
|
|
|
1173 |
> flush(stderr()); flush(stdout())
|
|
|
1174 |
>
|
|
|
1175 |
> ### Name: Titanic
|
|
|
1176 |
> ### Title: Survival of passengers on the Titanic
|
|
|
1177 |
> ### Aliases: Titanic
|
|
|
1178 |
> ### Keywords: datasets
|
|
|
1179 |
>
|
|
|
1180 |
> ### ** Examples
|
|
|
1181 |
>
|
|
|
1182 |
> require(graphics)
|
|
|
1183 |
> mosaicplot(Titanic, main = "Survival on the Titanic")
|
|
|
1184 |
> ## Higher survival rates in children?
|
|
|
1185 |
> apply(Titanic, c(3, 4), sum)
|
|
|
1186 |
Survived
|
|
|
1187 |
Age No Yes
|
|
|
1188 |
Child 52 57
|
|
|
1189 |
Adult 1438 654
|
|
|
1190 |
> ## Higher survival rates in females?
|
|
|
1191 |
> apply(Titanic, c(2, 4), sum)
|
|
|
1192 |
Survived
|
|
|
1193 |
Sex No Yes
|
|
|
1194 |
Male 1364 367
|
|
|
1195 |
Female 126 344
|
|
|
1196 |
> ## Use loglm() in package 'MASS' for further analysis ...
|
|
|
1197 |
>
|
|
|
1198 |
>
|
|
|
1199 |
>
|
|
|
1200 |
> cleanEx()
|
|
|
1201 |
> nameEx("ToothGrowth")
|
|
|
1202 |
> ### * ToothGrowth
|
|
|
1203 |
>
|
|
|
1204 |
> flush(stderr()); flush(stdout())
|
|
|
1205 |
>
|
|
|
1206 |
> ### Name: ToothGrowth
|
|
|
1207 |
> ### Title: The Effect of Vitamin C on Tooth Growth in Guinea Pigs
|
|
|
1208 |
> ### Aliases: ToothGrowth
|
|
|
1209 |
> ### Keywords: datasets
|
|
|
1210 |
>
|
|
|
1211 |
> ### ** Examples
|
|
|
1212 |
>
|
|
|
1213 |
> require(graphics)
|
|
|
1214 |
> coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
|
|
|
1215 |
+ xlab = "ToothGrowth data: length vs dose, given type of supplement")
|
|
|
1216 |
>
|
|
|
1217 |
>
|
|
|
1218 |
>
|
|
|
1219 |
> cleanEx()
|
|
|
1220 |
> nameEx("UCBAdmissions")
|
|
|
1221 |
> ### * UCBAdmissions
|
|
|
1222 |
>
|
|
|
1223 |
> flush(stderr()); flush(stdout())
|
|
|
1224 |
>
|
|
|
1225 |
> ### Name: UCBAdmissions
|
|
|
1226 |
> ### Title: Student Admissions at UC Berkeley
|
|
|
1227 |
> ### Aliases: UCBAdmissions
|
|
|
1228 |
> ### Keywords: datasets
|
|
|
1229 |
>
|
|
|
1230 |
> ### ** Examples
|
|
|
1231 |
>
|
|
|
1232 |
> require(graphics)
|
|
|
1233 |
> ## Data aggregated over departments
|
|
|
1234 |
> apply(UCBAdmissions, c(1, 2), sum)
|
|
|
1235 |
Gender
|
|
|
1236 |
Admit Male Female
|
|
|
1237 |
Admitted 1198 557
|
|
|
1238 |
Rejected 1493 1278
|
|
|
1239 |
> mosaicplot(apply(UCBAdmissions, c(1, 2), sum),
|
|
|
1240 |
+ main = "Student admissions at UC Berkeley")
|
|
|
1241 |
> ## Data for individual departments
|
|
|
1242 |
> opar <- par(mfrow = c(2, 3), oma = c(0, 0, 2, 0))
|
|
|
1243 |
> for(i in 1:6)
|
|
|
1244 |
+ mosaicplot(UCBAdmissions[,,i],
|
|
|
1245 |
+ xlab = "Admit", ylab = "Sex",
|
|
|
1246 |
+ main = paste("Department", LETTERS[i]))
|
|
|
1247 |
> mtext(expression(bold("Student admissions at UC Berkeley")),
|
|
|
1248 |
+ outer = TRUE, cex = 1.5)
|
|
|
1249 |
> par(opar)
|
|
|
1250 |
>
|
|
|
1251 |
>
|
|
|
1252 |
>
|
|
|
1253 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1254 |
> cleanEx()
|
|
|
1255 |
> nameEx("UKDriverDeaths")
|
|
|
1256 |
> ### * UKDriverDeaths
|
|
|
1257 |
>
|
|
|
1258 |
> flush(stderr()); flush(stdout())
|
|
|
1259 |
>
|
|
|
1260 |
> ### Name: UKDriverDeaths
|
|
|
1261 |
> ### Title: Road Casualties in Great Britain 1969-84
|
|
|
1262 |
> ### Aliases: UKDriverDeaths Seatbelts
|
|
|
1263 |
> ### Keywords: datasets
|
|
|
1264 |
>
|
|
|
1265 |
> ### ** Examples
|
|
|
1266 |
>
|
|
|
1267 |
> require(stats); require(graphics)
|
|
|
1268 |
> ## work with pre-seatbelt period to identify a model, use logs
|
|
|
1269 |
> work <- window(log10(UKDriverDeaths), end = 1982+11/12)
|
| 61169 |
ripley |
1270 |
> par(mfrow = c(3, 1))
|
| 56186 |
murdoch |
1271 |
> plot(work); acf(work); pacf(work)
|
| 61169 |
ripley |
1272 |
> par(mfrow = c(1, 1))
|
|
|
1273 |
> (fit <- arima(work, c(1, 0, 0), seasonal = list(order = c(1, 0, 0))))
|
| 56186 |
murdoch |
1274 |
|
|
|
1275 |
Call:
|
|
|
1276 |
arima(x = work, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0)))
|
|
|
1277 |
|
|
|
1278 |
Coefficients:
|
|
|
1279 |
ar1 sar1 intercept
|
|
|
1280 |
0.4378 0.6281 3.2274
|
|
|
1281 |
s.e. 0.0764 0.0637 0.0131
|
|
|
1282 |
|
|
|
1283 |
sigma^2 estimated as 0.00157: log likelihood = 300.85, aic = -593.7
|
|
|
1284 |
> z <- predict(fit, n.ahead = 24)
|
|
|
1285 |
> ts.plot(log10(UKDriverDeaths), z$pred, z$pred+2*z$se, z$pred-2*z$se,
|
| 61169 |
ripley |
1286 |
+ lty = c(1, 3, 2, 2), col = c("black", "red", "blue", "blue"))
|
| 56186 |
murdoch |
1287 |
>
|
|
|
1288 |
> ## now see the effect of the explanatory variables
|
|
|
1289 |
> X <- Seatbelts[, c("kms", "PetrolPrice", "law")]
|
|
|
1290 |
> X[, 1] <- log10(X[, 1]) - 4
|
| 61169 |
ripley |
1291 |
> arima(log10(Seatbelts[, "drivers"]), c(1, 0, 0),
|
|
|
1292 |
+ seasonal = list(order = c(1, 0, 0)), xreg = X)
|
| 56186 |
murdoch |
1293 |
|
|
|
1294 |
Call:
|
|
|
1295 |
arima(x = log10(Seatbelts[, "drivers"]), order = c(1, 0, 0), seasonal = list(order = c(1,
|
|
|
1296 |
0, 0)), xreg = X)
|
|
|
1297 |
|
|
|
1298 |
Coefficients:
|
|
|
1299 |
ar1 sar1 intercept kms PetrolPrice law
|
|
|
1300 |
0.3348 0.6672 3.3539 0.0082 -1.2224 -0.0963
|
|
|
1301 |
s.e. 0.0775 0.0612 0.0441 0.0902 0.3839 0.0166
|
|
|
1302 |
|
|
|
1303 |
sigma^2 estimated as 0.001476: log likelihood = 349.73, aic = -685.46
|
|
|
1304 |
>
|
|
|
1305 |
>
|
|
|
1306 |
>
|
|
|
1307 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1308 |
> cleanEx()
|
|
|
1309 |
> nameEx("UKLungDeaths")
|
|
|
1310 |
> ### * UKLungDeaths
|
|
|
1311 |
>
|
|
|
1312 |
> flush(stderr()); flush(stdout())
|
|
|
1313 |
>
|
|
|
1314 |
> ### Name: UKLungDeaths
|
|
|
1315 |
> ### Title: Monthly Deaths from Lung Diseases in the UK
|
|
|
1316 |
> ### Aliases: UKLungDeaths ldeaths fdeaths mdeaths
|
|
|
1317 |
> ### Keywords: datasets
|
|
|
1318 |
>
|
|
|
1319 |
> ### ** Examples
|
|
|
1320 |
>
|
|
|
1321 |
> require(stats); require(graphics) # for time
|
| 61435 |
ripley |
1322 |
> plot(ldeaths)
|
|
|
1323 |
> plot(mdeaths, fdeaths)
|
| 56186 |
murdoch |
1324 |
> ## Better labels:
|
|
|
1325 |
> yr <- floor(tt <- time(mdeaths))
|
|
|
1326 |
> plot(mdeaths, fdeaths,
|
| 61157 |
ripley |
1327 |
+ xy.labels = paste(month.abb[12*(tt - yr)], yr-1900, sep = "'"))
|
| 56186 |
murdoch |
1328 |
>
|
|
|
1329 |
>
|
|
|
1330 |
>
|
|
|
1331 |
> cleanEx()
|
|
|
1332 |
> nameEx("UKgas")
|
|
|
1333 |
> ### * UKgas
|
|
|
1334 |
>
|
|
|
1335 |
> flush(stderr()); flush(stdout())
|
|
|
1336 |
>
|
|
|
1337 |
> ### Name: UKgas
|
|
|
1338 |
> ### Title: UK Quarterly Gas Consumption
|
|
|
1339 |
> ### Aliases: UKgas
|
|
|
1340 |
> ### Keywords: datasets
|
|
|
1341 |
>
|
|
|
1342 |
> ### ** Examples
|
|
|
1343 |
>
|
|
|
1344 |
> ## maybe str(UKgas) ; plot(UKgas) ...
|
|
|
1345 |
>
|
|
|
1346 |
>
|
|
|
1347 |
>
|
|
|
1348 |
> cleanEx()
|
|
|
1349 |
> nameEx("USArrests")
|
|
|
1350 |
> ### * USArrests
|
|
|
1351 |
>
|
|
|
1352 |
> flush(stderr()); flush(stdout())
|
|
|
1353 |
>
|
|
|
1354 |
> ### Name: USArrests
|
|
|
1355 |
> ### Title: Violent Crime Rates by US State
|
|
|
1356 |
> ### Aliases: USArrests
|
|
|
1357 |
> ### Keywords: datasets
|
|
|
1358 |
>
|
|
|
1359 |
> ### ** Examples
|
|
|
1360 |
>
|
|
|
1361 |
> require(graphics)
|
|
|
1362 |
> pairs(USArrests, panel = panel.smooth, main = "USArrests data")
|
|
|
1363 |
>
|
|
|
1364 |
>
|
|
|
1365 |
>
|
|
|
1366 |
> cleanEx()
|
|
|
1367 |
> nameEx("USJudgeRatings")
|
|
|
1368 |
> ### * USJudgeRatings
|
|
|
1369 |
>
|
|
|
1370 |
> flush(stderr()); flush(stdout())
|
|
|
1371 |
>
|
|
|
1372 |
> ### Name: USJudgeRatings
|
|
|
1373 |
> ### Title: Lawyers' Ratings of State Judges in the US Superior Court
|
|
|
1374 |
> ### Aliases: USJudgeRatings
|
|
|
1375 |
> ### Keywords: datasets
|
|
|
1376 |
>
|
|
|
1377 |
> ### ** Examples
|
|
|
1378 |
>
|
|
|
1379 |
> require(graphics)
|
|
|
1380 |
> pairs(USJudgeRatings, main = "USJudgeRatings data")
|
|
|
1381 |
>
|
|
|
1382 |
>
|
|
|
1383 |
>
|
|
|
1384 |
> cleanEx()
|
|
|
1385 |
> nameEx("USPersonalExpenditure")
|
|
|
1386 |
> ### * USPersonalExpenditure
|
|
|
1387 |
>
|
|
|
1388 |
> flush(stderr()); flush(stdout())
|
|
|
1389 |
>
|
|
|
1390 |
> ### Name: USPersonalExpenditure
|
|
|
1391 |
> ### Title: Personal Expenditure Data
|
|
|
1392 |
> ### Aliases: USPersonalExpenditure
|
|
|
1393 |
> ### Keywords: datasets
|
|
|
1394 |
>
|
|
|
1395 |
> ### ** Examples
|
|
|
1396 |
>
|
|
|
1397 |
> require(stats) # for medpolish
|
|
|
1398 |
> USPersonalExpenditure
|
|
|
1399 |
1940 1945 1950 1955 1960
|
|
|
1400 |
Food and Tobacco 22.200 44.500 59.60 73.2 86.80
|
|
|
1401 |
Household Operation 10.500 15.500 29.00 36.5 46.20
|
|
|
1402 |
Medical and Health 3.530 5.760 9.71 14.0 21.10
|
|
|
1403 |
Personal Care 1.040 1.980 2.45 3.4 5.40
|
|
|
1404 |
Private Education 0.341 0.974 1.80 2.6 3.64
|
|
|
1405 |
> medpolish(log10(USPersonalExpenditure))
|
| 61435 |
ripley |
1406 |
1: 1.126317
|
|
|
1407 |
2: 1.032421
|
|
|
1408 |
Final: 1.032421
|
| 56186 |
murdoch |
1409 |
|
|
|
1410 |
Median Polish Results (Dataset: "log10(USPersonalExpenditure)")
|
|
|
1411 |
|
| 61435 |
ripley |
1412 |
Overall: 0.9872192
|
| 56186 |
murdoch |
1413 |
|
|
|
1414 |
Row Effects:
|
|
|
1415 |
Food and Tobacco Household Operation Medical and Health Personal Care
|
|
|
1416 |
0.7880270 0.4327608 0.0000000 -0.5606543
|
|
|
1417 |
Private Education
|
|
|
1418 |
-0.7319467
|
|
|
1419 |
|
|
|
1420 |
Column Effects:
|
|
|
1421 |
1940 1945 1950 1955 1960
|
|
|
1422 |
-0.4288933 -0.2267967 0.0000000 0.1423128 0.3058289
|
|
|
1423 |
|
|
|
1424 |
Residuals:
|
|
|
1425 |
1940 1945 1950 1955 1960
|
|
|
1426 |
Food and Tobacco 0.000000 0.0999105 0.000000 -0.053048 -0.142555
|
|
|
1427 |
Household Operation 0.030103 -0.0028516 0.042418 0.000000 -0.061167
|
|
|
1428 |
Medical and Health -0.010551 0.0000000 0.000000 0.016596 0.031234
|
|
|
1429 |
Personal Care 0.019362 0.0968971 -0.037399 -0.037399 0.000000
|
|
|
1430 |
Private Education -0.293625 -0.0399168 0.000000 0.017388 0.000000
|
|
|
1431 |
|
|
|
1432 |
>
|
|
|
1433 |
>
|
|
|
1434 |
>
|
|
|
1435 |
> cleanEx()
|
|
|
1436 |
> nameEx("VADeaths")
|
|
|
1437 |
> ### * VADeaths
|
|
|
1438 |
>
|
|
|
1439 |
> flush(stderr()); flush(stdout())
|
|
|
1440 |
>
|
|
|
1441 |
> ### Name: VADeaths
|
|
|
1442 |
> ### Title: Death Rates in Virginia (1940)
|
|
|
1443 |
> ### Aliases: VADeaths
|
|
|
1444 |
> ### Keywords: datasets
|
|
|
1445 |
>
|
|
|
1446 |
> ### ** Examples
|
|
|
1447 |
>
|
|
|
1448 |
> require(stats); require(graphics)
|
|
|
1449 |
> n <- length(dr <- c(VADeaths))
|
|
|
1450 |
> nam <- names(VADeaths)
|
|
|
1451 |
> d.VAD <- data.frame(
|
|
|
1452 |
+ Drate = dr,
|
| 61157 |
ripley |
1453 |
+ age = rep(ordered(rownames(VADeaths)), length.out = n),
|
| 61169 |
ripley |
1454 |
+ gender = gl(2, 5, n, labels = c("M", "F")),
|
|
|
1455 |
+ site = gl(2, 10, labels = c("rural", "urban")))
|
| 56186 |
murdoch |
1456 |
> coplot(Drate ~ as.numeric(age) | gender * site, data = d.VAD,
|
|
|
1457 |
+ panel = panel.smooth, xlab = "VADeaths data - Given: gender")
|
|
|
1458 |
> summary(aov.VAD <- aov(Drate ~ .^2, data = d.VAD))
|
| 57044 |
ripley |
1459 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1460 |
age 4 6288 1572.1 590.858 8.55e-06 ***
|
|
|
1461 |
gender 1 648 647.5 243.361 9.86e-05 ***
|
|
|
1462 |
site 1 77 76.8 28.876 0.00579 **
|
|
|
1463 |
age:gender 4 86 21.6 8.100 0.03358 *
|
|
|
1464 |
age:site 4 43 10.6 3.996 0.10414
|
|
|
1465 |
gender:site 1 73 73.0 27.422 0.00636 **
|
|
|
1466 |
Residuals 4 11 2.7
|
| 56186 |
murdoch |
1467 |
---
|
| 61435 |
ripley |
1468 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 61169 |
ripley |
1469 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
1470 |
> plot(aov.VAD)
|
|
|
1471 |
> par(opar)
|
|
|
1472 |
>
|
|
|
1473 |
>
|
|
|
1474 |
>
|
|
|
1475 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1476 |
> cleanEx()
|
|
|
1477 |
> nameEx("WWWusage")
|
|
|
1478 |
> ### * WWWusage
|
|
|
1479 |
>
|
|
|
1480 |
> flush(stderr()); flush(stdout())
|
|
|
1481 |
>
|
|
|
1482 |
> ### Name: WWWusage
|
|
|
1483 |
> ### Title: Internet Usage per Minute
|
|
|
1484 |
> ### Aliases: WWWusage
|
|
|
1485 |
> ### Keywords: datasets
|
|
|
1486 |
>
|
|
|
1487 |
> ### ** Examples
|
|
|
1488 |
>
|
|
|
1489 |
> require(graphics)
|
|
|
1490 |
> work <- diff(WWWusage)
|
| 61169 |
ripley |
1491 |
> par(mfrow = c(2, 1)); plot(WWWusage); plot(work)
|
| 56186 |
murdoch |
1492 |
> ## Not run:
|
|
|
1493 |
> ##D require(stats)
|
| 61157 |
ripley |
1494 |
> ##D aics <- matrix(, 6, 6, dimnames = list(p = 0:5, q = 0:5))
|
| 61169 |
ripley |
1495 |
> ##D for(q in 1:5) aics[1, 1+q] <- arima(WWWusage, c(0, 1, q),
|
| 56186 |
murdoch |
1496 |
> ##D optim.control = list(maxit = 500))$aic
|
|
|
1497 |
> ##D for(p in 1:5)
|
| 61169 |
ripley |
1498 |
> ##D for(q in 0:5) aics[1+p, 1+q] <- arima(WWWusage, c(p, 1, q),
|
| 56186 |
murdoch |
1499 |
> ##D optim.control = list(maxit = 500))$aic
|
| 61157 |
ripley |
1500 |
> ##D round(aics - min(aics, na.rm = TRUE), 2)
|
| 56186 |
murdoch |
1501 |
> ## End(Not run)
|
|
|
1502 |
>
|
|
|
1503 |
>
|
|
|
1504 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1505 |
> cleanEx()
|
|
|
1506 |
> nameEx("WorldPhones")
|
|
|
1507 |
> ### * WorldPhones
|
|
|
1508 |
>
|
|
|
1509 |
> flush(stderr()); flush(stdout())
|
|
|
1510 |
>
|
|
|
1511 |
> ### Name: WorldPhones
|
|
|
1512 |
> ### Title: The World's Telephones
|
|
|
1513 |
> ### Aliases: WorldPhones
|
|
|
1514 |
> ### Keywords: datasets
|
|
|
1515 |
>
|
|
|
1516 |
> ### ** Examples
|
|
|
1517 |
>
|
|
|
1518 |
> require(graphics)
|
|
|
1519 |
> matplot(rownames(WorldPhones), WorldPhones, type = "b", log = "y",
|
|
|
1520 |
+ xlab = "Year", ylab = "Number of telephones (1000's)")
|
| 61435 |
ripley |
1521 |
> legend(1951.5, 80000, colnames(WorldPhones), col = 1:6, lty = 1:5,
|
| 56186 |
murdoch |
1522 |
+ pch = rep(21, 7))
|
|
|
1523 |
> title(main = "World phones data: log scale for response")
|
|
|
1524 |
>
|
|
|
1525 |
>
|
|
|
1526 |
>
|
|
|
1527 |
> cleanEx()
|
|
|
1528 |
> nameEx("ability.cov")
|
|
|
1529 |
> ### * ability.cov
|
|
|
1530 |
>
|
|
|
1531 |
> flush(stderr()); flush(stdout())
|
|
|
1532 |
>
|
|
|
1533 |
> ### Name: ability.cov
|
|
|
1534 |
> ### Title: Ability and Intelligence Tests
|
|
|
1535 |
> ### Aliases: ability.cov
|
|
|
1536 |
> ### Keywords: datasets
|
|
|
1537 |
>
|
|
|
1538 |
> ### ** Examples
|
|
|
1539 |
>
|
|
|
1540 |
>
|
|
|
1541 |
> cleanEx()
|
|
|
1542 |
> nameEx("airmiles")
|
|
|
1543 |
> ### * airmiles
|
|
|
1544 |
>
|
|
|
1545 |
> flush(stderr()); flush(stdout())
|
|
|
1546 |
>
|
|
|
1547 |
> ### Name: airmiles
|
|
|
1548 |
> ### Title: Passenger Miles on Commercial US Airlines, 1937-1960
|
|
|
1549 |
> ### Aliases: airmiles
|
|
|
1550 |
> ### Keywords: datasets
|
|
|
1551 |
>
|
|
|
1552 |
> ### ** Examples
|
|
|
1553 |
>
|
|
|
1554 |
> require(graphics)
|
|
|
1555 |
> plot(airmiles, main = "airmiles data",
|
|
|
1556 |
+ xlab = "Passenger-miles flown by U.S. commercial airlines", col = 4)
|
|
|
1557 |
>
|
|
|
1558 |
>
|
|
|
1559 |
>
|
|
|
1560 |
> cleanEx()
|
|
|
1561 |
> nameEx("airquality")
|
|
|
1562 |
> ### * airquality
|
|
|
1563 |
>
|
|
|
1564 |
> flush(stderr()); flush(stdout())
|
|
|
1565 |
>
|
|
|
1566 |
> ### Name: airquality
|
|
|
1567 |
> ### Title: New York Air Quality Measurements
|
|
|
1568 |
> ### Aliases: airquality
|
|
|
1569 |
> ### Keywords: datasets
|
|
|
1570 |
>
|
|
|
1571 |
> ### ** Examples
|
|
|
1572 |
>
|
|
|
1573 |
> require(graphics)
|
|
|
1574 |
> pairs(airquality, panel = panel.smooth, main = "airquality data")
|
|
|
1575 |
>
|
|
|
1576 |
>
|
|
|
1577 |
>
|
|
|
1578 |
> cleanEx()
|
|
|
1579 |
> nameEx("anscombe")
|
|
|
1580 |
> ### * anscombe
|
|
|
1581 |
>
|
|
|
1582 |
> flush(stderr()); flush(stdout())
|
|
|
1583 |
>
|
|
|
1584 |
> ### Name: anscombe
|
|
|
1585 |
> ### Title: Anscombe's Quartet of 'Identical' Simple Linear Regressions
|
|
|
1586 |
> ### Aliases: anscombe
|
|
|
1587 |
> ### Keywords: datasets
|
|
|
1588 |
>
|
|
|
1589 |
> ### ** Examples
|
|
|
1590 |
>
|
|
|
1591 |
> require(stats); require(graphics)
|
|
|
1592 |
> summary(anscombe)
|
|
|
1593 |
x1 x2 x3 x4 y1
|
|
|
1594 |
Min. : 4.0 Min. : 4.0 Min. : 4.0 Min. : 8 Min. : 4.260
|
|
|
1595 |
1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 6.5 1st Qu.: 8 1st Qu.: 6.315
|
|
|
1596 |
Median : 9.0 Median : 9.0 Median : 9.0 Median : 8 Median : 7.580
|
|
|
1597 |
Mean : 9.0 Mean : 9.0 Mean : 9.0 Mean : 9 Mean : 7.501
|
|
|
1598 |
3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.:11.5 3rd Qu.: 8 3rd Qu.: 8.570
|
|
|
1599 |
Max. :14.0 Max. :14.0 Max. :14.0 Max. :19 Max. :10.840
|
|
|
1600 |
y2 y3 y4
|
|
|
1601 |
Min. :3.100 Min. : 5.39 Min. : 5.250
|
|
|
1602 |
1st Qu.:6.695 1st Qu.: 6.25 1st Qu.: 6.170
|
|
|
1603 |
Median :8.140 Median : 7.11 Median : 7.040
|
|
|
1604 |
Mean :7.501 Mean : 7.50 Mean : 7.501
|
|
|
1605 |
3rd Qu.:8.950 3rd Qu.: 7.98 3rd Qu.: 8.190
|
|
|
1606 |
Max. :9.260 Max. :12.74 Max. :12.500
|
|
|
1607 |
>
|
|
|
1608 |
> ##-- now some "magic" to do the 4 regressions in a loop:
|
|
|
1609 |
> ff <- y ~ x
|
| 57821 |
maechler |
1610 |
> mods <- setNames(as.list(1:4), paste0("lm", 1:4))
|
| 56186 |
murdoch |
1611 |
> for(i in 1:4) {
|
| 57821 |
maechler |
1612 |
+ ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
|
|
|
1613 |
+ ## or ff[[2]] <- as.name(paste0("y", i))
|
|
|
1614 |
+ ## ff[[3]] <- as.name(paste0("x", i))
|
| 61157 |
ripley |
1615 |
+ mods[[i]] <- lmi <- lm(ff, data = anscombe)
|
| 56186 |
murdoch |
1616 |
+ print(anova(lmi))
|
|
|
1617 |
+ }
|
|
|
1618 |
Analysis of Variance Table
|
|
|
1619 |
|
|
|
1620 |
Response: y1
|
|
|
1621 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1622 |
x1 1 27.510 27.5100 17.99 0.00217 **
|
|
|
1623 |
Residuals 9 13.763 1.5292
|
|
|
1624 |
---
|
| 61435 |
ripley |
1625 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1626 |
Analysis of Variance Table
|
|
|
1627 |
|
|
|
1628 |
Response: y2
|
|
|
1629 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1630 |
x2 1 27.500 27.5000 17.966 0.002179 **
|
|
|
1631 |
Residuals 9 13.776 1.5307
|
|
|
1632 |
---
|
| 61435 |
ripley |
1633 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1634 |
Analysis of Variance Table
|
|
|
1635 |
|
|
|
1636 |
Response: y3
|
|
|
1637 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1638 |
x3 1 27.470 27.4700 17.972 0.002176 **
|
|
|
1639 |
Residuals 9 13.756 1.5285
|
|
|
1640 |
---
|
| 61435 |
ripley |
1641 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1642 |
Analysis of Variance Table
|
|
|
1643 |
|
|
|
1644 |
Response: y4
|
|
|
1645 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1646 |
x4 1 27.490 27.4900 18.003 0.002165 **
|
|
|
1647 |
Residuals 9 13.742 1.5269
|
|
|
1648 |
---
|
| 61435 |
ripley |
1649 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1650 |
>
|
|
|
1651 |
> ## See how close they are (numerically!)
|
| 57821 |
maechler |
1652 |
> sapply(mods, coef)
|
|
|
1653 |
lm1 lm2 lm3 lm4
|
| 56186 |
murdoch |
1654 |
(Intercept) 3.0000909 3.000909 3.0024545 3.0017273
|
|
|
1655 |
x1 0.5000909 0.500000 0.4997273 0.4999091
|
| 57821 |
maechler |
1656 |
> lapply(mods, function(fm) coef(summary(fm)))
|
|
|
1657 |
$lm1
|
| 56186 |
murdoch |
1658 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1659 |
(Intercept) 3.0000909 1.1247468 2.667348 0.025734051
|
|
|
1660 |
x1 0.5000909 0.1179055 4.241455 0.002169629
|
|
|
1661 |
|
| 57821 |
maechler |
1662 |
$lm2
|
| 56186 |
murdoch |
1663 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1664 |
(Intercept) 3.000909 1.1253024 2.666758 0.025758941
|
|
|
1665 |
x2 0.500000 0.1179637 4.238590 0.002178816
|
|
|
1666 |
|
| 57821 |
maechler |
1667 |
$lm3
|
| 56186 |
murdoch |
1668 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1669 |
(Intercept) 3.0024545 1.1244812 2.670080 0.025619109
|
|
|
1670 |
x3 0.4997273 0.1178777 4.239372 0.002176305
|
|
|
1671 |
|
| 57821 |
maechler |
1672 |
$lm4
|
| 56186 |
murdoch |
1673 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1674 |
(Intercept) 3.0017273 1.1239211 2.670763 0.025590425
|
|
|
1675 |
x4 0.4999091 0.1178189 4.243028 0.002164602
|
|
|
1676 |
|
|
|
1677 |
>
|
|
|
1678 |
> ## Now, do what you should have done in the first place: PLOTS
|
| 61169 |
ripley |
1679 |
> op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma = c(0, 0, 2, 0))
|
| 56186 |
murdoch |
1680 |
> for(i in 1:4) {
|
| 57821 |
maechler |
1681 |
+ ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
|
| 61157 |
ripley |
1682 |
+ plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
|
| 61169 |
ripley |
1683 |
+ xlim = c(3, 19), ylim = c(3, 13))
|
| 61157 |
ripley |
1684 |
+ abline(mods[[i]], col = "blue")
|
| 56186 |
murdoch |
1685 |
+ }
|
| 61157 |
ripley |
1686 |
> mtext("Anscombe's 4 Regression data sets", outer = TRUE, cex = 1.5)
|
| 56186 |
murdoch |
1687 |
> par(op)
|
|
|
1688 |
>
|
|
|
1689 |
>
|
|
|
1690 |
>
|
|
|
1691 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1692 |
> cleanEx()
|
|
|
1693 |
> nameEx("attenu")
|
|
|
1694 |
> ### * attenu
|
|
|
1695 |
>
|
|
|
1696 |
> flush(stderr()); flush(stdout())
|
|
|
1697 |
>
|
|
|
1698 |
> ### Name: attenu
|
|
|
1699 |
> ### Title: The Joyner-Boore Attenuation Data
|
|
|
1700 |
> ### Aliases: attenu
|
|
|
1701 |
> ### Keywords: datasets
|
|
|
1702 |
>
|
|
|
1703 |
> ### ** Examples
|
|
|
1704 |
>
|
|
|
1705 |
> require(graphics)
|
|
|
1706 |
> ## check the data class of the variables
|
|
|
1707 |
> sapply(attenu, data.class)
|
|
|
1708 |
event mag station dist accel
|
|
|
1709 |
"numeric" "numeric" "factor" "numeric" "numeric"
|
|
|
1710 |
> summary(attenu)
|
|
|
1711 |
event mag station dist
|
|
|
1712 |
Min. : 1.00 Min. :5.000 117 : 5 Min. : 0.50
|
|
|
1713 |
1st Qu.: 9.00 1st Qu.:5.300 1028 : 4 1st Qu.: 11.32
|
|
|
1714 |
Median :18.00 Median :6.100 113 : 4 Median : 23.40
|
|
|
1715 |
Mean :14.74 Mean :6.084 112 : 3 Mean : 45.60
|
|
|
1716 |
3rd Qu.:20.00 3rd Qu.:6.600 135 : 3 3rd Qu.: 47.55
|
|
|
1717 |
Max. :23.00 Max. :7.700 (Other):147 Max. :370.00
|
|
|
1718 |
NA's : 16
|
|
|
1719 |
accel
|
|
|
1720 |
Min. :0.00300
|
|
|
1721 |
1st Qu.:0.04425
|
|
|
1722 |
Median :0.11300
|
|
|
1723 |
Mean :0.15422
|
|
|
1724 |
3rd Qu.:0.21925
|
|
|
1725 |
Max. :0.81000
|
|
|
1726 |
|
|
|
1727 |
> pairs(attenu, main = "attenu data")
|
|
|
1728 |
> coplot(accel ~ dist | as.factor(event), data = attenu, show.given = FALSE)
|
|
|
1729 |
> coplot(log(accel) ~ log(dist) | as.factor(event),
|
|
|
1730 |
+ data = attenu, panel = panel.smooth, show.given = FALSE)
|
|
|
1731 |
>
|
|
|
1732 |
>
|
|
|
1733 |
>
|
|
|
1734 |
> cleanEx()
|
|
|
1735 |
> nameEx("attitude")
|
|
|
1736 |
> ### * attitude
|
|
|
1737 |
>
|
|
|
1738 |
> flush(stderr()); flush(stdout())
|
|
|
1739 |
>
|
|
|
1740 |
> ### Name: attitude
|
|
|
1741 |
> ### Title: The Chatterjee-Price Attitude Data
|
|
|
1742 |
> ### Aliases: attitude
|
|
|
1743 |
> ### Keywords: datasets
|
|
|
1744 |
>
|
|
|
1745 |
> ### ** Examples
|
|
|
1746 |
>
|
|
|
1747 |
> require(stats); require(graphics)
|
|
|
1748 |
> pairs(attitude, main = "attitude data")
|
|
|
1749 |
> summary(attitude)
|
|
|
1750 |
rating complaints privileges learning raises
|
|
|
1751 |
Min. :40.00 Min. :37.0 Min. :30.00 Min. :34.00 Min. :43.00
|
|
|
1752 |
1st Qu.:58.75 1st Qu.:58.5 1st Qu.:45.00 1st Qu.:47.00 1st Qu.:58.25
|
|
|
1753 |
Median :65.50 Median :65.0 Median :51.50 Median :56.50 Median :63.50
|
|
|
1754 |
Mean :64.63 Mean :66.6 Mean :53.13 Mean :56.37 Mean :64.63
|
|
|
1755 |
3rd Qu.:71.75 3rd Qu.:77.0 3rd Qu.:62.50 3rd Qu.:66.75 3rd Qu.:71.00
|
|
|
1756 |
Max. :85.00 Max. :90.0 Max. :83.00 Max. :75.00 Max. :88.00
|
|
|
1757 |
critical advance
|
|
|
1758 |
Min. :49.00 Min. :25.00
|
|
|
1759 |
1st Qu.:69.25 1st Qu.:35.00
|
|
|
1760 |
Median :77.50 Median :41.00
|
|
|
1761 |
Mean :74.77 Mean :42.93
|
|
|
1762 |
3rd Qu.:80.00 3rd Qu.:47.75
|
|
|
1763 |
Max. :92.00 Max. :72.00
|
|
|
1764 |
> summary(fm1 <- lm(rating ~ ., data = attitude))
|
|
|
1765 |
|
|
|
1766 |
Call:
|
|
|
1767 |
lm(formula = rating ~ ., data = attitude)
|
|
|
1768 |
|
|
|
1769 |
Residuals:
|
|
|
1770 |
Min 1Q Median 3Q Max
|
|
|
1771 |
-10.9418 -4.3555 0.3158 5.5425 11.5990
|
|
|
1772 |
|
|
|
1773 |
Coefficients:
|
|
|
1774 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1775 |
(Intercept) 10.78708 11.58926 0.931 0.361634
|
|
|
1776 |
complaints 0.61319 0.16098 3.809 0.000903 ***
|
|
|
1777 |
privileges -0.07305 0.13572 -0.538 0.595594
|
|
|
1778 |
learning 0.32033 0.16852 1.901 0.069925 .
|
|
|
1779 |
raises 0.08173 0.22148 0.369 0.715480
|
|
|
1780 |
critical 0.03838 0.14700 0.261 0.796334
|
|
|
1781 |
advance -0.21706 0.17821 -1.218 0.235577
|
|
|
1782 |
---
|
| 61435 |
ripley |
1783 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1784 |
|
|
|
1785 |
Residual standard error: 7.068 on 23 degrees of freedom
|
| 61435 |
ripley |
1786 |
Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
|
|
|
1787 |
F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05
|
| 56186 |
murdoch |
1788 |
|
|
|
1789 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
1790 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
1791 |
> plot(fm1)
|
|
|
1792 |
> summary(fm2 <- lm(rating ~ complaints, data = attitude))
|
|
|
1793 |
|
|
|
1794 |
Call:
|
|
|
1795 |
lm(formula = rating ~ complaints, data = attitude)
|
|
|
1796 |
|
|
|
1797 |
Residuals:
|
|
|
1798 |
Min 1Q Median 3Q Max
|
|
|
1799 |
-12.8799 -5.9905 0.1783 6.2978 9.6294
|
|
|
1800 |
|
|
|
1801 |
Coefficients:
|
|
|
1802 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1803 |
(Intercept) 14.37632 6.61999 2.172 0.0385 *
|
|
|
1804 |
complaints 0.75461 0.09753 7.737 1.99e-08 ***
|
|
|
1805 |
---
|
| 61435 |
ripley |
1806 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1807 |
|
|
|
1808 |
Residual standard error: 6.993 on 28 degrees of freedom
|
| 61435 |
ripley |
1809 |
Multiple R-squared: 0.6813, Adjusted R-squared: 0.6699
|
|
|
1810 |
F-statistic: 59.86 on 1 and 28 DF, p-value: 1.988e-08
|
| 56186 |
murdoch |
1811 |
|
|
|
1812 |
> plot(fm2)
|
|
|
1813 |
> par(opar)
|
|
|
1814 |
>
|
|
|
1815 |
>
|
|
|
1816 |
>
|
|
|
1817 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1818 |
> cleanEx()
|
|
|
1819 |
> nameEx("beavers")
|
|
|
1820 |
> ### * beavers
|
|
|
1821 |
>
|
|
|
1822 |
> flush(stderr()); flush(stdout())
|
|
|
1823 |
>
|
|
|
1824 |
> ### Name: beavers
|
|
|
1825 |
> ### Title: Body Temperature Series of Two Beavers
|
|
|
1826 |
> ### Aliases: beavers beaver1 beaver2
|
|
|
1827 |
> ### Keywords: datasets
|
|
|
1828 |
>
|
|
|
1829 |
> ### ** Examples
|
|
|
1830 |
>
|
|
|
1831 |
> require(graphics)
|
|
|
1832 |
> (yl <- range(beaver1$temp, beaver2$temp))
|
|
|
1833 |
[1] 36.33 38.35
|
|
|
1834 |
>
|
|
|
1835 |
> beaver.plot <- function(bdat, ...) {
|
|
|
1836 |
+ nam <- deparse(substitute(bdat))
|
|
|
1837 |
+ with(bdat, {
|
|
|
1838 |
+ # Hours since start of day:
|
|
|
1839 |
+ hours <- time %/% 100 + 24*(day - day[1]) + (time %% 100)/60
|
|
|
1840 |
+ plot (hours, temp, type = "l", ...,
|
|
|
1841 |
+ main = paste(nam, "body temperature"))
|
|
|
1842 |
+ abline(h = 37.5, col = "gray", lty = 2)
|
|
|
1843 |
+ is.act <- activ == 1
|
|
|
1844 |
+ points(hours[is.act], temp[is.act], col = 2, cex = .8)
|
|
|
1845 |
+ })
|
|
|
1846 |
+ }
|
| 61169 |
ripley |
1847 |
> op <- par(mfrow = c(2, 1), mar = c(3, 3, 4, 2), mgp = 0.9 * 2:0)
|
| 56186 |
murdoch |
1848 |
> beaver.plot(beaver1, ylim = yl)
|
|
|
1849 |
> beaver.plot(beaver2, ylim = yl)
|
|
|
1850 |
> par(op)
|
|
|
1851 |
>
|
|
|
1852 |
>
|
|
|
1853 |
>
|
|
|
1854 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1855 |
> cleanEx()
|
|
|
1856 |
> nameEx("cars")
|
|
|
1857 |
> ### * cars
|
|
|
1858 |
>
|
|
|
1859 |
> flush(stderr()); flush(stdout())
|
|
|
1860 |
>
|
|
|
1861 |
> ### Name: cars
|
|
|
1862 |
> ### Title: Speed and Stopping Distances of Cars
|
|
|
1863 |
> ### Aliases: cars
|
|
|
1864 |
> ### Keywords: datasets
|
|
|
1865 |
>
|
|
|
1866 |
> ### ** Examples
|
|
|
1867 |
>
|
|
|
1868 |
> require(stats); require(graphics)
|
|
|
1869 |
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
|
|
|
1870 |
+ las = 1)
|
|
|
1871 |
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
|
|
|
1872 |
> title(main = "cars data")
|
|
|
1873 |
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
|
|
|
1874 |
+ las = 1, log = "xy")
|
|
|
1875 |
> title(main = "cars data (logarithmic scales)")
|
|
|
1876 |
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
|
|
|
1877 |
> summary(fm1 <- lm(log(dist) ~ log(speed), data = cars))
|
|
|
1878 |
|
|
|
1879 |
Call:
|
|
|
1880 |
lm(formula = log(dist) ~ log(speed), data = cars)
|
|
|
1881 |
|
|
|
1882 |
Residuals:
|
|
|
1883 |
Min 1Q Median 3Q Max
|
|
|
1884 |
-1.00215 -0.24578 -0.02898 0.20717 0.88289
|
|
|
1885 |
|
|
|
1886 |
Coefficients:
|
|
|
1887 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
1888 |
(Intercept) -0.7297 0.3758 -1.941 0.0581 .
|
|
|
1889 |
log(speed) 1.6024 0.1395 11.484 2.26e-15 ***
|
|
|
1890 |
---
|
| 61435 |
ripley |
1891 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1892 |
|
|
|
1893 |
Residual standard error: 0.4053 on 48 degrees of freedom
|
| 61435 |
ripley |
1894 |
Multiple R-squared: 0.7331, Adjusted R-squared: 0.7276
|
|
|
1895 |
F-statistic: 131.9 on 1 and 48 DF, p-value: 2.259e-15
|
| 56186 |
murdoch |
1896 |
|
|
|
1897 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
1898 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
1899 |
> plot(fm1)
|
|
|
1900 |
> par(opar)
|
|
|
1901 |
>
|
|
|
1902 |
> ## An example of polynomial regression
|
|
|
1903 |
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
|
|
|
1904 |
+ las = 1, xlim = c(0, 25))
|
|
|
1905 |
> d <- seq(0, 25, length.out = 200)
|
|
|
1906 |
> for(degree in 1:4) {
|
|
|
1907 |
+ fm <- lm(dist ~ poly(speed, degree), data = cars)
|
| 61157 |
ripley |
1908 |
+ assign(paste("cars", degree, sep = "."), fm)
|
|
|
1909 |
+ lines(d, predict(fm, data.frame(speed = d)), col = degree)
|
| 56186 |
murdoch |
1910 |
+ }
|
|
|
1911 |
> anova(cars.1, cars.2, cars.3, cars.4)
|
|
|
1912 |
Analysis of Variance Table
|
|
|
1913 |
|
|
|
1914 |
Model 1: dist ~ poly(speed, degree)
|
|
|
1915 |
Model 2: dist ~ poly(speed, degree)
|
|
|
1916 |
Model 3: dist ~ poly(speed, degree)
|
|
|
1917 |
Model 4: dist ~ poly(speed, degree)
|
|
|
1918 |
Res.Df RSS Df Sum of Sq F Pr(>F)
|
|
|
1919 |
1 48 11354
|
|
|
1920 |
2 47 10825 1 528.81 2.3108 0.1355
|
|
|
1921 |
3 46 10634 1 190.35 0.8318 0.3666
|
|
|
1922 |
4 45 10298 1 336.55 1.4707 0.2316
|
|
|
1923 |
>
|
|
|
1924 |
>
|
|
|
1925 |
>
|
|
|
1926 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1927 |
> cleanEx()
|
|
|
1928 |
> nameEx("chickwts")
|
|
|
1929 |
> ### * chickwts
|
|
|
1930 |
>
|
|
|
1931 |
> flush(stderr()); flush(stdout())
|
|
|
1932 |
>
|
|
|
1933 |
> ### Name: chickwts
|
|
|
1934 |
> ### Title: Chicken Weights by Feed Type
|
|
|
1935 |
> ### Aliases: chickwts
|
|
|
1936 |
> ### Keywords: datasets
|
|
|
1937 |
>
|
|
|
1938 |
> ### ** Examples
|
|
|
1939 |
>
|
|
|
1940 |
> require(stats); require(graphics)
|
|
|
1941 |
> boxplot(weight ~ feed, data = chickwts, col = "lightgray",
|
|
|
1942 |
+ varwidth = TRUE, notch = TRUE, main = "chickwt data",
|
|
|
1943 |
+ ylab = "Weight at six weeks (gm)")
|
|
|
1944 |
Warning in bxp(list(stats = c(216, 271.5, 342, 373.5, 404, 108, 136, 151.5, :
|
|
|
1945 |
some notches went outside hinges ('box'): maybe set notch=FALSE
|
|
|
1946 |
> anova(fm1 <- lm(weight ~ feed, data = chickwts))
|
|
|
1947 |
Analysis of Variance Table
|
|
|
1948 |
|
|
|
1949 |
Response: weight
|
|
|
1950 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
1951 |
feed 5 231129 46226 15.365 5.936e-10 ***
|
|
|
1952 |
Residuals 65 195556 3009
|
|
|
1953 |
---
|
| 61435 |
ripley |
1954 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
1955 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
1956 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
1957 |
> plot(fm1)
|
|
|
1958 |
> par(opar)
|
|
|
1959 |
>
|
|
|
1960 |
>
|
|
|
1961 |
>
|
|
|
1962 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
1963 |
> cleanEx()
|
|
|
1964 |
> nameEx("co2")
|
|
|
1965 |
> ### * co2
|
|
|
1966 |
>
|
|
|
1967 |
> flush(stderr()); flush(stdout())
|
|
|
1968 |
>
|
|
|
1969 |
> ### Name: co2
|
|
|
1970 |
> ### Title: Mauna Loa Atmospheric CO2 Concentration
|
|
|
1971 |
> ### Aliases: co2
|
|
|
1972 |
> ### Keywords: datasets
|
|
|
1973 |
>
|
|
|
1974 |
> ### ** Examples
|
|
|
1975 |
>
|
|
|
1976 |
> require(graphics)
|
|
|
1977 |
> plot(co2, ylab = expression("Atmospheric concentration of CO"[2]),
|
|
|
1978 |
+ las = 1)
|
|
|
1979 |
> title(main = "co2 data set")
|
|
|
1980 |
>
|
|
|
1981 |
>
|
|
|
1982 |
>
|
|
|
1983 |
> cleanEx()
|
|
|
1984 |
> nameEx("crimtab")
|
|
|
1985 |
> ### * crimtab
|
|
|
1986 |
>
|
|
|
1987 |
> flush(stderr()); flush(stdout())
|
|
|
1988 |
>
|
|
|
1989 |
> ### Name: crimtab
|
|
|
1990 |
> ### Title: Student's 3000 Criminals Data
|
|
|
1991 |
> ### Aliases: crimtab
|
|
|
1992 |
> ### Keywords: datasets
|
|
|
1993 |
>
|
|
|
1994 |
> ### ** Examples
|
|
|
1995 |
>
|
|
|
1996 |
> require(stats)
|
|
|
1997 |
> dim(crimtab)
|
|
|
1998 |
[1] 42 22
|
|
|
1999 |
> utils::str(crimtab)
|
|
|
2000 |
'table' int [1:42, 1:22] 0 0 0 0 0 0 1 0 0 0 ...
|
|
|
2001 |
- attr(*, "dimnames")=List of 2
|
|
|
2002 |
..$ : chr [1:42] "9.4" "9.5" "9.6" "9.7" ...
|
|
|
2003 |
..$ : chr [1:22] "142.24" "144.78" "147.32" "149.86" ...
|
|
|
2004 |
> ## for nicer printing:
|
|
|
2005 |
> local({cT <- crimtab
|
| 61169 |
ripley |
2006 |
+ colnames(cT) <- substring(colnames(cT), 2, 3)
|
| 56186 |
murdoch |
2007 |
+ print(cT, zero.print = " ")
|
|
|
2008 |
+ })
|
|
|
2009 |
42 44 47 49 52 54 57 60 62 65 67 70 72 75 77 80 82 85 87 90 93 95
|
|
|
2010 |
9.4
|
|
|
2011 |
9.5 1
|
|
|
2012 |
9.6
|
|
|
2013 |
9.7
|
|
|
2014 |
9.8 1
|
|
|
2015 |
9.9 1 1 1
|
|
|
2016 |
10 1 1 2 2 1
|
|
|
2017 |
10.1 1 3 1 1 1
|
|
|
2018 |
10.2 2 2 2 1 2 1
|
|
|
2019 |
10.3 1 1 3 2 2 3 5
|
|
|
2020 |
10.4 1 1 2 3 3 4 3 3
|
|
|
2021 |
10.5 1 3 7 6 4 3 1 3 1 1
|
|
|
2022 |
10.6 1 4 5 9 14 6 3 1 1
|
|
|
2023 |
10.7 1 2 4 9 14 16 15 7 3 1 2
|
|
|
2024 |
10.8 2 5 6 14 27 10 7 1 2 1
|
|
|
2025 |
10.9 2 6 14 24 27 14 10 4 1
|
|
|
2026 |
11 2 6 12 15 31 37 27 17 10 6
|
|
|
2027 |
11.1 3 3 12 22 26 24 26 24 7 4 1
|
|
|
2028 |
11.2 3 2 7 21 30 38 29 27 20 4 1 1
|
|
|
2029 |
11.3 1 5 10 24 26 39 26 24 7 2
|
|
|
2030 |
11.4 3 4 9 29 56 58 26 22 10 11
|
|
|
2031 |
11.5 5 11 17 33 57 38 34 25 11 2
|
|
|
2032 |
11.6 2 1 4 13 37 39 48 38 27 12 2 2 1
|
|
|
2033 |
11.7 2 9 17 30 37 48 45 24 9 9 2
|
|
|
2034 |
11.8 1 2 11 15 35 41 34 29 10 5 1
|
|
|
2035 |
11.9 1 1 2 12 10 27 32 35 19 10 9 3 1
|
|
|
2036 |
12 1 4 8 19 42 39 22 16 8 2 2
|
|
|
2037 |
12.1 2 4 13 22 28 15 27 10 4 1
|
|
|
2038 |
12.2 1 2 5 6 23 17 16 11 8 1 1
|
|
|
2039 |
12.3 4 8 10 13 20 23 6 5
|
|
|
2040 |
12.4 1 1 1 2 7 12 4 7 7 1 1
|
|
|
2041 |
12.5 1 1 3 12 11 8 6 8 2
|
|
|
2042 |
12.6 1 3 5 7 8 6 3 1 1
|
|
|
2043 |
12.7 1 1 7 5 5 8 2 2
|
|
|
2044 |
12.8 1 2 3 1 8 5 3 1 1
|
|
|
2045 |
12.9 1 2 2 1 1
|
|
|
2046 |
13 3 1 1 2 1
|
|
|
2047 |
13.1 1 1
|
|
|
2048 |
13.2 1 1 1 3
|
|
|
2049 |
13.3 1 1
|
|
|
2050 |
13.4
|
|
|
2051 |
13.5 1
|
|
|
2052 |
>
|
|
|
2053 |
> ## Repeat Student's experiment:
|
|
|
2054 |
>
|
|
|
2055 |
> # 1) Reconstitute 3000 raw data for heights in inches and rounded to
|
|
|
2056 |
> # nearest integer as in Student's paper:
|
|
|
2057 |
>
|
|
|
2058 |
> (heIn <- round(as.numeric(colnames(crimtab)) / 2.54))
|
|
|
2059 |
[1] 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
|
|
|
2060 |
> d.hei <- data.frame(height = rep(heIn, colSums(crimtab)))
|
|
|
2061 |
>
|
|
|
2062 |
> # 2) shuffle the data:
|
|
|
2063 |
>
|
|
|
2064 |
> set.seed(1)
|
|
|
2065 |
> d.hei <- d.hei[sample(1:3000), , drop = FALSE]
|
|
|
2066 |
>
|
|
|
2067 |
> # 3) Make 750 samples each of size 4:
|
|
|
2068 |
>
|
|
|
2069 |
> d.hei$sample <- as.factor(rep(1:750, each = 4))
|
|
|
2070 |
>
|
|
|
2071 |
> # 4) Compute the means and standard deviations (n) for the 750 samples:
|
|
|
2072 |
>
|
|
|
2073 |
> h.mean <- with(d.hei, tapply(height, sample, FUN = mean))
|
|
|
2074 |
> h.sd <- with(d.hei, tapply(height, sample, FUN = sd)) * sqrt(3/4)
|
|
|
2075 |
>
|
|
|
2076 |
> # 5) Compute the difference between the mean of each sample and
|
|
|
2077 |
> # the mean of the population and then divide by the
|
|
|
2078 |
> # standard deviation of the sample:
|
|
|
2079 |
>
|
|
|
2080 |
> zobs <- (h.mean - mean(d.hei[,"height"]))/h.sd
|
|
|
2081 |
>
|
|
|
2082 |
> # 6) Replace infinite values by +/- 6 as in Student's paper:
|
|
|
2083 |
>
|
|
|
2084 |
> zobs[infZ <- is.infinite(zobs)] # 3 of them
|
|
|
2085 |
73 312 674
|
|
|
2086 |
-Inf Inf -Inf
|
|
|
2087 |
> zobs[infZ] <- 6 * sign(zobs[infZ])
|
|
|
2088 |
>
|
|
|
2089 |
> # 7) Plot the distribution:
|
|
|
2090 |
>
|
|
|
2091 |
> require(grDevices); require(graphics)
|
|
|
2092 |
> hist(x = zobs, probability = TRUE, xlab = "Student's z",
|
|
|
2093 |
+ col = grey(0.8), border = grey(0.5),
|
|
|
2094 |
+ main = "Distribution of Student's z score for 'crimtab' data")
|
|
|
2095 |
>
|
|
|
2096 |
>
|
|
|
2097 |
>
|
|
|
2098 |
> cleanEx()
|
|
|
2099 |
> nameEx("discoveries")
|
|
|
2100 |
> ### * discoveries
|
|
|
2101 |
>
|
|
|
2102 |
> flush(stderr()); flush(stdout())
|
|
|
2103 |
>
|
|
|
2104 |
> ### Name: discoveries
|
|
|
2105 |
> ### Title: Yearly Numbers of Important Discoveries
|
|
|
2106 |
> ### Aliases: discoveries
|
|
|
2107 |
> ### Keywords: datasets
|
|
|
2108 |
>
|
|
|
2109 |
> ### ** Examples
|
|
|
2110 |
>
|
|
|
2111 |
> require(graphics)
|
|
|
2112 |
> plot(discoveries, ylab = "Number of important discoveries",
|
|
|
2113 |
+ las = 1)
|
|
|
2114 |
> title(main = "discoveries data set")
|
|
|
2115 |
>
|
|
|
2116 |
>
|
|
|
2117 |
>
|
|
|
2118 |
> cleanEx()
|
|
|
2119 |
> nameEx("esoph")
|
|
|
2120 |
> ### * esoph
|
|
|
2121 |
>
|
|
|
2122 |
> flush(stderr()); flush(stdout())
|
|
|
2123 |
>
|
|
|
2124 |
> ### Name: esoph
|
|
|
2125 |
> ### Title: Smoking, Alcohol and (O)esophageal Cancer
|
|
|
2126 |
> ### Aliases: esoph
|
|
|
2127 |
> ### Keywords: datasets
|
|
|
2128 |
>
|
|
|
2129 |
> ### ** Examples
|
|
|
2130 |
>
|
|
|
2131 |
> require(stats)
|
|
|
2132 |
> require(graphics) # for mosaicplot
|
|
|
2133 |
> summary(esoph)
|
|
|
2134 |
agegp alcgp tobgp ncases ncontrols
|
|
|
2135 |
25-34:15 0-39g/day:23 0-9g/day:24 Min. : 0.000 Min. : 1.00
|
|
|
2136 |
35-44:15 40-79 :23 10-19 :24 1st Qu.: 0.000 1st Qu.: 3.00
|
|
|
2137 |
45-54:16 80-119 :21 20-29 :20 Median : 1.000 Median : 6.00
|
|
|
2138 |
55-64:16 120+ :21 30+ :20 Mean : 2.273 Mean :11.08
|
|
|
2139 |
65-74:15 3rd Qu.: 4.000 3rd Qu.:14.00
|
|
|
2140 |
75+ :11 Max. :17.000 Max. :60.00
|
|
|
2141 |
> ## effects of alcohol, tobacco and interaction, age-adjusted
|
|
|
2142 |
> model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp,
|
|
|
2143 |
+ data = esoph, family = binomial())
|
|
|
2144 |
> anova(model1)
|
|
|
2145 |
Analysis of Deviance Table
|
|
|
2146 |
|
|
|
2147 |
Model: binomial, link: logit
|
|
|
2148 |
|
|
|
2149 |
Response: cbind(ncases, ncontrols)
|
|
|
2150 |
|
|
|
2151 |
Terms added sequentially (first to last)
|
|
|
2152 |
|
|
|
2153 |
|
|
|
2154 |
Df Deviance Resid. Df Resid. Dev
|
|
|
2155 |
NULL 87 227.241
|
|
|
2156 |
agegp 5 88.128 82 139.112
|
|
|
2157 |
tobgp 3 19.085 79 120.028
|
|
|
2158 |
alcgp 3 66.054 76 53.973
|
|
|
2159 |
tobgp:alcgp 9 6.489 67 47.484
|
|
|
2160 |
> ## Try a linear effect of alcohol and tobacco
|
|
|
2161 |
> model2 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp)
|
|
|
2162 |
+ + unclass(alcgp),
|
|
|
2163 |
+ data = esoph, family = binomial())
|
|
|
2164 |
> summary(model2)
|
|
|
2165 |
|
|
|
2166 |
Call:
|
|
|
2167 |
glm(formula = cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) +
|
|
|
2168 |
unclass(alcgp), family = binomial(), data = esoph)
|
|
|
2169 |
|
|
|
2170 |
Deviance Residuals:
|
|
|
2171 |
Min 1Q Median 3Q Max
|
|
|
2172 |
-1.7628 -0.6426 -0.2709 0.3043 2.0421
|
|
|
2173 |
|
|
|
2174 |
Coefficients:
|
|
|
2175 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
2176 |
(Intercept) -4.01097 0.31224 -12.846 < 2e-16 ***
|
|
|
2177 |
agegp.L 2.96113 0.65092 4.549 5.39e-06 ***
|
|
|
2178 |
agegp.Q -1.33735 0.58918 -2.270 0.02322 *
|
|
|
2179 |
agegp.C 0.15292 0.44792 0.341 0.73281
|
|
|
2180 |
agegp^4 0.06668 0.30776 0.217 0.82848
|
|
|
2181 |
agegp^5 -0.20288 0.19523 -1.039 0.29872
|
|
|
2182 |
unclass(tobgp) 0.26162 0.08198 3.191 0.00142 **
|
|
|
2183 |
unclass(alcgp) 0.65308 0.08452 7.727 1.10e-14 ***
|
|
|
2184 |
---
|
| 61435 |
ripley |
2185 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2186 |
|
|
|
2187 |
(Dispersion parameter for binomial family taken to be 1)
|
|
|
2188 |
|
|
|
2189 |
Null deviance: 227.241 on 87 degrees of freedom
|
|
|
2190 |
Residual deviance: 59.277 on 80 degrees of freedom
|
|
|
2191 |
AIC: 222.76
|
|
|
2192 |
|
|
|
2193 |
Number of Fisher Scoring iterations: 6
|
|
|
2194 |
|
|
|
2195 |
> ## Re-arrange data for a mosaic plot
|
|
|
2196 |
> ttt <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
|
| 57020 |
murdoch |
2197 |
> o <- with(esoph, order(tobgp, alcgp, agegp))
|
|
|
2198 |
> ttt[ttt == 1] <- esoph$ncases[o]
|
| 56186 |
murdoch |
2199 |
> tt1 <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
|
| 57020 |
murdoch |
2200 |
> tt1[tt1 == 1] <- esoph$ncontrols[o]
|
| 56186 |
murdoch |
2201 |
> tt <- array(c(ttt, tt1), c(dim(ttt),2),
|
|
|
2202 |
+ c(dimnames(ttt), list(c("Cancer", "control"))))
|
|
|
2203 |
> mosaicplot(tt, main = "esoph data set", color = TRUE)
|
|
|
2204 |
>
|
|
|
2205 |
>
|
|
|
2206 |
>
|
|
|
2207 |
> cleanEx()
|
|
|
2208 |
> nameEx("euro")
|
|
|
2209 |
> ### * euro
|
|
|
2210 |
>
|
|
|
2211 |
> flush(stderr()); flush(stdout())
|
|
|
2212 |
>
|
|
|
2213 |
> ### Name: euro
|
|
|
2214 |
> ### Title: Conversion Rates of Euro Currencies
|
|
|
2215 |
> ### Aliases: euro euro.cross
|
|
|
2216 |
> ### Keywords: datasets
|
|
|
2217 |
>
|
|
|
2218 |
> ### ** Examples
|
|
|
2219 |
>
|
|
|
2220 |
> cbind(euro)
|
|
|
2221 |
euro
|
|
|
2222 |
ATS 13.760300
|
|
|
2223 |
BEF 40.339900
|
|
|
2224 |
DEM 1.955830
|
|
|
2225 |
ESP 166.386000
|
|
|
2226 |
FIM 5.945730
|
|
|
2227 |
FRF 6.559570
|
|
|
2228 |
IEP 0.787564
|
|
|
2229 |
ITL 1936.270000
|
|
|
2230 |
LUF 40.339900
|
|
|
2231 |
NLG 2.203710
|
|
|
2232 |
PTE 200.482000
|
|
|
2233 |
>
|
|
|
2234 |
> ## These relations hold:
|
| 61169 |
ripley |
2235 |
> euro == signif(euro, 6) # [6 digit precision in Euro's definition]
|
| 56186 |
murdoch |
2236 |
ATS BEF DEM ESP FIM FRF IEP ITL LUF NLG PTE
|
|
|
2237 |
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
|
|
|
2238 |
> all(euro.cross == outer(1/euro, euro))
|
|
|
2239 |
[1] TRUE
|
|
|
2240 |
>
|
|
|
2241 |
> ## Convert 20 Euro to Belgian Franc
|
|
|
2242 |
> 20 * euro["BEF"]
|
|
|
2243 |
BEF
|
|
|
2244 |
806.798
|
|
|
2245 |
> ## Convert 20 Austrian Schilling to Euro
|
|
|
2246 |
> 20 / euro["ATS"]
|
|
|
2247 |
ATS
|
|
|
2248 |
1.453457
|
|
|
2249 |
> ## Convert 20 Spanish Pesetas to Italian Lira
|
|
|
2250 |
> 20 * euro.cross["ESP", "ITL"]
|
|
|
2251 |
[1] 232.7443
|
|
|
2252 |
>
|
|
|
2253 |
> require(graphics)
|
|
|
2254 |
> dotchart(euro,
|
|
|
2255 |
+ main = "euro data: 1 Euro in currency unit")
|
|
|
2256 |
> dotchart(1/euro,
|
|
|
2257 |
+ main = "euro data: 1 currency unit in Euros")
|
|
|
2258 |
> dotchart(log(euro, 10),
|
|
|
2259 |
+ main = "euro data: log10(1 Euro in currency unit)")
|
|
|
2260 |
>
|
|
|
2261 |
>
|
|
|
2262 |
>
|
|
|
2263 |
> cleanEx()
|
|
|
2264 |
> nameEx("faithful")
|
|
|
2265 |
> ### * faithful
|
|
|
2266 |
>
|
|
|
2267 |
> flush(stderr()); flush(stdout())
|
|
|
2268 |
>
|
|
|
2269 |
> ### Name: faithful
|
|
|
2270 |
> ### Title: Old Faithful Geyser Data
|
|
|
2271 |
> ### Aliases: faithful
|
|
|
2272 |
> ### Keywords: datasets
|
|
|
2273 |
>
|
|
|
2274 |
> ### ** Examples
|
|
|
2275 |
>
|
|
|
2276 |
> require(stats); require(graphics)
|
|
|
2277 |
> f.tit <- "faithful data: Eruptions of Old Faithful"
|
|
|
2278 |
>
|
|
|
2279 |
> ne60 <- round(e60 <- 60 * faithful$eruptions)
|
|
|
2280 |
> all.equal(e60, ne60) # relative diff. ~ 1/10000
|
|
|
2281 |
[1] "Mean relative difference: 9.515332e-05"
|
|
|
2282 |
> table(zapsmall(abs(e60 - ne60))) # 0, 0.02 or 0.04
|
|
|
2283 |
|
|
|
2284 |
|
|
|
2285 |
106 163 3
|
|
|
2286 |
> faithful$better.eruptions <- ne60 / 60
|
|
|
2287 |
> te <- table(ne60)
|
|
|
2288 |
> te[te >= 4] # (too) many multiples of 5 !
|
|
|
2289 |
ne60
|
|
|
2290 |
105 108 110 112 113 120 216 230 240 245 249 250 255 260 261 262 265 270 272 275
|
|
|
2291 |
6 4 7 8 4 4 4 5 6 5 4 4 4 5 4 4 4 8 5 4
|
|
|
2292 |
276 282 288
|
|
|
2293 |
4 6 6
|
| 61157 |
ripley |
2294 |
> plot(names(te), te, type = "h", main = f.tit, xlab = "Eruption time (sec)")
|
| 56186 |
murdoch |
2295 |
>
|
|
|
2296 |
> plot(faithful[, -3], main = f.tit,
|
|
|
2297 |
+ xlab = "Eruption time (min)",
|
|
|
2298 |
+ ylab = "Waiting time to next eruption (min)")
|
|
|
2299 |
> lines(lowess(faithful$eruptions, faithful$waiting, f = 2/3, iter = 3),
|
|
|
2300 |
+ col = "red")
|
|
|
2301 |
>
|
|
|
2302 |
>
|
|
|
2303 |
>
|
|
|
2304 |
> cleanEx()
|
|
|
2305 |
> nameEx("freeny")
|
|
|
2306 |
> ### * freeny
|
|
|
2307 |
>
|
|
|
2308 |
> flush(stderr()); flush(stdout())
|
|
|
2309 |
>
|
|
|
2310 |
> ### Name: freeny
|
|
|
2311 |
> ### Title: Freeny's Revenue Data
|
|
|
2312 |
> ### Aliases: freeny freeny.x freeny.y
|
|
|
2313 |
> ### Keywords: datasets
|
|
|
2314 |
>
|
|
|
2315 |
> ### ** Examples
|
|
|
2316 |
>
|
|
|
2317 |
> require(stats); require(graphics)
|
|
|
2318 |
> summary(freeny)
|
|
|
2319 |
y lag.quarterly.revenue price.index income.level
|
|
|
2320 |
Min. :8.791 Min. :8.791 Min. :4.278 Min. :5.821
|
|
|
2321 |
1st Qu.:9.045 1st Qu.:9.020 1st Qu.:4.392 1st Qu.:5.948
|
|
|
2322 |
Median :9.314 Median :9.284 Median :4.510 Median :6.061
|
|
|
2323 |
Mean :9.306 Mean :9.281 Mean :4.496 Mean :6.039
|
|
|
2324 |
3rd Qu.:9.591 3rd Qu.:9.561 3rd Qu.:4.605 3rd Qu.:6.139
|
|
|
2325 |
Max. :9.794 Max. :9.775 Max. :4.710 Max. :6.200
|
|
|
2326 |
market.potential
|
|
|
2327 |
Min. :12.97
|
|
|
2328 |
1st Qu.:13.01
|
|
|
2329 |
Median :13.07
|
|
|
2330 |
Mean :13.07
|
|
|
2331 |
3rd Qu.:13.12
|
|
|
2332 |
Max. :13.17
|
|
|
2333 |
> pairs(freeny, main = "freeny data")
|
|
|
2334 |
> # gives warning: freeny$y has class "ts"
|
|
|
2335 |
>
|
|
|
2336 |
> summary(fm1 <- lm(y ~ ., data = freeny))
|
|
|
2337 |
|
|
|
2338 |
Call:
|
|
|
2339 |
lm(formula = y ~ ., data = freeny)
|
|
|
2340 |
|
|
|
2341 |
Residuals:
|
|
|
2342 |
Min 1Q Median 3Q Max
|
|
|
2343 |
-0.0259426 -0.0101033 0.0003824 0.0103236 0.0267124
|
|
|
2344 |
|
|
|
2345 |
Coefficients:
|
|
|
2346 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
2347 |
(Intercept) -10.4726 6.0217 -1.739 0.0911 .
|
|
|
2348 |
lag.quarterly.revenue 0.1239 0.1424 0.870 0.3904
|
|
|
2349 |
price.index -0.7542 0.1607 -4.693 4.28e-05 ***
|
|
|
2350 |
income.level 0.7675 0.1339 5.730 1.93e-06 ***
|
|
|
2351 |
market.potential 1.3306 0.5093 2.613 0.0133 *
|
|
|
2352 |
---
|
| 61435 |
ripley |
2353 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2354 |
|
|
|
2355 |
Residual standard error: 0.01473 on 34 degrees of freedom
|
| 61435 |
ripley |
2356 |
Multiple R-squared: 0.9981, Adjusted R-squared: 0.9978
|
|
|
2357 |
F-statistic: 4354 on 4 and 34 DF, p-value: < 2.2e-16
|
| 56186 |
murdoch |
2358 |
|
|
|
2359 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
2360 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
2361 |
> plot(fm1)
|
|
|
2362 |
> par(opar)
|
|
|
2363 |
>
|
|
|
2364 |
>
|
|
|
2365 |
>
|
|
|
2366 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
2367 |
> cleanEx()
|
|
|
2368 |
> nameEx("infert")
|
|
|
2369 |
> ### * infert
|
|
|
2370 |
>
|
|
|
2371 |
> flush(stderr()); flush(stdout())
|
|
|
2372 |
>
|
|
|
2373 |
> ### Name: infert
|
|
|
2374 |
> ### Title: Infertility after Spontaneous and Induced Abortion
|
|
|
2375 |
> ### Aliases: infert
|
|
|
2376 |
> ### Keywords: datasets
|
|
|
2377 |
>
|
|
|
2378 |
> ### ** Examples
|
|
|
2379 |
>
|
|
|
2380 |
> require(stats)
|
| 61157 |
ripley |
2381 |
> model1 <- glm(case ~ spontaneous+induced, data = infert, family = binomial())
|
| 56186 |
murdoch |
2382 |
> summary(model1)
|
|
|
2383 |
|
|
|
2384 |
Call:
|
|
|
2385 |
glm(formula = case ~ spontaneous + induced, family = binomial(),
|
|
|
2386 |
data = infert)
|
|
|
2387 |
|
|
|
2388 |
Deviance Residuals:
|
|
|
2389 |
Min 1Q Median 3Q Max
|
|
|
2390 |
-1.6678 -0.8360 -0.5772 0.9030 1.9362
|
|
|
2391 |
|
|
|
2392 |
Coefficients:
|
|
|
2393 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
2394 |
(Intercept) -1.7079 0.2677 -6.380 1.78e-10 ***
|
|
|
2395 |
spontaneous 1.1972 0.2116 5.657 1.54e-08 ***
|
|
|
2396 |
induced 0.4181 0.2056 2.033 0.042 *
|
|
|
2397 |
---
|
| 61435 |
ripley |
2398 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2399 |
|
|
|
2400 |
(Dispersion parameter for binomial family taken to be 1)
|
|
|
2401 |
|
|
|
2402 |
Null deviance: 316.17 on 247 degrees of freedom
|
|
|
2403 |
Residual deviance: 279.61 on 245 degrees of freedom
|
|
|
2404 |
AIC: 285.61
|
|
|
2405 |
|
|
|
2406 |
Number of Fisher Scoring iterations: 4
|
|
|
2407 |
|
|
|
2408 |
> ## adjusted for other potential confounders:
|
|
|
2409 |
> summary(model2 <- glm(case ~ age+parity+education+spontaneous+induced,
|
| 61157 |
ripley |
2410 |
+ data = infert, family = binomial()))
|
| 56186 |
murdoch |
2411 |
|
|
|
2412 |
Call:
|
|
|
2413 |
glm(formula = case ~ age + parity + education + spontaneous +
|
|
|
2414 |
induced, family = binomial(), data = infert)
|
|
|
2415 |
|
|
|
2416 |
Deviance Residuals:
|
|
|
2417 |
Min 1Q Median 3Q Max
|
|
|
2418 |
-1.7603 -0.8162 -0.4956 0.8349 2.6536
|
|
|
2419 |
|
|
|
2420 |
Coefficients:
|
|
|
2421 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
2422 |
(Intercept) -1.14924 1.41220 -0.814 0.4158
|
|
|
2423 |
age 0.03958 0.03120 1.269 0.2046
|
|
|
2424 |
parity -0.82828 0.19649 -4.215 2.49e-05 ***
|
|
|
2425 |
education6-11yrs -1.04424 0.79255 -1.318 0.1876
|
|
|
2426 |
education12+ yrs -1.40321 0.83416 -1.682 0.0925 .
|
|
|
2427 |
spontaneous 2.04591 0.31016 6.596 4.21e-11 ***
|
|
|
2428 |
induced 1.28876 0.30146 4.275 1.91e-05 ***
|
|
|
2429 |
---
|
| 61435 |
ripley |
2430 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2431 |
|
|
|
2432 |
(Dispersion parameter for binomial family taken to be 1)
|
|
|
2433 |
|
|
|
2434 |
Null deviance: 316.17 on 247 degrees of freedom
|
|
|
2435 |
Residual deviance: 257.80 on 241 degrees of freedom
|
|
|
2436 |
AIC: 271.8
|
|
|
2437 |
|
|
|
2438 |
Number of Fisher Scoring iterations: 4
|
|
|
2439 |
|
|
|
2440 |
> ## Really should be analysed by conditional logistic regression
|
|
|
2441 |
> ## which is in the survival package
|
|
|
2442 |
> if(require(survival)){
|
| 61157 |
ripley |
2443 |
+ model3 <- clogit(case ~ spontaneous+induced+strata(stratum), data = infert)
|
| 56186 |
murdoch |
2444 |
+ print(summary(model3))
|
| 61169 |
ripley |
2445 |
+ detach() # survival (conflicts)
|
| 56186 |
murdoch |
2446 |
+ }
|
|
|
2447 |
Loading required package: survival
|
|
|
2448 |
Loading required package: splines
|
|
|
2449 |
Call:
|
|
|
2450 |
coxph(formula = Surv(rep(1, 248L), case) ~ spontaneous + induced +
|
|
|
2451 |
strata(stratum), data = infert, method = "exact")
|
|
|
2452 |
|
|
|
2453 |
n= 248, number of events= 83
|
|
|
2454 |
|
|
|
2455 |
coef exp(coef) se(coef) z Pr(>|z|)
|
|
|
2456 |
spontaneous 1.9859 7.2854 0.3524 5.635 1.75e-08 ***
|
|
|
2457 |
induced 1.4090 4.0919 0.3607 3.906 9.38e-05 ***
|
|
|
2458 |
---
|
| 61435 |
ripley |
2459 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2460 |
|
|
|
2461 |
exp(coef) exp(-coef) lower .95 upper .95
|
|
|
2462 |
spontaneous 7.285 0.1373 3.651 14.536
|
|
|
2463 |
induced 4.092 0.2444 2.018 8.298
|
|
|
2464 |
|
|
|
2465 |
Rsquare= 0.193 (max possible= 0.519 )
|
|
|
2466 |
Likelihood ratio test= 53.15 on 2 df, p=2.869e-12
|
|
|
2467 |
Wald test = 31.84 on 2 df, p=1.221e-07
|
|
|
2468 |
Score (logrank) test = 48.44 on 2 df, p=3.032e-11
|
|
|
2469 |
|
|
|
2470 |
>
|
|
|
2471 |
>
|
|
|
2472 |
>
|
|
|
2473 |
> cleanEx()
|
|
|
2474 |
|
|
|
2475 |
detaching ‘package:splines’
|
|
|
2476 |
|
|
|
2477 |
> nameEx("iris")
|
|
|
2478 |
> ### * iris
|
|
|
2479 |
>
|
|
|
2480 |
> flush(stderr()); flush(stdout())
|
|
|
2481 |
>
|
|
|
2482 |
> ### Name: iris
|
|
|
2483 |
> ### Title: Edgar Anderson's Iris Data
|
|
|
2484 |
> ### Aliases: iris iris3
|
|
|
2485 |
> ### Keywords: datasets
|
|
|
2486 |
>
|
|
|
2487 |
> ### ** Examples
|
|
|
2488 |
>
|
|
|
2489 |
> dni3 <- dimnames(iris3)
|
| 61157 |
ripley |
2490 |
> ii <- data.frame(matrix(aperm(iris3, c(1,3,2)), ncol = 4,
|
| 56186 |
murdoch |
2491 |
+ dimnames = list(NULL, sub(" L.",".Length",
|
|
|
2492 |
+ sub(" W.",".Width", dni3[[2]])))),
|
| 61157 |
ripley |
2493 |
+ Species = gl(3, 50, labels = sub("S", "s", sub("V", "v", dni3[[3]]))))
|
| 56186 |
murdoch |
2494 |
> all.equal(ii, iris) # TRUE
|
|
|
2495 |
[1] TRUE
|
|
|
2496 |
>
|
|
|
2497 |
>
|
|
|
2498 |
>
|
|
|
2499 |
> cleanEx()
|
|
|
2500 |
> nameEx("islands")
|
|
|
2501 |
> ### * islands
|
|
|
2502 |
>
|
|
|
2503 |
> flush(stderr()); flush(stdout())
|
|
|
2504 |
>
|
|
|
2505 |
> ### Name: islands
|
|
|
2506 |
> ### Title: Areas of the World's Major Landmasses
|
|
|
2507 |
> ### Aliases: islands
|
|
|
2508 |
> ### Keywords: datasets
|
|
|
2509 |
>
|
|
|
2510 |
> ### ** Examples
|
|
|
2511 |
>
|
|
|
2512 |
> require(graphics)
|
|
|
2513 |
> dotchart(log(islands, 10),
|
|
|
2514 |
+ main = "islands data: log10(area) (log10(sq. miles))")
|
|
|
2515 |
> dotchart(log(islands[order(islands)], 10),
|
|
|
2516 |
+ main = "islands data: log10(area) (log10(sq. miles))")
|
|
|
2517 |
>
|
|
|
2518 |
>
|
|
|
2519 |
>
|
|
|
2520 |
> cleanEx()
|
|
|
2521 |
> nameEx("longley")
|
|
|
2522 |
> ### * longley
|
|
|
2523 |
>
|
|
|
2524 |
> flush(stderr()); flush(stdout())
|
|
|
2525 |
>
|
|
|
2526 |
> ### Name: longley
|
|
|
2527 |
> ### Title: Longley's Economic Regression Data
|
|
|
2528 |
> ### Aliases: longley
|
|
|
2529 |
> ### Keywords: datasets
|
|
|
2530 |
>
|
|
|
2531 |
> ### ** Examples
|
|
|
2532 |
>
|
|
|
2533 |
> require(stats); require(graphics)
|
|
|
2534 |
> ## give the data set in the form it is used in S-PLUS:
|
|
|
2535 |
> longley.x <- data.matrix(longley[, 1:6])
|
|
|
2536 |
> longley.y <- longley[, "Employed"]
|
|
|
2537 |
> pairs(longley, main = "longley data")
|
|
|
2538 |
> summary(fm1 <- lm(Employed ~ ., data = longley))
|
|
|
2539 |
|
|
|
2540 |
Call:
|
|
|
2541 |
lm(formula = Employed ~ ., data = longley)
|
|
|
2542 |
|
|
|
2543 |
Residuals:
|
|
|
2544 |
Min 1Q Median 3Q Max
|
|
|
2545 |
-0.41011 -0.15767 -0.02816 0.10155 0.45539
|
|
|
2546 |
|
|
|
2547 |
Coefficients:
|
|
|
2548 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
2549 |
(Intercept) -3.482e+03 8.904e+02 -3.911 0.003560 **
|
|
|
2550 |
GNP.deflator 1.506e-02 8.492e-02 0.177 0.863141
|
|
|
2551 |
GNP -3.582e-02 3.349e-02 -1.070 0.312681
|
|
|
2552 |
Unemployed -2.020e-02 4.884e-03 -4.136 0.002535 **
|
|
|
2553 |
Armed.Forces -1.033e-02 2.143e-03 -4.822 0.000944 ***
|
|
|
2554 |
Population -5.110e-02 2.261e-01 -0.226 0.826212
|
|
|
2555 |
Year 1.829e+00 4.555e-01 4.016 0.003037 **
|
|
|
2556 |
---
|
| 61435 |
ripley |
2557 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2558 |
|
|
|
2559 |
Residual standard error: 0.3049 on 9 degrees of freedom
|
| 61435 |
ripley |
2560 |
Multiple R-squared: 0.9955, Adjusted R-squared: 0.9925
|
|
|
2561 |
F-statistic: 330.3 on 6 and 9 DF, p-value: 4.984e-10
|
| 56186 |
murdoch |
2562 |
|
|
|
2563 |
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
|
|
|
2564 |
+ mar = c(4.1, 4.1, 2.1, 1.1))
|
|
|
2565 |
> plot(fm1)
|
|
|
2566 |
> par(opar)
|
|
|
2567 |
>
|
|
|
2568 |
>
|
|
|
2569 |
>
|
|
|
2570 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
2571 |
> cleanEx()
|
|
|
2572 |
> nameEx("morley")
|
|
|
2573 |
> ### * morley
|
|
|
2574 |
>
|
|
|
2575 |
> flush(stderr()); flush(stdout())
|
|
|
2576 |
>
|
|
|
2577 |
> ### Name: morley
|
| 59017 |
ripley |
2578 |
> ### Title: Michelson Speed of Light Data
|
| 56186 |
murdoch |
2579 |
> ### Aliases: morley
|
|
|
2580 |
> ### Keywords: datasets
|
|
|
2581 |
>
|
|
|
2582 |
> ### ** Examples
|
|
|
2583 |
>
|
|
|
2584 |
> require(stats); require(graphics)
|
| 58992 |
maechler |
2585 |
> michelson <- transform(morley,
|
|
|
2586 |
+ Expt = factor(Expt), Run = factor(Run))
|
| 61169 |
ripley |
2587 |
> xtabs(~ Expt + Run, data = michelson) # 5 x 20 balanced (two-way)
|
| 56186 |
murdoch |
2588 |
Run
|
|
|
2589 |
Expt 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
|
|
|
2590 |
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2591 |
2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2592 |
3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2593 |
4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
|
|
2594 |
5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
|
| 58992 |
maechler |
2595 |
> plot(Speed ~ Expt, data = michelson,
|
| 56186 |
murdoch |
2596 |
+ main = "Speed of Light Data", xlab = "Experiment No.")
|
| 58992 |
maechler |
2597 |
> fm <- aov(Speed ~ Run + Expt, data = michelson)
|
| 56186 |
murdoch |
2598 |
> summary(fm)
|
| 57044 |
ripley |
2599 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
2600 |
Run 19 113344 5965 1.105 0.36321
|
|
|
2601 |
Expt 4 94514 23629 4.378 0.00307 **
|
|
|
2602 |
Residuals 76 410166 5397
|
| 56186 |
murdoch |
2603 |
---
|
| 61435 |
ripley |
2604 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2605 |
> fm0 <- update(fm, . ~ . - Run)
|
|
|
2606 |
> anova(fm0, fm)
|
|
|
2607 |
Analysis of Variance Table
|
|
|
2608 |
|
|
|
2609 |
Model 1: Speed ~ Expt
|
|
|
2610 |
Model 2: Speed ~ Run + Expt
|
|
|
2611 |
Res.Df RSS Df Sum of Sq F Pr(>F)
|
|
|
2612 |
1 95 523510
|
|
|
2613 |
2 76 410166 19 113344 1.1053 0.3632
|
|
|
2614 |
>
|
|
|
2615 |
>
|
|
|
2616 |
>
|
|
|
2617 |
> cleanEx()
|
|
|
2618 |
> nameEx("mtcars")
|
|
|
2619 |
> ### * mtcars
|
|
|
2620 |
>
|
|
|
2621 |
> flush(stderr()); flush(stdout())
|
|
|
2622 |
>
|
|
|
2623 |
> ### Name: mtcars
|
|
|
2624 |
> ### Title: Motor Trend Car Road Tests
|
|
|
2625 |
> ### Aliases: mtcars
|
|
|
2626 |
> ### Keywords: datasets
|
|
|
2627 |
>
|
|
|
2628 |
> ### ** Examples
|
|
|
2629 |
>
|
|
|
2630 |
> require(graphics)
|
|
|
2631 |
> pairs(mtcars, main = "mtcars data")
|
|
|
2632 |
> coplot(mpg ~ disp | as.factor(cyl), data = mtcars,
|
|
|
2633 |
+ panel = panel.smooth, rows = 1)
|
|
|
2634 |
>
|
|
|
2635 |
>
|
|
|
2636 |
>
|
|
|
2637 |
> cleanEx()
|
|
|
2638 |
> nameEx("nhtemp")
|
|
|
2639 |
> ### * nhtemp
|
|
|
2640 |
>
|
|
|
2641 |
> flush(stderr()); flush(stdout())
|
|
|
2642 |
>
|
|
|
2643 |
> ### Name: nhtemp
|
|
|
2644 |
> ### Title: Average Yearly Temperatures in New Haven
|
|
|
2645 |
> ### Aliases: nhtemp
|
|
|
2646 |
> ### Keywords: datasets
|
|
|
2647 |
>
|
|
|
2648 |
> ### ** Examples
|
|
|
2649 |
>
|
|
|
2650 |
> require(stats); require(graphics)
|
|
|
2651 |
> plot(nhtemp, main = "nhtemp data",
|
|
|
2652 |
+ ylab = "Mean annual temperature in New Haven, CT (deg. F)")
|
|
|
2653 |
>
|
|
|
2654 |
>
|
|
|
2655 |
>
|
|
|
2656 |
> cleanEx()
|
|
|
2657 |
> nameEx("nottem")
|
|
|
2658 |
> ### * nottem
|
|
|
2659 |
>
|
|
|
2660 |
> flush(stderr()); flush(stdout())
|
|
|
2661 |
>
|
|
|
2662 |
> ### Name: nottem
|
|
|
2663 |
> ### Title: Average Monthly Temperatures at Nottingham, 1920-1939
|
|
|
2664 |
> ### Aliases: nottem
|
|
|
2665 |
> ### Keywords: datasets
|
|
|
2666 |
>
|
|
|
2667 |
> ### ** Examples
|
|
|
2668 |
>
|
|
|
2669 |
> ## Not run:
|
|
|
2670 |
> ##D require(stats); require(graphics)
|
| 61157 |
ripley |
2671 |
> ##D nott <- window(nottem, end = c(1936,12))
|
| 61169 |
ripley |
2672 |
> ##D fit <- arima(nott, order = c(1,0,0), list(order = c(2,1,0), period = 12))
|
| 61157 |
ripley |
2673 |
> ##D nott.fore <- predict(fit, n.ahead = 36)
|
| 56186 |
murdoch |
2674 |
> ##D ts.plot(nott, nott.fore$pred, nott.fore$pred+2*nott.fore$se,
|
| 61157 |
ripley |
2675 |
> ##D nott.fore$pred-2*nott.fore$se, gpars = list(col = c(1,1,4,4)))
|
| 56186 |
murdoch |
2676 |
> ## End(Not run)
|
|
|
2677 |
>
|
|
|
2678 |
>
|
|
|
2679 |
> cleanEx()
|
|
|
2680 |
> nameEx("occupationalStatus")
|
|
|
2681 |
> ### * occupationalStatus
|
|
|
2682 |
>
|
|
|
2683 |
> flush(stderr()); flush(stdout())
|
|
|
2684 |
>
|
|
|
2685 |
> ### Name: occupationalStatus
|
|
|
2686 |
> ### Title: Occupational Status of Fathers and their Sons
|
|
|
2687 |
> ### Aliases: occupationalStatus
|
|
|
2688 |
> ### Keywords: datasets
|
|
|
2689 |
>
|
|
|
2690 |
> ### ** Examples
|
|
|
2691 |
>
|
|
|
2692 |
> require(stats); require(graphics)
|
|
|
2693 |
>
|
|
|
2694 |
> plot(occupationalStatus)
|
|
|
2695 |
>
|
|
|
2696 |
> ## Fit a uniform association model separating diagonal effects
|
|
|
2697 |
> Diag <- as.factor(diag(1:8))
|
|
|
2698 |
> Rscore <- scale(as.numeric(row(occupationalStatus)), scale = FALSE)
|
|
|
2699 |
> Cscore <- scale(as.numeric(col(occupationalStatus)), scale = FALSE)
|
|
|
2700 |
> modUnif <- glm(Freq ~ origin + destination + Diag + Rscore:Cscore,
|
|
|
2701 |
+ family = poisson, data = occupationalStatus)
|
|
|
2702 |
>
|
|
|
2703 |
> summary(modUnif)
|
|
|
2704 |
|
|
|
2705 |
Call:
|
|
|
2706 |
glm(formula = Freq ~ origin + destination + Diag + Rscore:Cscore,
|
|
|
2707 |
family = poisson, data = occupationalStatus)
|
|
|
2708 |
|
|
|
2709 |
Deviance Residuals:
|
|
|
2710 |
Min 1Q Median 3Q Max
|
|
|
2711 |
-2.6521 -0.6267 0.0000 0.5913 2.0964
|
|
|
2712 |
|
|
|
2713 |
Coefficients:
|
|
|
2714 |
Estimate Std. Error z value Pr(>|z|)
|
|
|
2715 |
(Intercept) 0.568592 0.183358 3.101 0.001929 **
|
|
|
2716 |
origin2 0.431314 0.149415 2.887 0.003893 **
|
|
|
2717 |
origin3 1.461862 0.131141 11.147 < 2e-16 ***
|
|
|
2718 |
origin4 1.788731 0.126588 14.130 < 2e-16 ***
|
|
|
2719 |
origin5 0.441011 0.144844 3.045 0.002329 **
|
|
|
2720 |
origin6 2.491735 0.121219 20.556 < 2e-16 ***
|
|
|
2721 |
origin7 1.127536 0.129032 8.738 < 2e-16 ***
|
|
|
2722 |
origin8 0.796445 0.131863 6.040 1.54e-09 ***
|
|
|
2723 |
destination2 0.873202 0.166902 5.232 1.68e-07 ***
|
|
|
2724 |
destination3 1.813992 0.153861 11.790 < 2e-16 ***
|
|
|
2725 |
destination4 2.082515 0.150887 13.802 < 2e-16 ***
|
|
|
2726 |
destination5 1.366383 0.155590 8.782 < 2e-16 ***
|
|
|
2727 |
destination6 2.816369 0.146054 19.283 < 2e-16 ***
|
|
|
2728 |
destination7 1.903918 0.147810 12.881 < 2e-16 ***
|
|
|
2729 |
destination8 1.398585 0.151658 9.222 < 2e-16 ***
|
|
|
2730 |
Diag1 1.665495 0.237383 7.016 2.28e-12 ***
|
|
|
2731 |
Diag2 0.959681 0.212122 4.524 6.06e-06 ***
|
|
|
2732 |
Diag3 0.021750 0.156530 0.139 0.889490
|
|
|
2733 |
Diag4 0.226399 0.124364 1.820 0.068689 .
|
|
|
2734 |
Diag5 0.808646 0.229754 3.520 0.000432 ***
|
|
|
2735 |
Diag6 0.132277 0.077191 1.714 0.086597 .
|
|
|
2736 |
Diag7 0.506709 0.115936 4.371 1.24e-05 ***
|
|
|
2737 |
Diag8 0.221880 0.134803 1.646 0.099771 .
|
|
|
2738 |
Rscore:Cscore 0.136974 0.007489 18.289 < 2e-16 ***
|
|
|
2739 |
---
|
| 61435 |
ripley |
2740 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2741 |
|
|
|
2742 |
(Dispersion parameter for poisson family taken to be 1)
|
|
|
2743 |
|
|
|
2744 |
Null deviance: 4679.004 on 63 degrees of freedom
|
|
|
2745 |
Residual deviance: 58.436 on 40 degrees of freedom
|
|
|
2746 |
AIC: 428.78
|
|
|
2747 |
|
|
|
2748 |
Number of Fisher Scoring iterations: 4
|
|
|
2749 |
|
|
|
2750 |
> plot(modUnif) # 4 plots, with warning about h_ii ~= 1
|
| 60841 |
ripley |
2751 |
Warning: not plotting observations with leverage one:
|
| 56186 |
murdoch |
2752 |
1, 10, 19, 28, 37, 46, 55, 64
|
| 60841 |
ripley |
2753 |
Warning: not plotting observations with leverage one:
|
| 56186 |
murdoch |
2754 |
1, 10, 19, 28, 37, 46, 55, 64
|
|
|
2755 |
>
|
|
|
2756 |
>
|
|
|
2757 |
>
|
|
|
2758 |
> cleanEx()
|
|
|
2759 |
> nameEx("precip")
|
|
|
2760 |
> ### * precip
|
|
|
2761 |
>
|
|
|
2762 |
> flush(stderr()); flush(stdout())
|
|
|
2763 |
>
|
|
|
2764 |
> ### Name: precip
|
|
|
2765 |
> ### Title: Annual Precipitation in US Cities
|
|
|
2766 |
> ### Aliases: precip
|
|
|
2767 |
> ### Keywords: datasets
|
|
|
2768 |
>
|
|
|
2769 |
> ### ** Examples
|
|
|
2770 |
>
|
|
|
2771 |
> require(graphics)
|
|
|
2772 |
> dotchart(precip[order(precip)], main = "precip data")
|
|
|
2773 |
> title(sub = "Average annual precipitation (in.)")
|
|
|
2774 |
>
|
|
|
2775 |
>
|
|
|
2776 |
>
|
|
|
2777 |
> cleanEx()
|
|
|
2778 |
> nameEx("presidents")
|
|
|
2779 |
> ### * presidents
|
|
|
2780 |
>
|
|
|
2781 |
> flush(stderr()); flush(stdout())
|
|
|
2782 |
>
|
|
|
2783 |
> ### Name: presidents
|
|
|
2784 |
> ### Title: Quarterly Approval Ratings of US Presidents
|
|
|
2785 |
> ### Aliases: presidents
|
|
|
2786 |
> ### Keywords: datasets
|
|
|
2787 |
>
|
|
|
2788 |
> ### ** Examples
|
|
|
2789 |
>
|
|
|
2790 |
> require(stats); require(graphics)
|
|
|
2791 |
> plot(presidents, las = 1, ylab = "Approval rating (%)",
|
|
|
2792 |
+ main = "presidents data")
|
|
|
2793 |
>
|
|
|
2794 |
>
|
|
|
2795 |
>
|
|
|
2796 |
> cleanEx()
|
|
|
2797 |
> nameEx("pressure")
|
|
|
2798 |
> ### * pressure
|
|
|
2799 |
>
|
|
|
2800 |
> flush(stderr()); flush(stdout())
|
|
|
2801 |
>
|
|
|
2802 |
> ### Name: pressure
|
|
|
2803 |
> ### Title: Vapor Pressure of Mercury as a Function of Temperature
|
|
|
2804 |
> ### Aliases: pressure
|
|
|
2805 |
> ### Keywords: datasets
|
|
|
2806 |
>
|
|
|
2807 |
> ### ** Examples
|
|
|
2808 |
>
|
|
|
2809 |
> require(graphics)
|
|
|
2810 |
> plot(pressure, xlab = "Temperature (deg C)",
|
|
|
2811 |
+ ylab = "Pressure (mm of Hg)",
|
|
|
2812 |
+ main = "pressure data: Vapor Pressure of Mercury")
|
|
|
2813 |
> plot(pressure, xlab = "Temperature (deg C)", log = "y",
|
|
|
2814 |
+ ylab = "Pressure (mm of Hg)",
|
|
|
2815 |
+ main = "pressure data: Vapor Pressure of Mercury")
|
|
|
2816 |
>
|
|
|
2817 |
>
|
|
|
2818 |
>
|
|
|
2819 |
> cleanEx()
|
|
|
2820 |
> nameEx("quakes")
|
|
|
2821 |
> ### * quakes
|
|
|
2822 |
>
|
|
|
2823 |
> flush(stderr()); flush(stdout())
|
|
|
2824 |
>
|
|
|
2825 |
> ### Name: quakes
|
|
|
2826 |
> ### Title: Locations of Earthquakes off Fiji
|
|
|
2827 |
> ### Aliases: quakes
|
|
|
2828 |
> ### Keywords: datasets
|
|
|
2829 |
>
|
|
|
2830 |
> ### ** Examples
|
|
|
2831 |
>
|
|
|
2832 |
> require(graphics)
|
| 61157 |
ripley |
2833 |
> pairs(quakes, main = "Fiji Earthquakes, N = 1000", cex.main = 1.2, pch = ".")
|
| 56186 |
murdoch |
2834 |
>
|
|
|
2835 |
>
|
|
|
2836 |
>
|
|
|
2837 |
> cleanEx()
|
|
|
2838 |
> nameEx("randu")
|
|
|
2839 |
> ### * randu
|
|
|
2840 |
>
|
|
|
2841 |
> flush(stderr()); flush(stdout())
|
|
|
2842 |
>
|
|
|
2843 |
> ### Name: randu
|
|
|
2844 |
> ### Title: Random Numbers from Congruential Generator RANDU
|
|
|
2845 |
> ### Aliases: randu
|
|
|
2846 |
> ### Keywords: datasets
|
|
|
2847 |
>
|
|
|
2848 |
> ### ** Examples
|
|
|
2849 |
>
|
|
|
2850 |
> ## Not run:
|
|
|
2851 |
> ##D ## We could re-generate the dataset by the following R code
|
|
|
2852 |
> ##D seed <- as.double(1)
|
|
|
2853 |
> ##D RANDU <- function() {
|
|
|
2854 |
> ##D seed <<- ((2^16 + 3) * seed) %% (2^31)
|
|
|
2855 |
> ##D seed/(2^31)
|
|
|
2856 |
> ##D }
|
|
|
2857 |
> ##D for(i in 1:400) {
|
|
|
2858 |
> ##D U <- c(RANDU(), RANDU(), RANDU(), RANDU(), RANDU())
|
|
|
2859 |
> ##D print(round(U[1:3], 6))
|
|
|
2860 |
> ##D }
|
|
|
2861 |
> ## End(Not run)
|
|
|
2862 |
>
|
|
|
2863 |
>
|
|
|
2864 |
>
|
|
|
2865 |
> cleanEx()
|
|
|
2866 |
> nameEx("sleep")
|
|
|
2867 |
> ### * sleep
|
|
|
2868 |
>
|
|
|
2869 |
> flush(stderr()); flush(stdout())
|
|
|
2870 |
>
|
|
|
2871 |
> ### Name: sleep
|
|
|
2872 |
> ### Title: Student's Sleep Data
|
|
|
2873 |
> ### Aliases: sleep
|
|
|
2874 |
> ### Keywords: datasets
|
|
|
2875 |
>
|
|
|
2876 |
> ### ** Examples
|
|
|
2877 |
>
|
|
|
2878 |
> require(stats)
|
|
|
2879 |
> ## Student's paired t-test
|
|
|
2880 |
> with(sleep,
|
|
|
2881 |
+ t.test(extra[group == 1],
|
|
|
2882 |
+ extra[group == 2], paired = TRUE))
|
|
|
2883 |
|
|
|
2884 |
Paired t-test
|
|
|
2885 |
|
| 61430 |
ripley |
2886 |
data: extra[group == 1] and extra[group == 2]
|
| 56186 |
murdoch |
2887 |
t = -4.0621, df = 9, p-value = 0.002833
|
| 61430 |
ripley |
2888 |
alternative hypothesis: true difference in means is not equal to 0
|
| 56186 |
murdoch |
2889 |
95 percent confidence interval:
|
| 61430 |
ripley |
2890 |
-2.4598858 -0.7001142
|
| 56186 |
murdoch |
2891 |
sample estimates:
|
|
|
2892 |
mean of the differences
|
|
|
2893 |
-1.58
|
|
|
2894 |
|
|
|
2895 |
>
|
| 59811 |
maechler |
2896 |
> ## The sleep *prolongations*
|
|
|
2897 |
> sleep1 <- with(sleep, extra[group == 2] - extra[group == 1])
|
|
|
2898 |
> summary(sleep1)
|
|
|
2899 |
Min. 1st Qu. Median Mean 3rd Qu. Max.
|
|
|
2900 |
0.00 1.05 1.30 1.58 1.70 4.60
|
|
|
2901 |
> stripchart(sleep1, method = "stack", xlab = "hours",
|
|
|
2902 |
+ main = "Sleep prolongation (n = 10)")
|
| 61157 |
ripley |
2903 |
> boxplot(sleep1, horizontal = TRUE, add = TRUE,
|
|
|
2904 |
+ at = .6, pars = list(boxwex = 0.5, staplewex = 0.25))
|
| 56186 |
murdoch |
2905 |
>
|
|
|
2906 |
>
|
| 59811 |
maechler |
2907 |
>
|
| 56186 |
murdoch |
2908 |
> cleanEx()
|
|
|
2909 |
> nameEx("stackloss")
|
|
|
2910 |
> ### * stackloss
|
|
|
2911 |
>
|
|
|
2912 |
> flush(stderr()); flush(stdout())
|
|
|
2913 |
>
|
|
|
2914 |
> ### Name: stackloss
|
|
|
2915 |
> ### Title: Brownlee's Stack Loss Plant Data
|
|
|
2916 |
> ### Aliases: stackloss stack.loss stack.x
|
|
|
2917 |
> ### Keywords: datasets
|
|
|
2918 |
>
|
|
|
2919 |
> ### ** Examples
|
|
|
2920 |
>
|
|
|
2921 |
> require(stats)
|
|
|
2922 |
> summary(lm.stack <- lm(stack.loss ~ stack.x))
|
|
|
2923 |
|
|
|
2924 |
Call:
|
|
|
2925 |
lm(formula = stack.loss ~ stack.x)
|
|
|
2926 |
|
|
|
2927 |
Residuals:
|
|
|
2928 |
Min 1Q Median 3Q Max
|
|
|
2929 |
-7.2377 -1.7117 -0.4551 2.3614 5.6978
|
|
|
2930 |
|
|
|
2931 |
Coefficients:
|
|
|
2932 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
2933 |
(Intercept) -39.9197 11.8960 -3.356 0.00375 **
|
|
|
2934 |
stack.xAir.Flow 0.7156 0.1349 5.307 5.8e-05 ***
|
|
|
2935 |
stack.xWater.Temp 1.2953 0.3680 3.520 0.00263 **
|
|
|
2936 |
stack.xAcid.Conc. -0.1521 0.1563 -0.973 0.34405
|
|
|
2937 |
---
|
| 61435 |
ripley |
2938 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
2939 |
|
|
|
2940 |
Residual standard error: 3.243 on 17 degrees of freedom
|
| 61435 |
ripley |
2941 |
Multiple R-squared: 0.9136, Adjusted R-squared: 0.8983
|
|
|
2942 |
F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09
|
| 56186 |
murdoch |
2943 |
|
|
|
2944 |
>
|
|
|
2945 |
>
|
|
|
2946 |
>
|
|
|
2947 |
> cleanEx()
|
|
|
2948 |
> nameEx("sunspot.month")
|
|
|
2949 |
> ### * sunspot.month
|
|
|
2950 |
>
|
|
|
2951 |
> flush(stderr()); flush(stdout())
|
|
|
2952 |
>
|
|
|
2953 |
> ### Name: sunspot.month
|
|
|
2954 |
> ### Title: Monthly Sunspot Data, 1749-1997
|
|
|
2955 |
> ### Aliases: sunspot.month
|
|
|
2956 |
> ### Keywords: datasets
|
|
|
2957 |
>
|
|
|
2958 |
> ### ** Examples
|
|
|
2959 |
>
|
|
|
2960 |
> require(stats); require(graphics)
|
| 61435 |
ripley |
2961 |
> ## Compare the monthly series
|
| 56186 |
murdoch |
2962 |
> plot (sunspot.month, main = "sunspot.month [stats]", col = 2)
|
|
|
2963 |
> lines(sunspots) # "very barely" see something
|
|
|
2964 |
>
|
|
|
2965 |
> ## Now look at the difference :
|
|
|
2966 |
> all(tsp(sunspots) [c(1,3)] ==
|
|
|
2967 |
+ tsp(sunspot.month)[c(1,3)]) ## Start & Periodicity are the same
|
|
|
2968 |
[1] TRUE
|
|
|
2969 |
> n1 <- length(sunspots)
|
|
|
2970 |
> table(eq <- sunspots == sunspot.month[1:n1]) #> 132 are different !
|
|
|
2971 |
|
|
|
2972 |
FALSE TRUE
|
|
|
2973 |
132 2688
|
| 61435 |
ripley |
2974 |
> i <- which(!eq)
|
| 56186 |
murdoch |
2975 |
> rug(time(eq)[i])
|
|
|
2976 |
> s1 <- sunspots[i] ; s2 <- sunspot.month[i]
|
|
|
2977 |
> cbind(i = i, sunspots = s1, ss.month = s2,
|
|
|
2978 |
+ perc.diff = round(100*2*abs(s1-s2)/(s1+s2), 1))
|
|
|
2979 |
i sunspots ss.month perc.diff
|
|
|
2980 |
[1,] 55 22.2 22.0 0.9
|
|
|
2981 |
[2,] 838 31.7 31.6 0.3
|
|
|
2982 |
[3,] 841 32.5 32.8 0.9
|
|
|
2983 |
[4,] 862 9.0 8.9 1.1
|
|
|
2984 |
[5,] 864 9.7 9.1 6.4
|
|
|
2985 |
[6,] 866 4.3 4.2 2.4
|
|
|
2986 |
[7,] 876 0.0 0.2 200.0
|
|
|
2987 |
[8,] 901 21.6 21.7 0.5
|
|
|
2988 |
[9,] 917 15.4 15.5 0.6
|
|
|
2989 |
[10,] 920 25.4 25.7 1.2
|
|
|
2990 |
[11,] 943 42.9 42.3 1.4
|
|
|
2991 |
[12,] 946 57.2 56.1 1.9
|
|
|
2992 |
[13,] 955 54.3 54.2 0.2
|
|
|
2993 |
[14,] 960 46.6 46.9 0.6
|
|
|
2994 |
[15,] 965 67.5 67.4 0.1
|
|
|
2995 |
[16,] 968 78.3 77.6 0.9
|
|
|
2996 |
[17,] 976 107.1 106.3 0.7
|
|
|
2997 |
[18,] 988 54.6 54.5 0.2
|
|
|
2998 |
[19,] 992 54.9 55.0 0.2
|
|
|
2999 |
[20,] 994 46.2 46.3 0.2
|
|
|
3000 |
[21,] 998 55.5 55.6 0.2
|
|
|
3001 |
[22,] 1003 13.9 14.0 0.7
|
|
|
3002 |
[23,] 1047 98.1 98.2 0.1
|
|
|
3003 |
[24,] 1061 111.3 111.7 0.4
|
|
|
3004 |
[25,] 1081 107.6 105.6 1.9
|
|
|
3005 |
[26,] 1087 84.7 84.8 0.1
|
|
|
3006 |
[27,] 1090 90.8 90.9 0.1
|
|
|
3007 |
[28,] 1092 63.6 63.7 0.2
|
|
|
3008 |
[29,] 1095 55.5 67.8 20.0
|
|
|
3009 |
[30,] 1102 49.8 55.0 9.9
|
|
|
3010 |
[31,] 1105 24.0 24.1 0.4
|
|
|
3011 |
[32,] 1108 42.6 40.2 5.8
|
|
|
3012 |
[33,] 1109 67.4 67.5 0.1
|
|
|
3013 |
[34,] 1113 35.1 36.5 3.9
|
|
|
3014 |
[35,] 1124 26.5 26.6 0.4
|
|
|
3015 |
[36,] 1125 18.5 18.4 0.5
|
|
|
3016 |
[37,] 1132 8.8 9.5 7.7
|
|
|
3017 |
[38,] 1145 12.0 11.6 3.4
|
|
|
3018 |
[39,] 1149 6.9 7.0 1.4
|
|
|
3019 |
[40,] 1156 56.9 57.0 0.2
|
|
|
3020 |
[41,] 1168 69.2 69.3 0.1
|
|
|
3021 |
[42,] 1185 161.2 160.9 0.2
|
|
|
3022 |
[43,] 1191 108.9 108.6 0.3
|
|
|
3023 |
[44,] 1194 123.8 129.0 4.1
|
|
|
3024 |
[45,] 1196 132.5 132.6 0.1
|
|
|
3025 |
[46,] 1200 159.9 159.5 0.3
|
|
|
3026 |
[47,] 1201 156.7 157.0 0.2
|
|
|
3027 |
[48,] 1203 96.5 96.2 0.3
|
|
|
3028 |
[49,] 1206 81.2 81.1 0.1
|
|
|
3029 |
[50,] 1208 61.3 67.7 9.9
|
|
|
3030 |
[51,] 1211 99.7 99.0 0.7
|
|
|
3031 |
[52,] 1224 60.0 61.0 1.7
|
|
|
3032 |
[53,] 1235 50.9 51.0 0.2
|
|
|
3033 |
[54,] 1238 67.5 66.4 1.6
|
|
|
3034 |
[55,] 1243 42.0 42.1 0.2
|
|
|
3035 |
[56,] 1264 26.4 26.5 0.4
|
|
|
3036 |
[57,] 1270 12.7 12.6 0.8
|
|
|
3037 |
[58,] 1272 21.4 21.6 0.9
|
|
|
3038 |
[59,] 1282 9.7 9.6 1.0
|
|
|
3039 |
[60,] 1283 4.3 4.2 2.4
|
|
|
3040 |
[61,] 1290 5.0 5.2 3.9
|
|
|
3041 |
[62,] 1301 29.2 28.6 2.1
|
|
|
3042 |
[63,] 1333 81.5 82.4 1.1
|
|
|
3043 |
[64,] 1334 88.0 88.3 0.3
|
|
|
3044 |
[65,] 1346 77.8 77.7 0.1
|
|
|
3045 |
[66,] 1350 87.8 88.1 0.3
|
|
|
3046 |
[67,] 1366 42.0 41.9 0.2
|
|
|
3047 |
[68,] 1424 4.9 4.8 2.1
|
|
|
3048 |
[69,] 1427 9.3 9.6 3.2
|
|
|
3049 |
[70,] 1430 15.8 15.7 0.6
|
|
|
3050 |
[71,] 1435 28.6 29.0 1.4
|
|
|
3051 |
[72,] 1437 43.8 47.2 7.5
|
|
|
3052 |
[73,] 1438 61.7 61.6 0.2
|
|
|
3053 |
[74,] 1442 59.3 59.9 1.0
|
|
|
3054 |
[75,] 1445 104.0 103.9 0.1
|
|
|
3055 |
[76,] 1450 59.4 59.3 0.2
|
|
|
3056 |
[77,] 1451 77.4 78.1 0.9
|
|
|
3057 |
[78,] 1455 159.4 157.6 1.1
|
|
|
3058 |
[79,] 1472 110.0 110.1 0.1
|
|
|
3059 |
[80,] 1476 90.3 90.4 0.1
|
|
|
3060 |
[81,] 1486 103.5 102.6 0.9
|
|
|
3061 |
[82,] 1497 47.5 47.1 0.8
|
|
|
3062 |
[83,] 1498 47.4 47.1 0.6
|
|
|
3063 |
[84,] 1514 22.2 21.5 3.2
|
|
|
3064 |
[85,] 1527 31.2 30.6 1.9
|
|
|
3065 |
[86,] 1539 11.7 11.9 1.7
|
|
|
3066 |
[87,] 1541 21.2 21.6 1.9
|
|
|
3067 |
[88,] 1542 13.4 14.2 5.8
|
|
|
3068 |
[89,] 1543 5.9 6.0 1.7
|
|
|
3069 |
[90,] 1545 16.4 16.9 3.0
|
|
|
3070 |
[91,] 1547 14.5 14.2 2.1
|
|
|
3071 |
[92,] 1548 2.3 2.2 4.4
|
|
|
3072 |
[93,] 1550 6.0 6.6 9.5
|
|
|
3073 |
[94,] 1553 5.8 5.9 1.7
|
|
|
3074 |
[95,] 1561 0.8 1.0 22.2
|
|
|
3075 |
[96,] 1571 12.9 13.1 1.5
|
|
|
3076 |
[97,] 1572 7.2 7.3 1.4
|
|
|
3077 |
[98,] 1574 27.5 27.2 1.1
|
|
|
3078 |
[99,] 1575 19.5 19.3 1.0
|
|
|
3079 |
[100,] 1576 19.3 19.5 1.0
|
|
|
3080 |
[101,] 1588 51.7 51.6 0.2
|
|
|
3081 |
[102,] 1592 58.0 58.4 0.7
|
|
|
3082 |
[103,] 1594 64.0 64.4 0.6
|
|
|
3083 |
[104,] 1598 69.3 69.5 0.3
|
|
|
3084 |
[105,] 1599 67.5 66.8 1.0
|
|
|
3085 |
[106,] 1613 32.1 31.5 1.9
|
|
|
3086 |
[107,] 1614 76.5 76.3 0.3
|
|
|
3087 |
[108,] 1623 86.8 87.5 0.8
|
|
|
3088 |
[109,] 1643 33.3 30.9 7.5
|
|
|
3089 |
[110,] 1656 12.4 13.0 4.7
|
|
|
3090 |
[111,] 1683 7.0 6.7 4.4
|
|
|
3091 |
[112,] 1687 9.7 9.4 3.1
|
|
|
3092 |
[113,] 1712 33.2 33.0 0.6
|
|
|
3093 |
[114,] 1716 32.3 32.5 0.6
|
|
|
3094 |
[115,] 1723 76.8 76.5 0.4
|
|
|
3095 |
[116,] 1734 88.2 89.9 1.9
|
|
|
3096 |
[117,] 1735 88.8 88.6 0.2
|
|
|
3097 |
[118,] 1738 79.7 80.0 0.4
|
|
|
3098 |
[119,] 1774 28.4 28.7 1.1
|
|
|
3099 |
[120,] 1837 5.2 5.5 5.6
|
|
|
3100 |
[121,] 2126 70.0 69.9 0.1
|
|
|
3101 |
[122,] 2153 76.9 77.0 0.1
|
|
|
3102 |
[123,] 2162 64.1 62.8 2.0
|
|
|
3103 |
[124,] 2174 49.2 49.9 1.4
|
|
|
3104 |
[125,] 2233 18.9 18.6 1.6
|
|
|
3105 |
[126,] 2315 38.3 38.4 0.3
|
|
|
3106 |
[127,] 2378 113.4 133.4 16.2
|
|
|
3107 |
[128,] 2427 59.9 55.9 6.9
|
|
|
3108 |
[129,] 2552 55.9 55.8 0.2
|
|
|
3109 |
[130,] 2556 40.0 39.9 0.3
|
|
|
3110 |
[131,] 2694 39.5 37.5 5.2
|
|
|
3111 |
[132,] 2790 90.0 90.9 1.0
|
|
|
3112 |
>
|
|
|
3113 |
>
|
|
|
3114 |
>
|
|
|
3115 |
> cleanEx()
|
|
|
3116 |
> nameEx("sunspots")
|
|
|
3117 |
> ### * sunspots
|
|
|
3118 |
>
|
|
|
3119 |
> flush(stderr()); flush(stdout())
|
|
|
3120 |
>
|
|
|
3121 |
> ### Name: sunspots
|
|
|
3122 |
> ### Title: Monthly Sunspot Numbers, 1749-1983
|
|
|
3123 |
> ### Aliases: sunspots
|
|
|
3124 |
> ### Keywords: datasets
|
|
|
3125 |
>
|
|
|
3126 |
> ### ** Examples
|
|
|
3127 |
>
|
|
|
3128 |
> require(graphics)
|
|
|
3129 |
> plot(sunspots, main = "sunspots data", xlab = "Year",
|
|
|
3130 |
+ ylab = "Monthly sunspot numbers")
|
|
|
3131 |
>
|
|
|
3132 |
>
|
|
|
3133 |
>
|
|
|
3134 |
> cleanEx()
|
|
|
3135 |
> nameEx("swiss")
|
|
|
3136 |
> ### * swiss
|
|
|
3137 |
>
|
|
|
3138 |
> flush(stderr()); flush(stdout())
|
|
|
3139 |
>
|
|
|
3140 |
> ### Name: swiss
|
|
|
3141 |
> ### Title: Swiss Fertility and Socioeconomic Indicators (1888) Data
|
|
|
3142 |
> ### Aliases: swiss
|
|
|
3143 |
> ### Keywords: datasets
|
|
|
3144 |
>
|
|
|
3145 |
> ### ** Examples
|
|
|
3146 |
>
|
|
|
3147 |
> require(stats); require(graphics)
|
|
|
3148 |
> pairs(swiss, panel = panel.smooth, main = "swiss data",
|
|
|
3149 |
+ col = 3 + (swiss$Catholic > 50))
|
|
|
3150 |
> summary(lm(Fertility ~ . , data = swiss))
|
|
|
3151 |
|
|
|
3152 |
Call:
|
|
|
3153 |
lm(formula = Fertility ~ ., data = swiss)
|
|
|
3154 |
|
|
|
3155 |
Residuals:
|
|
|
3156 |
Min 1Q Median 3Q Max
|
|
|
3157 |
-15.2743 -5.2617 0.5032 4.1198 15.3213
|
|
|
3158 |
|
|
|
3159 |
Coefficients:
|
|
|
3160 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3161 |
(Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
|
|
|
3162 |
Agriculture -0.17211 0.07030 -2.448 0.01873 *
|
|
|
3163 |
Examination -0.25801 0.25388 -1.016 0.31546
|
|
|
3164 |
Education -0.87094 0.18303 -4.758 2.43e-05 ***
|
|
|
3165 |
Catholic 0.10412 0.03526 2.953 0.00519 **
|
|
|
3166 |
Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
|
|
|
3167 |
---
|
| 61435 |
ripley |
3168 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3169 |
|
|
|
3170 |
Residual standard error: 7.165 on 41 degrees of freedom
|
| 61435 |
ripley |
3171 |
Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
|
|
|
3172 |
F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10
|
| 56186 |
murdoch |
3173 |
|
|
|
3174 |
>
|
|
|
3175 |
>
|
|
|
3176 |
>
|
|
|
3177 |
> cleanEx()
|
|
|
3178 |
> nameEx("trees")
|
|
|
3179 |
> ### * trees
|
|
|
3180 |
>
|
|
|
3181 |
> flush(stderr()); flush(stdout())
|
|
|
3182 |
>
|
|
|
3183 |
> ### Name: trees
|
|
|
3184 |
> ### Title: Girth, Height and Volume for Black Cherry Trees
|
|
|
3185 |
> ### Aliases: trees
|
|
|
3186 |
> ### Keywords: datasets
|
|
|
3187 |
>
|
|
|
3188 |
> ### ** Examples
|
|
|
3189 |
>
|
|
|
3190 |
> require(stats); require(graphics)
|
|
|
3191 |
> pairs(trees, panel = panel.smooth, main = "trees data")
|
|
|
3192 |
> plot(Volume ~ Girth, data = trees, log = "xy")
|
|
|
3193 |
> coplot(log(Volume) ~ log(Girth) | Height, data = trees,
|
|
|
3194 |
+ panel = panel.smooth)
|
|
|
3195 |
> summary(fm1 <- lm(log(Volume) ~ log(Girth), data = trees))
|
|
|
3196 |
|
|
|
3197 |
Call:
|
|
|
3198 |
lm(formula = log(Volume) ~ log(Girth), data = trees)
|
|
|
3199 |
|
|
|
3200 |
Residuals:
|
|
|
3201 |
Min 1Q Median 3Q Max
|
|
|
3202 |
-0.205999 -0.068702 0.001011 0.072585 0.247963
|
|
|
3203 |
|
|
|
3204 |
Coefficients:
|
|
|
3205 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3206 |
(Intercept) -2.35332 0.23066 -10.20 4.18e-11 ***
|
|
|
3207 |
log(Girth) 2.19997 0.08983 24.49 < 2e-16 ***
|
|
|
3208 |
---
|
| 61435 |
ripley |
3209 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3210 |
|
|
|
3211 |
Residual standard error: 0.115 on 29 degrees of freedom
|
| 61435 |
ripley |
3212 |
Multiple R-squared: 0.9539, Adjusted R-squared: 0.9523
|
|
|
3213 |
F-statistic: 599.7 on 1 and 29 DF, p-value: < 2.2e-16
|
| 56186 |
murdoch |
3214 |
|
|
|
3215 |
> summary(fm2 <- update(fm1, ~ . + log(Height), data = trees))
|
|
|
3216 |
|
|
|
3217 |
Call:
|
|
|
3218 |
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
|
|
|
3219 |
|
|
|
3220 |
Residuals:
|
|
|
3221 |
Min 1Q Median 3Q Max
|
|
|
3222 |
-0.168561 -0.048488 0.002431 0.063637 0.129223
|
|
|
3223 |
|
|
|
3224 |
Coefficients:
|
|
|
3225 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3226 |
(Intercept) -6.63162 0.79979 -8.292 5.06e-09 ***
|
|
|
3227 |
log(Girth) 1.98265 0.07501 26.432 < 2e-16 ***
|
|
|
3228 |
log(Height) 1.11712 0.20444 5.464 7.81e-06 ***
|
|
|
3229 |
---
|
| 61435 |
ripley |
3230 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3231 |
|
|
|
3232 |
Residual standard error: 0.08139 on 28 degrees of freedom
|
| 61435 |
ripley |
3233 |
Multiple R-squared: 0.9777, Adjusted R-squared: 0.9761
|
|
|
3234 |
F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16
|
| 56186 |
murdoch |
3235 |
|
|
|
3236 |
> step(fm2)
|
|
|
3237 |
Start: AIC=-152.69
|
|
|
3238 |
log(Volume) ~ log(Girth) + log(Height)
|
|
|
3239 |
|
| 57809 |
ripley |
3240 |
Df Sum of Sq RSS AIC
|
|
|
3241 |
<none> 0.1855 -152.685
|
|
|
3242 |
- log(Height) 1 0.1978 0.3832 -132.185
|
|
|
3243 |
- log(Girth) 1 4.6275 4.8130 -53.743
|
| 56186 |
murdoch |
3244 |
|
|
|
3245 |
Call:
|
|
|
3246 |
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
|
|
|
3247 |
|
|
|
3248 |
Coefficients:
|
|
|
3249 |
(Intercept) log(Girth) log(Height)
|
|
|
3250 |
-6.632 1.983 1.117
|
|
|
3251 |
|
|
|
3252 |
> ## i.e., Volume ~= c * Height * Girth^2 seems reasonable
|
|
|
3253 |
>
|
|
|
3254 |
>
|
|
|
3255 |
>
|
|
|
3256 |
> cleanEx()
|
|
|
3257 |
> nameEx("uspop")
|
|
|
3258 |
> ### * uspop
|
|
|
3259 |
>
|
|
|
3260 |
> flush(stderr()); flush(stdout())
|
|
|
3261 |
>
|
|
|
3262 |
> ### Name: uspop
|
|
|
3263 |
> ### Title: Populations Recorded by the US Census
|
|
|
3264 |
> ### Aliases: uspop
|
|
|
3265 |
> ### Keywords: datasets
|
|
|
3266 |
>
|
|
|
3267 |
> ### ** Examples
|
|
|
3268 |
>
|
|
|
3269 |
> require(graphics)
|
|
|
3270 |
> plot(uspop, log = "y", main = "uspop data", xlab = "Year",
|
|
|
3271 |
+ ylab = "U.S. Population (millions)")
|
|
|
3272 |
>
|
|
|
3273 |
>
|
|
|
3274 |
>
|
|
|
3275 |
> cleanEx()
|
|
|
3276 |
> nameEx("volcano")
|
|
|
3277 |
> ### * volcano
|
|
|
3278 |
>
|
|
|
3279 |
> flush(stderr()); flush(stdout())
|
|
|
3280 |
>
|
|
|
3281 |
> ### Name: volcano
|
|
|
3282 |
> ### Title: Topographic Information on Auckland's Maunga Whau Volcano
|
|
|
3283 |
> ### Aliases: volcano
|
|
|
3284 |
> ### Keywords: datasets
|
|
|
3285 |
>
|
|
|
3286 |
> ### ** Examples
|
|
|
3287 |
>
|
|
|
3288 |
> require(grDevices); require(graphics)
|
|
|
3289 |
> filled.contour(volcano, color.palette = terrain.colors, asp = 1)
|
|
|
3290 |
> title(main = "volcano data: filled contour map")
|
|
|
3291 |
>
|
|
|
3292 |
>
|
|
|
3293 |
>
|
|
|
3294 |
> cleanEx()
|
|
|
3295 |
> nameEx("warpbreaks")
|
|
|
3296 |
> ### * warpbreaks
|
|
|
3297 |
>
|
|
|
3298 |
> flush(stderr()); flush(stdout())
|
|
|
3299 |
>
|
|
|
3300 |
> ### Name: warpbreaks
|
|
|
3301 |
> ### Title: The Number of Breaks in Yarn during Weaving
|
|
|
3302 |
> ### Aliases: warpbreaks
|
|
|
3303 |
> ### Keywords: datasets
|
|
|
3304 |
>
|
|
|
3305 |
> ### ** Examples
|
|
|
3306 |
>
|
|
|
3307 |
> require(stats); require(graphics)
|
|
|
3308 |
> summary(warpbreaks)
|
|
|
3309 |
breaks wool tension
|
|
|
3310 |
Min. :10.00 A:27 L:18
|
|
|
3311 |
1st Qu.:18.25 B:27 M:18
|
|
|
3312 |
Median :26.00 H:18
|
|
|
3313 |
Mean :28.15
|
|
|
3314 |
3rd Qu.:34.00
|
|
|
3315 |
Max. :70.00
|
| 61169 |
ripley |
3316 |
> opar <- par(mfrow = c(1, 2), oma = c(0, 0, 1.1, 0))
|
| 56186 |
murdoch |
3317 |
> plot(breaks ~ tension, data = warpbreaks, col = "lightgray",
|
|
|
3318 |
+ varwidth = TRUE, subset = wool == "A", main = "Wool A")
|
|
|
3319 |
> plot(breaks ~ tension, data = warpbreaks, col = "lightgray",
|
|
|
3320 |
+ varwidth = TRUE, subset = wool == "B", main = "Wool B")
|
|
|
3321 |
> mtext("warpbreaks data", side = 3, outer = TRUE)
|
|
|
3322 |
> par(opar)
|
|
|
3323 |
> summary(fm1 <- lm(breaks ~ wool*tension, data = warpbreaks))
|
|
|
3324 |
|
|
|
3325 |
Call:
|
|
|
3326 |
lm(formula = breaks ~ wool * tension, data = warpbreaks)
|
|
|
3327 |
|
|
|
3328 |
Residuals:
|
|
|
3329 |
Min 1Q Median 3Q Max
|
|
|
3330 |
-19.5556 -6.8889 -0.6667 7.1944 25.4444
|
|
|
3331 |
|
|
|
3332 |
Coefficients:
|
|
|
3333 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3334 |
(Intercept) 44.556 3.647 12.218 2.43e-16 ***
|
|
|
3335 |
woolB -16.333 5.157 -3.167 0.002677 **
|
|
|
3336 |
tensionM -20.556 5.157 -3.986 0.000228 ***
|
|
|
3337 |
tensionH -20.000 5.157 -3.878 0.000320 ***
|
|
|
3338 |
woolB:tensionM 21.111 7.294 2.895 0.005698 **
|
|
|
3339 |
woolB:tensionH 10.556 7.294 1.447 0.154327
|
|
|
3340 |
---
|
| 61435 |
ripley |
3341 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3342 |
|
|
|
3343 |
Residual standard error: 10.94 on 48 degrees of freedom
|
| 61435 |
ripley |
3344 |
Multiple R-squared: 0.3778, Adjusted R-squared: 0.3129
|
|
|
3345 |
F-statistic: 5.828 on 5 and 48 DF, p-value: 0.0002772
|
| 56186 |
murdoch |
3346 |
|
|
|
3347 |
> anova(fm1)
|
|
|
3348 |
Analysis of Variance Table
|
|
|
3349 |
|
|
|
3350 |
Response: breaks
|
|
|
3351 |
Df Sum Sq Mean Sq F value Pr(>F)
|
|
|
3352 |
wool 1 450.7 450.67 3.7653 0.0582130 .
|
|
|
3353 |
tension 2 2034.3 1017.13 8.4980 0.0006926 ***
|
|
|
3354 |
wool:tension 2 1002.8 501.39 4.1891 0.0210442 *
|
|
|
3355 |
Residuals 48 5745.1 119.69
|
|
|
3356 |
---
|
| 61435 |
ripley |
3357 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3358 |
>
|
|
|
3359 |
>
|
|
|
3360 |
>
|
|
|
3361 |
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
|
|
|
3362 |
> cleanEx()
|
|
|
3363 |
> nameEx("women")
|
|
|
3364 |
> ### * women
|
|
|
3365 |
>
|
|
|
3366 |
> flush(stderr()); flush(stdout())
|
|
|
3367 |
>
|
|
|
3368 |
> ### Name: women
|
|
|
3369 |
> ### Title: Average Heights and Weights for American Women
|
|
|
3370 |
> ### Aliases: women
|
|
|
3371 |
> ### Keywords: datasets
|
|
|
3372 |
>
|
|
|
3373 |
> ### ** Examples
|
|
|
3374 |
>
|
|
|
3375 |
> require(graphics)
|
|
|
3376 |
> plot(women, xlab = "Height (in)", ylab = "Weight (lb)",
|
|
|
3377 |
+ main = "women data: American women aged 30-39")
|
|
|
3378 |
>
|
|
|
3379 |
>
|
|
|
3380 |
>
|
|
|
3381 |
> cleanEx()
|
|
|
3382 |
> nameEx("zCO2")
|
|
|
3383 |
> ### * zCO2
|
|
|
3384 |
>
|
|
|
3385 |
> flush(stderr()); flush(stdout())
|
|
|
3386 |
>
|
|
|
3387 |
> ### Name: CO2
|
|
|
3388 |
> ### Title: Carbon Dioxide Uptake in Grass Plants
|
|
|
3389 |
> ### Aliases: CO2
|
|
|
3390 |
> ### Keywords: datasets
|
|
|
3391 |
>
|
|
|
3392 |
> ### ** Examples
|
|
|
3393 |
>
|
|
|
3394 |
> require(stats); require(graphics)
|
|
|
3395 |
> ## Don't show:
|
|
|
3396 |
> options(show.nls.convergence=FALSE)
|
|
|
3397 |
> ## End Don't show
|
|
|
3398 |
> coplot(uptake ~ conc | Plant, data = CO2, show.given = FALSE, type = "b")
|
|
|
3399 |
> ## fit the data for the first plant
|
|
|
3400 |
> fm1 <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0),
|
| 61169 |
ripley |
3401 |
+ data = CO2, subset = Plant == "Qn1")
|
| 56186 |
murdoch |
3402 |
> summary(fm1)
|
|
|
3403 |
|
|
|
3404 |
Formula: uptake ~ SSasymp(conc, Asym, lrc, c0)
|
|
|
3405 |
|
|
|
3406 |
Parameters:
|
|
|
3407 |
Estimate Std. Error t value Pr(>|t|)
|
|
|
3408 |
Asym 38.1398 0.9164 41.620 1.99e-06 ***
|
|
|
3409 |
lrc -34.2766 18.9661 -1.807 0.145
|
|
|
3410 |
c0 -4.3806 0.2042 -21.457 2.79e-05 ***
|
|
|
3411 |
---
|
| 61435 |
ripley |
3412 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
|
| 56186 |
murdoch |
3413 |
|
| 62439 |
ripley |
3414 |
Residual standard error: 1.663 on 4 degrees of freedom
|
| 56186 |
murdoch |
3415 |
|
|
|
3416 |
> ## fit each plant separately
|
|
|
3417 |
> fmlist <- list()
|
|
|
3418 |
> for (pp in levels(CO2$Plant)) {
|
|
|
3419 |
+ fmlist[[pp]] <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0),
|
|
|
3420 |
+ data = CO2, subset = Plant == pp)
|
|
|
3421 |
+ }
|
|
|
3422 |
> ## check the coefficients by plant
|
| 61157 |
ripley |
3423 |
> print(sapply(fmlist, coef), digits = 3)
|
| 56186 |
murdoch |
3424 |
Qn1 Qn2 Qn3 Qc1 Qc3 Qc2 Mn3 Mn2 Mn1 Mc2 Mc3
|
|
|
3425 |
Asym 38.14 42.87 44.23 36.43 40.68 39.82 28.48 32.13 34.08 13.56 18.54
|
|
|
3426 |
lrc -34.28 -29.66 -37.63 -9.90 -11.54 -51.53 -17.37 -29.04 -8.81 -1.98 -136.11
|
|
|
3427 |
c0 -4.38 -4.67 -4.49 -4.86 -4.95 -4.46 -4.59 -4.47 -5.06 -4.56 -3.47
|
|
|
3428 |
Mc1
|
|
|
3429 |
Asym 21.79
|
|
|
3430 |
lrc 2.45
|
|
|
3431 |
c0 -5.14
|
|
|
3432 |
>
|
|
|
3433 |
>
|
|
|
3434 |
>
|
|
|
3435 |
> ### * <FOOTER>
|
|
|
3436 |
> ###
|
| 62439 |
ripley |
3437 |
> options(digits = 7L)
|
| 61787 |
ripley |
3438 |
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
|
| 63202 |
ripley |
3439 |
Time elapsed: 1.556 0.042 1.605 0 0
|
| 56186 |
murdoch |
3440 |
> grDevices::dev.off()
|
|
|
3441 |
null device
|
|
|
3442 |
1
|
|
|
3443 |
> ###
|
|
|
3444 |
> ### Local variables: ***
|
|
|
3445 |
> ### mode: outline-minor ***
|
|
|
3446 |
> ### outline-regexp: "\\(> \\)?### [*]+" ***
|
|
|
3447 |
> ### End: ***
|
|
|
3448 |
> quit('no')
|
The R Project SVN
R