#I have written the interpretation in my copy,that's why I didn't add any comments on my R ouput
#Problem 13
strength <- c(25.66, 29.15, 35.73,
28.00, 35.09, 39.56,
20.65, 29.79, 35.66)
speed <- c(6.42, 6.42, 6.42,
13.00, 13.00, 13.00,
27.00, 27.00, 27.00)
power <- c(40, 50, 60,
40, 50, 60,
40, 50, 60)
dat <- data.frame(speed, power, strength)
dat$x1 <- (dat$power - 50)/10
dat$x2 <- ifelse(dat$speed == 6.42, -1,
ifelse(dat$speed == 13.00, 0, 1))
#Full model
mod_full <- lm(strength ~ x1 + x2 + I(x1^2) + I(x2^2) + x1:x2, data = dat)
summary(mod_full)
##
## Call:
## lm(formula = strength ~ x1 + x2 + I(x1^2) + I(x2^2) + x1:x2,
## data = dat)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## 0.50722 -1.34111 0.83389 0.04556 0.56222 -0.60778 -0.55278 0.77889
## 9
## -0.22611
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.5278 0.9029 38.241 3.93e-05 ***
## x1 6.1067 0.4945 12.348 0.00114 **
## x2 -0.7400 0.4945 -1.496 0.23146
## I(x1^2) -0.4667 0.8566 -0.545 0.62377
## I(x2^2) -4.7767 0.8566 -5.577 0.01138 *
## x1:x2 1.2350 0.6057 2.039 0.13417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.211 on 3 degrees of freedom
## Multiple R-squared: 0.9845, Adjusted R-squared: 0.9586
## F-statistic: 38.05 on 5 and 3 DF, p-value: 0.006475
anova(mod_full)
## Analysis of Variance Table
##
## Response: strength
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 223.748 223.748 152.4813 0.001144 **
## x2 1 3.286 3.286 2.2391 0.231460
## I(x1^2) 1 0.436 0.436 0.2968 0.623770
## I(x2^2) 1 45.633 45.633 31.0983 0.011383 *
## x1:x2 1 6.101 6.101 4.1577 0.134171
## Residuals 3 4.402 1.467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Reduced model
mod_reduced <- lm(strength ~ x1 + x2 + I(x1^2) + I(x2^2), data = dat)
summary(mod_reduced)
##
## Call:
## lm(formula = strength ~ x1 + x2 + I(x1^2) + I(x2^2), data = dat)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## 1.74222 -1.34111 -0.40111 0.04556 0.56222 -0.60778 -1.78778 0.77889
## 9
## 1.00889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.5278 1.2078 28.588 8.91e-06 ***
## x1 6.1067 0.6615 9.231 0.000765 ***
## x2 -0.7400 0.6615 -1.119 0.325944
## I(x1^2) -0.4667 1.1458 -0.407 0.704655
## I(x2^2) -4.7767 1.1458 -4.169 0.014045 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.62 on 4 degrees of freedom
## Multiple R-squared: 0.963, Adjusted R-squared: 0.9259
## F-statistic: 26 on 4 and 4 DF, p-value: 0.004013
anova(mod_reduced)
## Analysis of Variance Table
##
## Response: strength
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 223.748 223.748 85.2127 0.0007654 ***
## x2 1 3.286 3.286 1.2513 0.3259444
## I(x1^2) 1 0.436 0.436 0.1659 0.7046552
## I(x2^2) 1 45.633 45.633 17.3790 0.0140447 *
## Residuals 4 10.503 2.626
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comparison of models
anova(mod_reduced, mod_full)
## Analysis of Variance Table
##
## Model 1: strength ~ x1 + x2 + I(x1^2) + I(x2^2)
## Model 2: strength ~ x1 + x2 + I(x1^2) + I(x2^2) + x1:x2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 4 10.5030
## 2 3 4.4021 1 6.1009 4.1577 0.1342
#model Diagnostics
par(mfrow = c(2,2))
plot(mod_reduced)
## hat values (leverages) are all = 0.5555556
## and there are no factor predictors; no plot no. 5
residuals(mod_reduced)
## 1 2 3 4 5 6
## 1.74222222 -1.34111111 -0.40111111 0.04555556 0.56222222 -0.60777778
## 7 8 9
## -1.78777778 0.77888889 1.00888889
fitted(mod_reduced)
## 1 2 3 4 5 6 7 8
## 23.91778 30.49111 36.13111 27.95444 34.52778 40.16778 22.43778 29.01111
## 9
## 34.65111
rstandard(mod_reduced)
## 1 2 3 4 5 6
## 1.61275056 -1.24144766 -0.37130290 0.04217014 0.52044119 -0.56261133
## 7 8 9
## -1.65492070 0.72100647 0.93391423
cooks.distance(mod_reduced)
## 1 2 3 4 5 6
## 0.6502410941 0.3852980723 0.0344664616 0.0004445801 0.0677147578 0.0791328758
## 7 8 9
## 0.6846906294 0.1299625823 0.2180489467
shapiro.test(residuals(mod_reduced))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod_reduced)
## W = 0.9789, p-value = 0.9584
#problem 28
plate <- factor(c(1,1,1, 2,2,2, 3,3,3, 4,4,4))
shape <- factor(c("A","C","D", "A","B","D", "A","B","C", "B","C","D"))
noise <- c(1.11,0.95,0.82,
1.70,1.22,0.97,
1.60,1.11,1.52,
1.22,1.54,1.18)
dat <- data.frame(plate, shape, noise)
dat
## plate shape noise
## 1 1 A 1.11
## 2 1 C 0.95
## 3 1 D 0.82
## 4 2 A 1.70
## 5 2 B 1.22
## 6 2 D 0.97
## 7 3 A 1.60
## 8 3 B 1.11
## 9 3 C 1.52
## 10 4 B 1.22
## 11 4 C 1.54
## 12 4 D 1.18
# Treatment means
aggregate(noise ~ shape, data = dat, mean)
## shape noise
## 1 A 1.470000
## 2 B 1.183333
## 3 C 1.336667
## 4 D 0.990000
#Block means
aggregate(noise ~ plate, data = dat, mean)
## plate noise
## 1 1 0.960000
## 2 2 1.296667
## 3 3 1.410000
## 4 4 1.313333
#additive incomplete block model
mod <- lm(noise ~ plate + shape, data = dat)
# Model summary
summary(mod)
##
## Call:
## lm(formula = noise ~ plate + shape, data = dat)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.07000 -0.07375 0.14375 0.08375 0.05875 -0.14250 -0.01375 -0.04875
## 9 10 11 12
## 0.06250 -0.01000 0.01125 -0.00125
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.18000 0.08938 13.202 4.45e-05 ***
## plate2 0.43625 0.10135 4.305 0.00768 **
## plate3 0.43375 0.10135 4.280 0.00786 **
## plate4 0.50500 0.10135 4.983 0.00417 **
## shapeB -0.45500 0.10135 -4.490 0.00646 **
## shapeC -0.15625 0.10135 -1.542 0.18378
## shapeD -0.50375 0.10135 -4.971 0.00421 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.117 on 5 degrees of freedom
## Multiple R-squared: 0.9223, Adjusted R-squared: 0.829
## F-statistic: 9.887 on 6 and 5 DF, p-value: 0.01184
#ANOVA table
anova(mod)
## Analysis of Variance Table
##
## Response: noise
## Df Sum Sq Mean Sq F value Pr(>F)
## plate 3 0.34737 0.115789 8.4548 0.02105 *
## shape 3 0.46506 0.155019 11.3194 0.01146 *
## Residuals 5 0.06847 0.013695
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Diagnostics
par(mfrow = c(2,2))
plot(mod)
# Residual
residuals(mod)
## 1 2 3 4 5 6 7 8
## -0.07000 -0.07375 0.14375 0.08375 0.05875 -0.14250 -0.01375 -0.04875
## 9 10 11 12
## 0.06250 -0.01000 0.01125 -0.00125
fitted(mod)
## 1 2 3 4 5 6 7 8 9 10
## 1.18000 1.02375 0.67625 1.61625 1.16125 1.11250 1.61375 1.15875 1.45750 1.23000
## 11 12
## 1.52875 1.18125
shapiro.test(residuals(mod))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod)
## W = 0.98581, p-value = 0.9975
#problem 38
WholePlot <- factor(c(1,1,1,1,
2,2,2,2,
3,3,3,3,
4,4,4,4,
5,5,5,5,
6,6,6,6))
T <- factor(c(1,1,1,1,
2,2,2,2,
3,3,3,3,
1,1,1,1,
2,2,2,2,
3,3,3,3))
C <- factor(c(2,3,1,4,
1,3,4,2,
3,1,2,4,
4,3,2,1,
4,1,3,2,
1,4,2,3))
Rep <- factor(c(1,1,1,1,
1,1,1,1,
1,1,1,1,
2,2,2,2,
2,2,2,2,
2,2,2,2))
Resistance <- c(73,83,67,89,
65,87,86,91,
147,155,127,212,
153,90,100,108,
150,140,121,142,
33,54,8,46)
dat <- data.frame(WholePlot, T, C, Rep, Resistance)
dat
## WholePlot T C Rep Resistance
## 1 1 1 2 1 73
## 2 1 1 3 1 83
## 3 1 1 1 1 67
## 4 1 1 4 1 89
## 5 2 2 1 1 65
## 6 2 2 3 1 87
## 7 2 2 4 1 86
## 8 2 2 2 1 91
## 9 3 3 3 1 147
## 10 3 3 1 1 155
## 11 3 3 2 1 127
## 12 3 3 4 1 212
## 13 4 1 4 2 153
## 14 4 1 3 2 90
## 15 4 1 2 2 100
## 16 4 1 1 2 108
## 17 5 2 4 2 150
## 18 5 2 1 2 140
## 19 5 2 3 2 121
## 20 5 2 2 2 142
## 21 6 3 1 2 33
## 22 6 3 4 2 54
## 23 6 3 2 2 8
## 24 6 3 3 2 46
fit <- aov(Resistance ~ Rep + T*C + Error(Rep/T), data = dat)
summary(fit)
##
## Error: Rep
## Df Sum Sq Mean Sq
## Rep 1 782 782
##
## Error: Rep:T
## Df Sum Sq Mean Sq F value Pr(>F)
## T 2 1022 511 0.026 0.975
## Residuals 2 39155 19578
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## C 3 4289 1429.7 5.891 0.0166 *
## T:C 6 2206 367.7 1.515 0.2759
## Residuals 9 2184 242.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aggregate(Resistance ~ T, data = dat, mean)
## T Resistance
## 1 1 95.375
## 2 2 110.250
## 3 3 97.750
aggregate(Resistance ~ C, data = dat, mean)
## C Resistance
## 1 1 94.66667
## 2 2 90.16667
## 3 3 95.66667
## 4 4 124.00000
aggregate(Resistance ~ T + C, data = dat, mean)
## T C Resistance
## 1 1 1 87.5
## 2 2 1 102.5
## 3 3 1 94.0
## 4 1 2 86.5
## 5 2 2 116.5
## 6 3 2 67.5
## 7 1 3 86.5
## 8 2 3 104.0
## 9 3 3 96.5
## 10 1 4 121.0
## 11 2 4 118.0
## 12 3 4 133.0