This report was automatically generated with the R package knitr (version 1.5).
library(faraway)
data(composite, package = "faraway")
composite
strength laser tape
1 25.66 40W slow
2 29.15 50W slow
3 35.73 60W slow
4 28.00 40W medium
5 35.09 50W medium
6 39.56 60W medium
7 20.65 40W fast
8 29.79 50W fast
9 35.66 60W fast
require(ggplot2)
ggplot(composite, aes(x = laser, y = strength, group = tape, linetype = tape)) +
geom_line() + theme(legend.position = "top", legend.direction = "horizontal")
ggplot(composite, aes(x = tape, y = strength, group = laser, linetype = laser)) +
geom_line() + theme(legend.position = "top", legend.direction = "horizontal")
lmod <- lm(strength ~ laser + tape + laser:tape, composite)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.66 NA NA NA
laser50W 3.49 NA NA NA
laser60W 10.07 NA NA NA
tapemedium 2.34 NA NA NA
tapefast -5.01 NA NA NA
laser50W:tapemedium 3.60 NA NA NA
laser60W:tapemedium 1.49 NA NA NA
laser50W:tapefast 5.65 NA NA NA
laser60W:tapefast 4.94 NA NA NA
n = 9, p = 9, Residual SE = NaN, R-Squared = 1
lmod <- lm(strength ~ laser + tape, composite)
coef(lmod)
(Intercept) laser50W laser60W tapemedium tapefast
23.918 6.573 12.213 4.037 -1.480
lasercoefs <- rep(c(0, 6.5733, 12.2133), 3)
tapecoefs <- rep(c(0, 4.0367, -1.48), each = 3)
tmod <- update(lmod, . ~ . + I(lasercoefs * tapecoefs))
anova(tmod)
Analysis of Variance Table
Response: strength
Df Sum Sq Mean Sq F value Pr(>F)
laser 2 224.2 112.1 36.82 0.0077
tape 2 48.9 24.5 8.03 0.0624
I(lasercoefs * tapecoefs) 1 1.4 1.4 0.45 0.5503
Residuals 3 9.1 3.0
anova(lmod)
Analysis of Variance Table
Response: strength
Df Sum Sq Mean Sq F value Pr(>F)
laser 2 224.2 112.1 42.69 0.002
tape 2 48.9 24.5 9.32 0.031
Residuals 4 10.5 2.6
composite$laser <- as.ordered(composite$laser)
composite$tape <- as.ordered(composite$tape)
lmod <- lm(strength ~ laser + tape, composite)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 31.032 0.540 57.45 5.5e-07
laser.L 8.636 0.936 9.23 0.00077
laser.Q -0.381 0.936 -0.41 0.70466
tape.L -1.047 0.936 -1.12 0.32594
tape.Q -3.900 0.936 -4.17 0.01404
n = 9, p = 5, Residual SE = 1.62, R-Squared = 0.96
round(model.matrix(lmod), 2)
(Intercept) laser.L laser.Q tape.L tape.Q
1 1 -0.71 0.41 -0.71 0.41
2 1 0.00 -0.82 -0.71 0.41
3 1 0.71 0.41 -0.71 0.41
4 1 -0.71 0.41 0.00 -0.82
5 1 0.00 -0.82 0.00 -0.82
6 1 0.71 0.41 0.00 -0.82
7 1 -0.71 0.41 0.71 0.41
8 1 0.00 -0.82 0.71 0.41
9 1 0.71 0.41 0.71 0.41
attr(,"assign")
[1] 0 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$laser
[1] "contr.poly"
attr(,"contrasts")$tape
[1] "contr.poly"
composite$Ntape <- rep(c(6.42, 13, 27), each = 3)
composite$Nlaser <- rep(c(40, 50, 60), 3)
lmodn <- lm(strength ~ Nlaser + poly(log(Ntape), 2), composite)
sumary(lmodn)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4989 3.0592 0.16 0.87684
Nlaser 0.6107 0.0604 10.11 0.00016
poly(log(Ntape), 2)1 -1.8814 1.4791 -1.27 0.25933
poly(log(Ntape), 2)2 -6.7364 1.4791 -4.55 0.00609
n = 9, p = 4, Residual SE = 1.48, R-Squared = 0.96
data(pvc, package = "faraway")
require(ggplot2)
p <- ggplot(pvc, aes(x = operator, y = psize)) + geom_point() + stat_summary(fun.y = "mean",
geom = "line", aes(group = resin))
op1means <- with(pvc[pvc$operator == 1, ], sapply(split(psize, resin), mean))
tdf <- data.frame(x = rep(0.9, 8), y = op1means, label = 1:8)
p + geom_text(data = tdf, aes(x = x, y = y, label = label))
ggplot(pvc, aes(x = resin, y = psize, shape = operator)) + geom_point() + stat_summary(fun.y = "mean",
geom = "line", aes(group = operator, linetype = operator)) + theme(legend.position = "top",
legend.direction = "horizontal")
lmod <- lm(psize ~ operator * resin, pvc)
anova(lmod)
Analysis of Variance Table
Response: psize
Df Sum Sq Mean Sq F value Pr(>F)
operator 2 20.7 10.4 7.01 0.004
resin 7 283.9 40.6 27.44 5.7e-10
operator:resin 14 14.3 1.0 0.69 0.760
Residuals 24 35.5 1.5
anova(lm(psize ~ operator + resin, pvc))
Analysis of Variance Table
Response: psize
Df Sum Sq Mean Sq F value Pr(>F)
operator 2 20.7 10.4 7.9 0.0014
resin 7 283.9 40.6 30.9 8.1e-14
Residuals 38 49.8 1.3
qqnorm(residuals(lmod), main = "")
qqline(residuals(lmod))
plot(fitted(lmod), residuals(lmod), xlab = "Fitted", ylab = "Residuals")
abline(h = 0)
plot(residuals(lmod) ~ operator, pvc, ylab = "Residuals")
pvce <- pvc[(1:24) * 2, ]
pvce$res <- sqrt(abs(residuals(lmod))[(1:24) * 2])
vmod <- lm(res ~ operator + resin, pvce)
anova(vmod)
Analysis of Variance Table
Response: res
Df Sum Sq Mean Sq F value Pr(>F)
operator 2 1.490 0.745 15.12 0.00032
resin 7 0.638 0.091 1.85 0.15545
Residuals 14 0.690 0.049
lmod <- lm(psize ~ operator + resin, pvc)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 36.240 0.523 69.34 < 2e-16
operator2 -0.263 0.405 -0.65 0.52059
operator3 -1.506 0.405 -3.72 0.00064
resin2 -1.033 0.661 -1.56 0.12630
resin3 -5.800 0.661 -8.77 1.1e-10
resin4 -6.183 0.661 -9.35 2.1e-11
resin5 -4.800 0.661 -7.26 1.1e-08
resin6 -5.450 0.661 -8.24 5.5e-10
resin7 -2.917 0.661 -4.41 8.2e-05
resin8 -0.183 0.661 -0.28 0.78302
n = 48, p = 10, Residual SE = 1.14, R-Squared = 0.86
TukeyHSD(aov(psize ~ operator + resin, data = pvc))
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = psize ~ operator + resin, data = pvc)
$operator
diff lwr upr p adj
2-1 -0.2625 -1.250 0.7247 0.7944
3-1 -1.5063 -2.493 -0.5190 0.0018
3-2 -1.2438 -2.231 -0.2565 0.0107
$resin
diff lwr upr p adj
2-1 -1.0333 -3.1523 1.0856 0.7683
3-1 -5.8000 -7.9189 -3.6811 0.0000
4-1 -6.1833 -8.3023 -4.0644 0.0000
5-1 -4.8000 -6.9189 -2.6811 0.0000
6-1 -5.4500 -7.5689 -3.3311 0.0000
7-1 -2.9167 -5.0356 -0.7977 0.0019
8-1 -0.1833 -2.3023 1.9356 1.0000
3-2 -4.7667 -6.8856 -2.6477 0.0000
4-2 -5.1500 -7.2689 -3.0311 0.0000
5-2 -3.7667 -5.8856 -1.6477 0.0000
6-2 -4.4167 -6.5356 -2.2977 0.0000
7-2 -1.8833 -4.0023 0.2356 0.1128
8-2 0.8500 -1.2689 2.9689 0.8985
4-3 -0.3833 -2.5023 1.7356 0.9989
5-3 1.0000 -1.1189 3.1189 0.7959
6-3 0.3500 -1.7689 2.4689 0.9994
7-3 2.8833 0.7644 5.0023 0.0022
8-3 5.6167 3.4977 7.7356 0.0000
5-4 1.3833 -0.7356 3.5023 0.4376
6-4 0.7333 -1.3856 2.8523 0.9508
7-4 3.2667 1.1477 5.3856 0.0004
8-4 6.0000 3.8811 8.1189 0.0000
6-5 -0.6500 -2.7689 1.4689 0.9741
7-5 1.8833 -0.2356 4.0023 0.1128
8-5 4.6167 2.4977 6.7356 0.0000
7-6 2.5333 0.4144 4.6523 0.0099
8-6 5.2667 3.1477 7.3856 0.0000
8-7 2.7333 0.6144 4.8523 0.0042
plot(breaks ~ wool, warpbreaks)
with(warpbreaks, interaction.plot(wool, tension, breaks))
require(ggplot2)
ggplot(warpbreaks, aes(x = wool, y = breaks, shape = tension)) + geom_point(position = position_jitter(width = 0.1)) +
stat_summary(fun.y = "mean", geom = "line", aes(group = tension, linetype = tension)) +
theme(legend.position = "top", legend.direction = "horizontal")
ggplot(warpbreaks, aes(x = tension, y = breaks, shape = wool)) + geom_point(position = position_jitter(width = 0.1)) +
stat_summary(fun.y = "mean", geom = "line", aes(group = wool, linetype = wool)) +
theme(legend.position = "top", legend.direction = "horizontal")
lmod <- lm(breaks ~ wool * tension, warpbreaks)
plot(residuals(lmod) ~ fitted(lmod), xlab = "Fitted", ylab = "Residuals")
abline(h = 0)
lmod <- lm(sqrt(breaks) ~ wool * tension, warpbreaks)
plot(residuals(lmod) ~ fitted(lmod), xlab = "Fitted", ylab = "Residuals")
abline(h = 0)
anova(lmod)
Analysis of Variance Table
Response: sqrt(breaks)
Df Sum Sq Mean Sq F value Pr(>F)
wool 1 2.9 2.90 3.02 0.08854
tension 2 15.9 7.95 8.28 0.00082
wool:tension 2 7.2 3.60 3.75 0.03067
Residuals 48 46.1 0.96
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.548 0.327 20.05 < 2e-16
woolB -1.309 0.462 -2.83 0.00669
tensionM -1.722 0.462 -3.73 0.00051
tensionH -1.691 0.462 -3.66 0.00062
woolB:tensionM 1.782 0.653 2.73 0.00887
woolB:tensionH 0.755 0.653 1.16 0.25335
n = 54, p = 6, Residual SE = 0.98, R-Squared = 0.36
lmod <- lm(sqrt(breaks) ~ wool:tension - 1, warpbreaks)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
woolA:tensionL 6.548 0.327 20.1 <2e-16
woolB:tensionL 5.238 0.327 16.0 <2e-16
woolA:tensionM 4.826 0.327 14.8 <2e-16
woolB:tensionM 5.299 0.327 16.2 <2e-16
woolA:tensionH 4.856 0.327 14.9 <2e-16
woolB:tensionH 4.302 0.327 13.2 <2e-16
n = 54, p = 6, Residual SE = 0.98, R-Squared = 0.97
TukeyHSD(aov(breaks ~ wool:tension, warpbreaks))
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = breaks ~ wool:tension, data = warpbreaks)
$`wool:tension`
diff lwr upr p adj
B:L-A:L -16.3333 -31.64 -1.0270 0.0302
A:M-A:L -20.5556 -35.86 -5.2492 0.0030
B:M-A:L -15.7778 -31.08 -0.4715 0.0398
A:H-A:L -20.0000 -35.31 -4.6937 0.0041
B:H-A:L -25.7778 -41.08 -10.4715 0.0001
A:M-B:L -4.2222 -19.53 11.0841 0.9627
B:M-B:L 0.5556 -14.75 15.8619 1.0000
A:H-B:L -3.6667 -18.97 11.6397 0.9797
B:H-B:L -9.4444 -24.75 5.8619 0.4561
B:M-A:M 4.7778 -10.53 20.0841 0.9377
A:H-A:M 0.5556 -14.75 15.8619 1.0000
B:H-A:M -5.2222 -20.53 10.0841 0.9115
A:H-B:M -4.2222 -19.53 11.0841 0.9627
B:H-B:M -10.0000 -25.31 5.3063 0.3919
B:H-A:H -5.7778 -21.08 9.5285 0.8706
data(speedo, package = "faraway")
speedo
h d l b j f n a i e m c k g o y
1 - - + - + + - - + + - + - - + 0.4850
2 + - - - - + + - - + + + + - - 0.5750
3 - + - - + - + - + - + + - + - 0.0875
4 + + + - - - - - - - - + + + + 0.1750
5 - - + + - - + - + + - - + + - 0.1950
6 + - - + + - - - - + + - - + + 0.1450
7 - + - + - + - - + - + - + - + 0.2250
8 + + + + + + + - - - - - - - - 0.1750
9 - - + - + + - + - - + - + + - 0.1250
10 + - - - - + + + + - - - - + + 0.1200
11 - + - - + - + + - + - - + - + 0.4550
12 + + + - - - - + + + + - - - - 0.5350
13 - - + + - - + + - - + + - - + 0.1700
14 + - - + + - - + + - - + + - - 0.2750
15 - + - + - + - + - + - + - + - 0.3425
16 + + + + + + + + + + + + + + + 0.5825
lmod <- lm(y ~ ., speedo)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.582500 NA NA NA
h- -0.062188 NA NA NA
d- -0.060938 NA NA NA
l- -0.027188 NA NA NA
b- 0.055937 NA NA NA
j- 0.000937 NA NA NA
f- -0.074062 NA NA NA
n- -0.006563 NA NA NA
a- -0.067813 NA NA NA
i- -0.042813 NA NA NA
e- -0.245312 NA NA NA
m- -0.027812 NA NA NA
c- -0.089687 NA NA NA
k- -0.068437 NA NA NA
g- 0.140312 NA NA NA
o- -0.005938 NA NA NA
n = 16, p = 16, Residual SE = NaN, R-Squared = 1
model.matrix(lmod)
(Intercept) h- d- l- b- j- f- n- a- i- e- m- c- k- g- o-
1 1 1 1 0 1 0 0 1 1 0 0 1 0 1 1 0
2 1 0 1 1 1 1 0 0 1 1 0 0 0 0 1 1
3 1 1 0 1 1 0 1 0 1 0 1 0 0 1 0 1
4 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0
5 1 1 1 0 0 1 1 0 1 0 0 1 1 0 0 1
6 1 0 1 1 0 0 1 1 1 1 0 0 1 1 0 0
7 1 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0
8 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
9 1 1 1 0 1 0 0 1 0 1 1 0 1 0 0 1
10 1 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0
11 1 1 0 1 1 0 1 0 0 1 0 1 1 0 1 0
12 1 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1
13 1 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0
14 1 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1
15 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
16 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
attr(,"assign")
[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
attr(,"contrasts")
attr(,"contrasts")$h
[1] "contr.treatment"
attr(,"contrasts")$d
[1] "contr.treatment"
attr(,"contrasts")$l
[1] "contr.treatment"
attr(,"contrasts")$b
[1] "contr.treatment"
attr(,"contrasts")$j
[1] "contr.treatment"
attr(,"contrasts")$f
[1] "contr.treatment"
attr(,"contrasts")$n
[1] "contr.treatment"
attr(,"contrasts")$a
[1] "contr.treatment"
attr(,"contrasts")$i
[1] "contr.treatment"
attr(,"contrasts")$e
[1] "contr.treatment"
attr(,"contrasts")$m
[1] "contr.treatment"
attr(,"contrasts")$c
[1] "contr.treatment"
attr(,"contrasts")$k
[1] "contr.treatment"
attr(,"contrasts")$g
[1] "contr.treatment"
attr(,"contrasts")$o
[1] "contr.treatment"
qqnorm(coef(lmod)[-1], pch = names(coef(lmod)[-1]))
halfnorm(coef(lmod)[-1], labs = names(coef(lmod)[-1]))
The R session information (including the OS info, R version and all packages used):
sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] graphics grDevices utils datasets methods stats base
other attached packages:
[1] faraway_1.0.6 knitr_1.5 ggplot2_0.9.3.1
loaded via a namespace (and not attached):
[1] colorspace_1.2-4 dichromat_2.0-0 digest_0.6.4
[4] evaluate_0.5.3 formatR_0.10 grid_3.1.0
[7] gtable_0.1.2 labeling_0.2 MASS_7.3-31
[10] munsell_0.4.2 plyr_1.8.1 proto_0.3-10
[13] RColorBrewer_1.0-5 Rcpp_0.11.1 reshape2_1.2.2
[16] scales_0.2.3 stringr_0.6.2 tools_3.1.0
Sys.time()
[1] "2014-06-16 14:02:34 BST"