The R Project SVN – R – Blame – /trunk/tests/Examples/datasets-Ex.Rout.save

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

The R Project SVN R

(root)/trunk/tests/Examples/datasets-Ex.Rout.save – Rev 77458