df_A <- 1
ms_A <- 0.0002
df_B <- 3
ss_B <- 180.378
df_AB <- 3
ss_AB <- 8.479
df_E <- 8
ss_E <- 158.797
ms_B <- ss_B / df_B
ms_AB <- ss_AB / df_AB
ms_E <- ss_E / df_E
F_A <- ms_A / ms_E
F_B <- ms_B / ms_E
F_AB <- ms_AB / ms_E
p_A <- pf(F_A, df_A, df_E, lower.tail = FALSE)
p_B <- pf(F_B, df_B, df_E, lower.tail = FALSE)
p_AB <- pf(F_AB, df_AB, df_E, lower.tail = FALSE)
anova_table <- data.frame(
Source = c("A", "B", "Interaction", "Error"),
DF = c(df_A, df_B, df_AB, df_E),
SS = c(ms_A * df_A, ss_B, ss_AB, ss_E),
MS = c(ms_A, ms_B, ms_AB, ms_E),
F = c(F_A, F_B, F_AB, NA),
P = c(p_A, p_B, p_AB, NA)
)
print(anova_table)
## Source DF SS MS F P
## 1 A 1 0.0002 0.000200 1.007576e-05 0.99754506
## 2 B 3 180.3780 60.126000 3.029075e+00 0.09334682
## 3 Interaction 3 8.4790 2.826333 1.423872e-01 0.93169691
## 4 Error 8 158.7970 19.849625 NA NA
\[DF_B=b-1 \]\[ b=3+1=4\]So factor B has 4 levels.
\[a=1+1=2 \]\[ ab=2\times4=8 \]\[ N=15+1=16 \]\[ n=\frac{16}{8}=2\]The experiment was replicated twice.
At the 0.05 significance level, neither factor A nor factor B has a statistically significant effect on the response. The interaction between A and B is also not significant. Therefore, there is insufficient evidence to conclude that either factor or their interaction affects the response variable.
Speed <- rep(c(125, 200), each = 8)
Feed <- rep(rep(c(0.015, 0.030, 0.045, 0.060),each = 2), 2)
Thrust <- c(
2.70, 2.78,
2.45, 2.49,
2.60, 2.72,
2.75, 2.86,
2.83, 2.86,
2.85, 2.80,
2.86, 2.87,
2.94, 2.88)
dat <- data.frame(Speed, Feed, Thrust)
dat$Speed <- as.factor(dat$Speed)
dat$Feed <- as.factor(dat$Feed)
dat
## Speed Feed Thrust
## 1 125 0.015 2.70
## 2 125 0.015 2.78
## 3 125 0.03 2.45
## 4 125 0.03 2.49
## 5 125 0.045 2.60
## 6 125 0.045 2.72
## 7 125 0.06 2.75
## 8 125 0.06 2.86
## 9 200 0.015 2.83
## 10 200 0.015 2.86
## 11 200 0.03 2.85
## 12 200 0.03 2.80
## 13 200 0.045 2.86
## 14 200 0.045 2.87
## 15 200 0.06 2.94
## 16 200 0.06 2.88
The full factorial model is: \[Y_{ijk}=\mu+\alpha_i+\beta_j+ (\alpha\beta)_{ij}+\epsilon_{ijk}\]
model <- aov(Thrust ~ Speed * Feed, data = dat)
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## Speed 1 0.14823 0.14823 57.010 6.61e-05 ***
## Feed 3 0.09250 0.03083 11.859 0.00258 **
## Speed:Feed 3 0.04187 0.01396 5.369 0.02557 *
## Residuals 8 0.02080 0.00260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, both drill speed and feed rate have significant effects on thrust force. More importantly, the interaction between drill speed and feed rate is also significant (p = 0.025567). Therefore, the effect of feed rate on thrust force depends on the drill speed.
aggregate(Thrust ~ Speed + Feed, data = dat, mean)
## Speed Feed Thrust
## 1 125 0.015 2.740
## 2 200 0.015 2.845
## 3 125 0.03 2.470
## 4 200 0.03 2.825
## 5 125 0.045 2.660
## 6 200 0.045 2.865
## 7 125 0.06 2.805
## 8 200 0.06 2.910
interaction.plot(Feed, Speed, Thrust)
The interaction plot shows that the two lines are not parallel, indicating an interaction between drill speed and feed rate. The effect of feed rate on thrust force depends on the drill speed.
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))
The residual plots can be used to check constant variance and model adequacy.
At the 0.05 significance level, drill speed, feed rate, and their interaction significantly affect thrust force. The lowest observed mean thrust force is 2.470, occurring at a drill speed of 125 and a feed rate of 0.030. Therefore, this combination provides the lowest average thrust force among the tested conditions.
Feed <- rep(c(0.20, 0.25, 0.30), each = 12)
Depth <- rep(rep(c(0.15, 0.18, 0.20, 0.25),each = 3), 3)
Finish <- c(74, 64, 60, 79, 68, 73, 82, 88, 92, 99, 104, 96,
92, 86, 88, 98, 104, 88, 99, 108, 95, 104, 110, 99,
99, 98, 102, 104, 99, 95, 108, 110, 99, 114, 111, 107)
dat <- data.frame(Feed, Depth, Finish)
dat$Feed <- as.factor(dat$Feed)
dat$Depth <- as.factor(dat$Depth)
dat
## Feed Depth Finish
## 1 0.2 0.15 74
## 2 0.2 0.15 64
## 3 0.2 0.15 60
## 4 0.2 0.18 79
## 5 0.2 0.18 68
## 6 0.2 0.18 73
## 7 0.2 0.2 82
## 8 0.2 0.2 88
## 9 0.2 0.2 92
## 10 0.2 0.25 99
## 11 0.2 0.25 104
## 12 0.2 0.25 96
## 13 0.25 0.15 92
## 14 0.25 0.15 86
## 15 0.25 0.15 88
## 16 0.25 0.18 98
## 17 0.25 0.18 104
## 18 0.25 0.18 88
## 19 0.25 0.2 99
## 20 0.25 0.2 108
## 21 0.25 0.2 95
## 22 0.25 0.25 104
## 23 0.25 0.25 110
## 24 0.25 0.25 99
## 25 0.3 0.15 99
## 26 0.3 0.15 98
## 27 0.3 0.15 102
## 28 0.3 0.18 104
## 29 0.3 0.18 99
## 30 0.3 0.18 95
## 31 0.3 0.2 108
## 32 0.3 0.2 110
## 33 0.3 0.2 99
## 34 0.3 0.25 114
## 35 0.3 0.25 111
## 36 0.3 0.25 107
The full factorial model is: \[Y_{ijk}=\mu+\alpha_i+\beta_j+ (\alpha\beta)_{ij}+\epsilon_{ijk}\]
model <- aov(Finish ~ Feed * Depth, data = dat)
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## Feed 2 3160.5 1580.2 55.018 1.09e-09 ***
## Depth 3 2125.1 708.4 24.663 1.65e-07 ***
## Feed:Depth 6 557.1 92.8 3.232 0.018 *
## Residuals 24 689.3 28.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, feed rate, depth of cut, and their interaction are all statistically significant. Therefore, both factors affect surface finish, and their effects depend on each other.
interaction.plot(Depth, Feed, Finish)
The interaction plot indicates that the effect of depth of cut on surface finish differs across feed rates. This agrees with the significant interaction found in the ANOVA.
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))
The residual plots should also be checked for constant variance and any systematic patterns. Overall, there is no strong evidence of a normality problem.
aggregate(Finish ~ Feed, data = dat, mean)
## Feed Finish
## 1 0.2 81.58333
## 2 0.25 97.58333
## 3 0.3 103.83333
The estimated mean surface finish increases as the feed rate increases.
summary(model)[[1]][["Pr(>F)"]]
## [1] 1.086046e-09 1.652000e-07 1.797302e-02 NA
Feed Rate: \(p=1.086\times10^{-9}\); Depth of Cut: \(p=1.652\times10^{-7}\); Interaction: \(p=0.01797\)
All three p-values are less than 0.05. Therefore, all effects are statistically significant.
At the 0.05 significance level, feed rate, depth of cut, and their interaction significantly affect surface finish. The interaction indicates that the effect of depth of cut depends on the feed rate. The estimated mean surface finish is highest at a feed rate of 0.30.
Feed <- rep(c(0.20, 0.25, 0.30), each = 12)
Depth <- rep(rep(c(0.15, 0.18, 0.20, 0.25),each = 3), 3)
Finish <- c(74, 64, 60, 79, 68, 73, 82, 88, 92, 99, 104, 96,
92, 86, 88, 98, 104, 88, 99, 108, 95, 104, 110, 99,
99, 98, 102, 104, 99, 95, 108, 110, 99, 114, 111, 107)
Block <- rep(1:3, times = 12)
dat <- data.frame(Feed, Depth, Block, Finish)
dat$Feed <- as.factor(dat$Feed)
dat$Depth <- as.factor(dat$Depth)
dat$Block <- as.factor(dat$Block)
dat
## Feed Depth Block Finish
## 1 0.2 0.15 1 74
## 2 0.2 0.15 2 64
## 3 0.2 0.15 3 60
## 4 0.2 0.18 1 79
## 5 0.2 0.18 2 68
## 6 0.2 0.18 3 73
## 7 0.2 0.2 1 82
## 8 0.2 0.2 2 88
## 9 0.2 0.2 3 92
## 10 0.2 0.25 1 99
## 11 0.2 0.25 2 104
## 12 0.2 0.25 3 96
## 13 0.25 0.15 1 92
## 14 0.25 0.15 2 86
## 15 0.25 0.15 3 88
## 16 0.25 0.18 1 98
## 17 0.25 0.18 2 104
## 18 0.25 0.18 3 88
## 19 0.25 0.2 1 99
## 20 0.25 0.2 2 108
## 21 0.25 0.2 3 95
## 22 0.25 0.25 1 104
## 23 0.25 0.25 2 110
## 24 0.25 0.25 3 99
## 25 0.3 0.15 1 99
## 26 0.3 0.15 2 98
## 27 0.3 0.15 3 102
## 28 0.3 0.18 1 104
## 29 0.3 0.18 2 99
## 30 0.3 0.18 3 95
## 31 0.3 0.2 1 108
## 32 0.3 0.2 2 110
## 33 0.3 0.2 3 99
## 34 0.3 0.25 1 114
## 35 0.3 0.25 2 111
## 36 0.3 0.25 3 107
model <- aov(Finish ~ Block + Feed * Depth, data = dat)
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 2 180.7 90.3 3.907 0.03532 *
## Feed 2 3160.5 1580.2 68.346 3.64e-10 ***
## Depth 3 2125.1 708.4 30.637 4.89e-08 ***
## Feed:Depth 6 557.1 92.8 4.015 0.00726 **
## Residuals 22 508.7 23.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, blocks, feed rate, depth of cut, and their interaction are all statistically significant.
anova_table <- anova(model)
MS_block <- anova_table["Block", "Mean Sq"]
MS_error <- anova_table["Residuals", "Mean Sq"]
var_block <- (MS_block - MS_error) / (3 * 4)
var_block
## [1] 5.60101
The estimated block variance component is 5.6010, indicating variation between blocks. Since the block effect is statistically significant (p = 0.03532), blocking was useful in this experiment.
Time <- rep(rep(c(12, 18), each = 12), 1)
Medium <- rep(rep(c(1, 2), each = 6), 2)
Growth <- c( 21, 22, 23, 28, 20, 26,
25, 26, 24, 25, 29, 27,
37, 39, 38, 38, 35, 36,
31, 34, 29, 33, 30, 35)
dat <- data.frame(Time, Medium, Growth)
dat$Time <- as.factor(dat$Time)
dat$Medium <- as.factor(dat$Medium)
dat
## Time Medium Growth
## 1 12 1 21
## 2 12 1 22
## 3 12 1 23
## 4 12 1 28
## 5 12 1 20
## 6 12 1 26
## 7 12 2 25
## 8 12 2 26
## 9 12 2 24
## 10 12 2 25
## 11 12 2 29
## 12 12 2 27
## 13 18 1 37
## 14 18 1 39
## 15 18 1 38
## 16 18 1 38
## 17 18 1 35
## 18 18 1 36
## 19 18 2 31
## 20 18 2 34
## 21 18 2 29
## 22 18 2 33
## 23 18 2 30
## 24 18 2 35
The full factorial model is: \[Y_{ijk}=\mu+\alpha_i+\beta_j+ (\alpha\beta)_{ij}+\epsilon_{ijk}\]
model <- aov(Growth ~ Time * Medium, data = dat)
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## Time 1 590.0 590.0 115.506 9.29e-10 ***
## Medium 1 9.4 9.4 1.835 0.190617
## Time:Medium 1 92.0 92.0 18.018 0.000397 ***
## Residuals 20 102.2 5.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, time and the interaction between time and culture medium are significant. Culture medium alone is not significant.
interaction.plot(Time, Medium, Growth)
The interaction plot shows that the effect of culture medium depends on incubation time. Medium 2 produces higher growth at 12 hours, while Medium 1 produces higher growth at 18 hours.
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))
The residual plots can also be used to check constant variance and model adequacy. Overall, there is no strong evidence of model inadequacy.
A <- rep(c(-1, 1, -1, 1), each = 4)
B <- rep(c(-1, -1, 1, 1), each = 4)
y <- c(14.037, 16.165, 13.972, 13.907,
13.880, 13.860, 14.032, 13.914,
14.821, 14.757, 14.843, 14.878,
14.888, 14.921, 14.415, 14.932)
dat <- data.frame(A, B, y)
dat
## A B y
## 1 -1 -1 14.037
## 2 -1 -1 16.165
## 3 -1 -1 13.972
## 4 -1 -1 13.907
## 5 1 -1 13.880
## 6 1 -1 13.860
## 7 1 -1 14.032
## 8 1 -1 13.914
## 9 -1 1 14.821
## 10 -1 1 14.757
## 11 -1 1 14.843
## 12 -1 1 14.878
## 13 1 1 14.888
## 14 1 1 14.921
## 15 1 1 14.415
## 16 1 1 14.932
model <- lm(y ~ A * B, data = dat)
coef(model)
## (Intercept) A B A:B
## 14.513875 -0.158625 0.293000 0.140750
effects <- 2 * coef(model)[-1]
effects
## A B A:B
## -0.31725 0.58600 0.28150
Factor B has the largest estimated effect (0.586). Factor A has a negative effect (-0.31725), while the interaction effect is positive (0.28150).
anova(model)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 0.4026 0.40259 1.2619 0.28327
## B 1 1.3736 1.37358 4.3054 0.06016 .
## A:B 1 0.3170 0.31697 0.9935 0.33856
## Residuals 12 3.8285 0.31904
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, none of the factors or interactions are statistically significant. Factor B has the smallest p-value (0.0602), but it is still greater than 0.05.
coef(model)
## (Intercept) A B A:B
## 14.513875 -0.158625 0.293000 0.140750
The estimated regression equation using coded variables is:
\[\hat{y}=14.513875-0.158625A +0.293B+0.14075AB\]
par(mfrow = c(2, 2))
plot(model)
## hat values (leverages) are all = 0.25
## and there are no factor predictors; no plot no. 5
par(mfrow = c(1, 1))
cbind(dat, Residual = residuals(model),
Standardized_Residual = rstandard(model))
## A B y Residual Standardized_Residual
## 1 -1 -1 14.037 -0.48325 -0.987912850
## 2 -1 -1 16.165 1.64475 3.362379016
## 3 -1 -1 13.972 -0.54825 -1.120793005
## 4 -1 -1 13.907 -0.61325 -1.253673161
## 5 1 -1 13.880 -0.04150 -0.084838869
## 6 1 -1 13.860 -0.06150 -0.125725070
## 7 1 -1 14.032 0.11050 0.225896265
## 8 1 -1 13.914 -0.00750 -0.015332326
## 9 -1 1 14.821 -0.00375 -0.007666163
## 10 -1 1 14.757 -0.06775 -0.138502008
## 11 -1 1 14.843 0.01825 0.037308659
## 12 -1 1 14.878 0.05325 0.108859512
## 13 1 1 14.888 0.09900 0.202386699
## 14 1 1 14.921 0.13200 0.269848932
## 15 1 1 14.415 -0.37400 -0.764571973
## 16 1 1 14.932 0.14300 0.292336343
The residual plots indicate a potential outlier. Observation 2 has a standardized residual of approximately 3.36, which is greater than 3. This observation corresponds to a thickness of 16.165 μm at A = -1 and B = -1. It should be investigated further.
The original measurement and experimental records should be checked to determine whether the unusual observation resulted from a measurement or recording error. If an assignable cause is identified, the observation may be corrected or excluded with justification. Otherwise, it should be retained, and the analysis should be reported with and without the observation for comparison.
A: Length of Putt (10 ft = −1, 30 ft = +1)
B: Type of Putter (Mallet = −1, Cavity back = +1)
C: Break of Putt (Straight = −1, Breaking = +1)
D: Slope of Putt (Level = −1, Downhill = +1)
A <- rep(c(-1, 1, -1, 1,
-1, 1, -1, 1,
-1, 1, -1, 1,
-1, 1, -1, 1), each = 7)
B <- rep(c(-1, -1, 1, 1,
-1, -1, 1, 1,
-1, -1, 1, 1,
-1, -1, 1, 1), each = 7)
C <- rep(c(-1, -1, -1, -1,
1, 1, 1, 1,
-1, -1, -1, -1,
1, 1, 1, 1), each = 7)
D <- rep(c(rep(-1, 8), rep(1, 8)), each = 7)
y <- c(10.0, 18.0, 14.0, 12.5, 19.0, 16.0, 18.5,
0.0, 16.5, 4.5, 17.5, 20.5, 17.5, 33.0,
4.0, 6.0, 1.0, 14.5, 12.0, 14.0, 5.0,
0.0, 10.0, 34.0, 11.0, 25.5, 21.5, 0.0,
0.0, 0.0, 18.5, 19.5, 16.0, 15.0, 11.0,
5.0, 20.5, 18.0, 20.0, 29.5, 19.0, 10.0,
6.5, 18.5, 7.5, 6.0, 0.0, 10.0, 0.0,
16.5, 4.5, 0.0, 23.5, 8.0, 8.0, 8.0,
4.5, 18.0, 14.5, 10.0, 8.0, 17.5, 6.0,
19.5, 18.0, 16.0, 5.5, 10.0, 7.0, 36.0,
15.0, 16.0, 8.5, 0.0, 0.5, 9.0, 3.0,
41.5, 39.0, 6.5, 3.5, 7.0, 8.5, 36.0,
8.0, 4.5, 6.5, 10.0, 13.0, 41.0, 14.0,
21.5, 10.5, 6.5, 0.0, 15.5, 24.0, 16.0,
0.0, 0.0, 0.0, 4.5, 1.0, 4.0, 6.5,
18.0, 5.0, 7.0, 10.0, 32.5, 18.5, 8.0)
dat <- data.frame(A, B, C, D, y)
model <- lm(y ~ A * B * C * D, data = dat)
effects <- 2 * coef(model)[-1]
effects
## A B C D A:B A:C
## 5.58035714 -3.86607143 -2.41964286 -0.08035714 2.93750000 -0.50892857
## B:C A:D B:D C:D A:B:C A:B:D
## -1.88392857 1.68750000 1.27678571 -0.38392857 -0.65178571 2.15178571
## A:C:D B:C:D A:B:C:D
## -1.04464286 -0.95535714 1.70535714
anova(model)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 871.9 871.93 10.1962 0.001904 **
## B 1 418.5 418.50 4.8939 0.029323 *
## C 1 163.9 163.93 1.9170 0.169401
## D 1 0.2 0.18 0.0021 0.963421
## A:B 1 241.6 241.61 2.8253 0.096039 .
## A:C 1 7.3 7.25 0.0848 0.771515
## B:C 1 99.4 99.38 1.1621 0.283733
## A:D 1 79.7 79.73 0.9324 0.336667
## B:D 1 45.6 45.65 0.5338 0.466808
## C:D 1 4.1 4.13 0.0483 0.826580
## A:B:C 1 11.9 11.90 0.1391 0.710001
## A:B:D 1 129.6 129.65 1.5160 0.221228
## A:C:D 1 30.6 30.56 0.3573 0.551410
## B:C:D 1 25.6 25.56 0.2988 0.585877
## A:B:C:D 1 81.4 81.43 0.9522 0.331604
## Residuals 96 8209.5 85.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, factors A (length of putt) and B (type of putter) are statistically significant. The other factors and interactions are not significant. A longer putt increases the average distance from the cup, while a cavity-back putter decreases it. Therefore, putt length and putter type are the important factors affecting putting performance.
model.reduced <- lm(y ~ A + B, data = dat)
par(mfrow = c(2, 2))
plot(model.reduced)
## hat values (leverages) are all = 0.02678571
## and there are no factor predictors; no plot no. 5
par(mfrow = c(1, 1))
The residual plots show one potentially unusual observation. This may indicate a possible outlier. However, there are no other obvious patterns suggesting serious model inadequacy. The unusual observation should be investigated further to determine whether it resulted from experimental or measurement error.
A <- c(-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1)
B <- rep(c(-1,-1,1,1), 4)
C <- rep(c(-1,-1,-1,-1,1,1,1,1), 2)
D <- rep(c(rep(-1,8), rep(1,8)))
y <- c(1.92, 11.28, 1.09, 5.75,
2.13, 9.53, 1.03, 5.35,
1.60, 11.73, 1.16, 4.68,
2.16, 9.11, 1.07, 5.30)
dat <- data.frame(A, B, C, D, y)
model.full <- lm(y ~ A * B * C * D, data = dat)
library(FrF2)
## Loading required package: DoE.base
## Loading required package: grid
## Loading required package: conf.design
## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design
##
## Attaching package: 'DoE.base'
## The following objects are masked from 'package:stats':
##
## aov, lm
## The following object is masked from 'package:graphics':
##
## plot.design
## The following object is masked from 'package:base':
##
## lengths
DanielPlot(model.full, half = TRUE)
The half normal plot indicates that A, B, and AB are the most important effects. Therefore, the tentative model includes A, B, and their interaction AB.
model <- lm(y ~ A + B + A:B, data = dat)
summary(model)
##
## Call:
## lm.default(formula = y ~ A + B + A:B, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.30250 -0.13125 0.01625 0.18500 1.31750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6806 0.1732 27.028 4.05e-12 ***
## A 3.1606 0.1732 18.251 4.05e-10 ***
## B -1.5019 0.1732 -8.672 1.63e-06 ***
## A:B -1.0694 0.1732 -6.175 4.76e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6927 on 12 degrees of freedom
## Multiple R-squared: 0.9738, Adjusted R-squared: 0.9673
## F-statistic: 148.8 on 3 and 12 DF, p-value: 9.328e-10
par(mfrow = c(2, 2))
plot(model)
## hat values (leverages) are all = 0.25
## and there are no factor predictors; no plot no. 5
par(mfrow = c(1, 1))
The model explains approximately 97.4% of the variation in resistivity. The residual plots should be examined for normality, constant variance, and unusual patterns.
dat$ln_y <- log(dat$y)
model.log.full <- lm(ln_y ~ A * B * C * D, data = dat)
effects.log <- 2 * coef(model.log.full)[-1]
effects.log
## A B C D A:B A:C
## 1.625740689 -0.628555107 -0.012817116 -0.036154779 -0.049369141 -0.079447400
## B:C A:D B:D C:D A:B:C A:B:D
## -0.008451593 -0.019156489 0.007417447 0.035560865 0.126868817 -0.059751921
## A:C:D B:C:D A:B:C:D
## -0.007480469 0.007531521 0.062644086
qqnorm(effects.log)
qqline(effects.log)
model.log <- lm(ln_y ~ A + B + A:B:C, data = dat)
par(mfrow = c(2, 2))
plot(model.log)
## hat values (leverages) are all = 0.25
## and there are no factor predictors; no plot no. 5
par(mfrow = c(1, 1))
The normal probability plot suggests that A, B, and ABC are important effects. The log transformation improves the residual behavior, particularly the variation across fitted values. Therefore, the transformation is useful for improving model adequacy.
coef(model.log)
## (Intercept) A B A:B:C
## 1.18541712 0.81287034 -0.31427755 0.06343441
predict(model.log, newdata = dat)
## 1 2 3 4 5 6
## 0.623389916 2.375999422 0.121703626 1.620575499 0.750258733 2.249130606
## 7 8 9 10 11 12
## -0.005165191 1.747444315 0.623389916 2.375999422 0.121703626 1.620575499
## 13 14 15 16
## 0.750258733 2.249130606 -0.005165191 1.747444315
The fitted regression equation is: \[\ln(\hat y)=1.1854+0.8129A-0.3143B+0.0634ABC\]
Converting back to the original resistivity scale: \[\hat y=\exp(1.1854+0.8129A-0.3143B+0.0634ABC)\] The equation can be used to predict resistivity using the coded factor settings. The logarithmic transformation helps improve the model’s residual behavior.
A <- rep(c(-1, 1), times = 16)
B <- rep(rep(c(-1, 1), each = 2), times = 8)
C <- rep(rep(c(-1, 1), each = 4), times = 4)
D <- rep(rep(c(-1, 1), each = 8), times = 2)
E <- rep(c(-1, 1), each = 16)
y <- c(
8.11, 5.56, 5.77, 5.82,
9.17, 7.80, 3.23, 5.69,
8.82, 14.23, 9.20, 8.94,
8.68, 11.49, 6.25, 9.12,
7.93, 5.00, 7.47, 12.00,
9.86, 3.65, 6.40, 11.61,
12.43, 17.55, 8.87, 25.38,
13.06, 18.85, 11.78, 26.05
)
dat <- data.frame(A, B, C, D, E, y)
library(FrF2)
model.full <- lm(y ~ A * B * C * D * E, data = dat)
DanielPlot(model.full, half = TRUE)
Based on the Daniel half-normal plot, the important effects are Main effects A, D, E; Two-factor interactions AB, AD, AE, BE, DE; Three-factor interactions ABE, ADE.
Factor C and the remaining interactions appear to be insignificant.
model <- lm(y ~ A * B * E + A * D * E, data = dat)
par(mfrow = c(2, 2))
plot(model)
## hat values (leverages) are all = 0.375
## and there are no factor predictors; no plot no. 5
par(mfrow = c(1, 1))
The residual plots show no strong indication of nonconstant variance or serious departures from normality. The model appears to be adequate.
dat4 <- dat[, c("A", "B", "D", "E", "y")]
model4 <- lm(y ~ A * B * D * E, data = dat4)
anova(model4)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 83.560 83.560 57.2328 1.136e-06 ***
## B 1 0.060 0.060 0.0414 0.8414184
## D 1 285.784 285.784 195.7422 2.161e-10 ***
## E 1 153.169 153.169 104.9099 1.966e-08 ***
## A:B 1 48.931 48.931 33.5142 2.767e-05 ***
## A:D 1 88.878 88.878 60.8751 7.661e-07 ***
## B:D 1 0.006 0.006 0.0040 0.9506177
## A:E 1 33.764 33.764 23.1257 0.0001928 ***
## B:E 1 52.711 52.711 36.1032 1.822e-05 ***
## D:E 1 61.799 61.799 42.3283 7.243e-06 ***
## A:B:D 1 3.816 3.816 2.6135 0.1255014
## A:B:E 1 44.959 44.959 30.7937 4.402e-05 ***
## A:D:E 1 26.010 26.010 17.8151 0.0006496 ***
## B:D:E 1 0.050 0.050 0.0345 0.8549347
## A:B:D:E 1 5.306 5.306 3.6340 0.0747350 .
## Residuals 16 23.360 1.460
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
After removing factor C, the experiment becomes a replicated \(2^4\) factorial design with two replicates. The ANOVA confirms that A, D, E, AB, AD, AE, BE, DE, ABE, and ADE are significant at the 0.05 level. The results are consistent with part (a).
settings <- expand.grid(
A = c(-1, 1),
B = c(-1, 1),
D = c(-1, 1),
E = c(-1, 1))
settings$Predicted <- predict(model4, newdata = settings)
settings[which.max(settings$Predicted), ]
## A B D E Predicted
## 16 1 1 1 1 25.715
The maximum predicted response is 25.715, obtained when all four factors are at their high levels (+1).
A: Length of Putt (10 ft = −1, 30 ft = +1)
B: Type of Putter (Mallet = −1, Cavity back = +1)
C: Break of Putt (Straight = −1, Breaking = +1)
D: Slope of Putt (Level = −1, Downhill = +1)
A <- rep(c(-1, 1, -1, 1,
-1, 1, -1, 1,
-1, 1, -1, 1,
-1, 1, -1, 1), each = 7)
B <- rep(c(-1, -1, 1, 1,
-1, -1, 1, 1,
-1, -1, 1, 1,
-1, -1, 1, 1), each = 7)
C <- rep(c(-1, -1, -1, -1,
1, 1, 1, 1,
-1, -1, -1, -1,
1, 1, 1, 1), each = 7)
D <- rep(c(rep(-1, 8), rep(1, 8)), each = 7)
Block <- factor(rep(1:7, times = 16))
y <- c(10.0, 18.0, 14.0, 12.5, 19.0, 16.0, 18.5,
0.0, 16.5, 4.5, 17.5, 20.5, 17.5, 33.0,
4.0, 6.0, 1.0, 14.5, 12.0, 14.0, 5.0,
0.0, 10.0, 34.0, 11.0, 25.5, 21.5, 0.0,
0.0, 0.0, 18.5, 19.5, 16.0, 15.0, 11.0,
5.0, 20.5, 18.0, 20.0, 29.5, 19.0, 10.0,
6.5, 18.5, 7.5, 6.0, 0.0, 10.0, 0.0,
16.5, 4.5, 0.0, 23.5, 8.0, 8.0, 8.0,
4.5, 18.0, 14.5, 10.0, 8.0, 17.5, 6.0,
19.5, 18.0, 16.0, 5.5, 10.0, 7.0, 36.0,
15.0, 16.0, 8.5, 0.0, 0.5, 9.0, 3.0,
41.5, 39.0, 6.5, 3.5, 7.0, 8.5, 36.0,
8.0, 4.5, 6.5, 10.0, 13.0, 41.0, 14.0,
21.5, 10.5, 6.5, 0.0, 15.5, 24.0, 16.0,
0.0, 0.0, 0.0, 4.5, 1.0, 4.0, 6.5,
18.0, 5.0, 7.0, 10.0, 32.5, 18.5, 8.0)
dat <- data.frame(A, B, C, D, Block, y)
model <- lm(y ~ Block + A * B * C * D, data = dat)
anova(model)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 6 392.7 65.45 0.7536 0.608214
## A 1 871.9 871.93 10.0391 0.002094 **
## B 1 418.5 418.50 4.8185 0.030729 *
## C 1 163.9 163.93 1.8874 0.172903
## D 1 0.2 0.18 0.0021 0.963710
## A:B 1 241.6 241.61 2.7818 0.098816 .
## A:C 1 7.3 7.25 0.0835 0.773274
## B:C 1 99.4 99.38 1.1442 0.287628
## A:D 1 79.7 79.73 0.9180 0.340558
## B:D 1 45.6 45.65 0.5255 0.470368
## C:D 1 4.1 4.13 0.0475 0.827931
## A:B:C 1 11.9 11.90 0.1370 0.712196
## A:B:D 1 129.6 129.65 1.4927 0.224990
## A:C:D 1 30.6 30.56 0.3518 0.554578
## B:C:D 1 25.6 25.56 0.2942 0.588856
## A:B:C:D 1 81.4 81.43 0.9376 0.335499
## Residuals 90 7816.8 86.85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, factors A (length of putt) and B (type of putter) are statistically significant. The other factors and interactions are not significant. The block effect is not significant (p = 0.6082), indicating that there is no strong evidence of variation among the seven replicates. Therefore, blocking does not appear to be useful in this experiment. The results are consistent with Problem 6.21, where factors A and B were also significant.
design <- expand.grid(
A = c(-1, 1), B = c(-1, 1),
C = c(-1, 1), D = c(-1, 1),
E = c(-1, 1), F = c(-1, 1))
design$G1 <- with(design, A * B * C * D)
design$G2 <- with(design, A * B * E * F)
design$Block <- interaction(design$G1, design$G2)
design
## A B C D E F G1 G2 Block
## 1 -1 -1 -1 -1 -1 -1 1 1 1.1
## 2 1 -1 -1 -1 -1 -1 -1 -1 -1.-1
## 3 -1 1 -1 -1 -1 -1 -1 -1 -1.-1
## 4 1 1 -1 -1 -1 -1 1 1 1.1
## 5 -1 -1 1 -1 -1 -1 -1 1 -1.1
## 6 1 -1 1 -1 -1 -1 1 -1 1.-1
## 7 -1 1 1 -1 -1 -1 1 -1 1.-1
## 8 1 1 1 -1 -1 -1 -1 1 -1.1
## 9 -1 -1 -1 1 -1 -1 -1 1 -1.1
## 10 1 -1 -1 1 -1 -1 1 -1 1.-1
## 11 -1 1 -1 1 -1 -1 1 -1 1.-1
## 12 1 1 -1 1 -1 -1 -1 1 -1.1
## 13 -1 -1 1 1 -1 -1 1 1 1.1
## 14 1 -1 1 1 -1 -1 -1 -1 -1.-1
## 15 -1 1 1 1 -1 -1 -1 -1 -1.-1
## 16 1 1 1 1 -1 -1 1 1 1.1
## 17 -1 -1 -1 -1 1 -1 1 -1 1.-1
## 18 1 -1 -1 -1 1 -1 -1 1 -1.1
## 19 -1 1 -1 -1 1 -1 -1 1 -1.1
## 20 1 1 -1 -1 1 -1 1 -1 1.-1
## 21 -1 -1 1 -1 1 -1 -1 -1 -1.-1
## 22 1 -1 1 -1 1 -1 1 1 1.1
## 23 -1 1 1 -1 1 -1 1 1 1.1
## 24 1 1 1 -1 1 -1 -1 -1 -1.-1
## 25 -1 -1 -1 1 1 -1 -1 -1 -1.-1
## 26 1 -1 -1 1 1 -1 1 1 1.1
## 27 -1 1 -1 1 1 -1 1 1 1.1
## 28 1 1 -1 1 1 -1 -1 -1 -1.-1
## 29 -1 -1 1 1 1 -1 1 -1 1.-1
## 30 1 -1 1 1 1 -1 -1 1 -1.1
## 31 -1 1 1 1 1 -1 -1 1 -1.1
## 32 1 1 1 1 1 -1 1 -1 1.-1
## 33 -1 -1 -1 -1 -1 1 1 -1 1.-1
## 34 1 -1 -1 -1 -1 1 -1 1 -1.1
## 35 -1 1 -1 -1 -1 1 -1 1 -1.1
## 36 1 1 -1 -1 -1 1 1 -1 1.-1
## 37 -1 -1 1 -1 -1 1 -1 -1 -1.-1
## 38 1 -1 1 -1 -1 1 1 1 1.1
## 39 -1 1 1 -1 -1 1 1 1 1.1
## 40 1 1 1 -1 -1 1 -1 -1 -1.-1
## 41 -1 -1 -1 1 -1 1 -1 -1 -1.-1
## 42 1 -1 -1 1 -1 1 1 1 1.1
## 43 -1 1 -1 1 -1 1 1 1 1.1
## 44 1 1 -1 1 -1 1 -1 -1 -1.-1
## 45 -1 -1 1 1 -1 1 1 -1 1.-1
## 46 1 -1 1 1 -1 1 -1 1 -1.1
## 47 -1 1 1 1 -1 1 -1 1 -1.1
## 48 1 1 1 1 -1 1 1 -1 1.-1
## 49 -1 -1 -1 -1 1 1 1 1 1.1
## 50 1 -1 -1 -1 1 1 -1 -1 -1.-1
## 51 -1 1 -1 -1 1 1 -1 -1 -1.-1
## 52 1 1 -1 -1 1 1 1 1 1.1
## 53 -1 -1 1 -1 1 1 -1 1 -1.1
## 54 1 -1 1 -1 1 1 1 -1 1.-1
## 55 -1 1 1 -1 1 1 1 -1 1.-1
## 56 1 1 1 -1 1 1 -1 1 -1.1
## 57 -1 -1 -1 1 1 1 -1 1 -1.1
## 58 1 -1 -1 1 1 1 1 -1 1.-1
## 59 -1 1 -1 1 1 1 1 -1 1.-1
## 60 1 1 -1 1 1 1 -1 1 -1.1
## 61 -1 -1 1 1 1 1 1 1 1.1
## 62 1 -1 1 1 1 1 -1 -1 -1.-1
## 63 -1 1 1 1 1 1 -1 -1 -1.-1
## 64 1 1 1 1 1 1 1 1 1.1
table(design$Block)
##
## -1.-1 1.-1 -1.1 1.1
## 16 16 16 16
The 64 experimental runs are divided into four blocks, with 16 runs in each block. The effects confounded with blocks are ABCD, ABEF, and CDEF.
design <- expand.grid(
A = c(-1, 1), B = c(-1, 1),
C = c(-1, 1), D = c(-1, 1),
E = c(-1, 1), F = c(-1, 1))
design$G1 <- with(design, A*B*C*D)
design$G2 <- with(design, A*C*E)
design$G3 <- with(design, A*B*E*F)
design$Block <- interaction(design$G1, design$G2, design$G3)
design
## A B C D E F G1 G2 G3 Block
## 1 -1 -1 -1 -1 -1 -1 1 -1 1 1.-1.1
## 2 1 -1 -1 -1 -1 -1 -1 1 -1 -1.1.-1
## 3 -1 1 -1 -1 -1 -1 -1 -1 -1 -1.-1.-1
## 4 1 1 -1 -1 -1 -1 1 1 1 1.1.1
## 5 -1 -1 1 -1 -1 -1 -1 1 1 -1.1.1
## 6 1 -1 1 -1 -1 -1 1 -1 -1 1.-1.-1
## 7 -1 1 1 -1 -1 -1 1 1 -1 1.1.-1
## 8 1 1 1 -1 -1 -1 -1 -1 1 -1.-1.1
## 9 -1 -1 -1 1 -1 -1 -1 -1 1 -1.-1.1
## 10 1 -1 -1 1 -1 -1 1 1 -1 1.1.-1
## 11 -1 1 -1 1 -1 -1 1 -1 -1 1.-1.-1
## 12 1 1 -1 1 -1 -1 -1 1 1 -1.1.1
## 13 -1 -1 1 1 -1 -1 1 1 1 1.1.1
## 14 1 -1 1 1 -1 -1 -1 -1 -1 -1.-1.-1
## 15 -1 1 1 1 -1 -1 -1 1 -1 -1.1.-1
## 16 1 1 1 1 -1 -1 1 -1 1 1.-1.1
## 17 -1 -1 -1 -1 1 -1 1 1 -1 1.1.-1
## 18 1 -1 -1 -1 1 -1 -1 -1 1 -1.-1.1
## 19 -1 1 -1 -1 1 -1 -1 1 1 -1.1.1
## 20 1 1 -1 -1 1 -1 1 -1 -1 1.-1.-1
## 21 -1 -1 1 -1 1 -1 -1 -1 -1 -1.-1.-1
## 22 1 -1 1 -1 1 -1 1 1 1 1.1.1
## 23 -1 1 1 -1 1 -1 1 -1 1 1.-1.1
## 24 1 1 1 -1 1 -1 -1 1 -1 -1.1.-1
## 25 -1 -1 -1 1 1 -1 -1 1 -1 -1.1.-1
## 26 1 -1 -1 1 1 -1 1 -1 1 1.-1.1
## 27 -1 1 -1 1 1 -1 1 1 1 1.1.1
## 28 1 1 -1 1 1 -1 -1 -1 -1 -1.-1.-1
## 29 -1 -1 1 1 1 -1 1 -1 -1 1.-1.-1
## 30 1 -1 1 1 1 -1 -1 1 1 -1.1.1
## 31 -1 1 1 1 1 -1 -1 -1 1 -1.-1.1
## 32 1 1 1 1 1 -1 1 1 -1 1.1.-1
## 33 -1 -1 -1 -1 -1 1 1 -1 -1 1.-1.-1
## 34 1 -1 -1 -1 -1 1 -1 1 1 -1.1.1
## 35 -1 1 -1 -1 -1 1 -1 -1 1 -1.-1.1
## 36 1 1 -1 -1 -1 1 1 1 -1 1.1.-1
## 37 -1 -1 1 -1 -1 1 -1 1 -1 -1.1.-1
## 38 1 -1 1 -1 -1 1 1 -1 1 1.-1.1
## 39 -1 1 1 -1 -1 1 1 1 1 1.1.1
## 40 1 1 1 -1 -1 1 -1 -1 -1 -1.-1.-1
## 41 -1 -1 -1 1 -1 1 -1 -1 -1 -1.-1.-1
## 42 1 -1 -1 1 -1 1 1 1 1 1.1.1
## 43 -1 1 -1 1 -1 1 1 -1 1 1.-1.1
## 44 1 1 -1 1 -1 1 -1 1 -1 -1.1.-1
## 45 -1 -1 1 1 -1 1 1 1 -1 1.1.-1
## 46 1 -1 1 1 -1 1 -1 -1 1 -1.-1.1
## 47 -1 1 1 1 -1 1 -1 1 1 -1.1.1
## 48 1 1 1 1 -1 1 1 -1 -1 1.-1.-1
## 49 -1 -1 -1 -1 1 1 1 1 1 1.1.1
## 50 1 -1 -1 -1 1 1 -1 -1 -1 -1.-1.-1
## 51 -1 1 -1 -1 1 1 -1 1 -1 -1.1.-1
## 52 1 1 -1 -1 1 1 1 -1 1 1.-1.1
## 53 -1 -1 1 -1 1 1 -1 -1 1 -1.-1.1
## 54 1 -1 1 -1 1 1 1 1 -1 1.1.-1
## 55 -1 1 1 -1 1 1 1 -1 -1 1.-1.-1
## 56 1 1 1 -1 1 1 -1 1 1 -1.1.1
## 57 -1 -1 -1 1 1 1 -1 1 1 -1.1.1
## 58 1 -1 -1 1 1 1 1 -1 -1 1.-1.-1
## 59 -1 1 -1 1 1 1 1 1 -1 1.1.-1
## 60 1 1 -1 1 1 1 -1 -1 1 -1.-1.1
## 61 -1 -1 1 1 1 1 1 -1 1 1.-1.1
## 62 1 -1 1 1 1 1 -1 1 -1 -1.1.-1
## 63 -1 1 1 1 1 1 -1 -1 -1 -1.-1.-1
## 64 1 1 1 1 1 1 1 1 1 1.1.1
table(design$Block)
##
## -1.-1.-1 1.-1.-1 -1.1.-1 1.1.-1 -1.-1.1 1.-1.1 -1.1.1 1.1.1
## 8 8 8 8 8 8 8 8
The 64 runs are divided into 8 blocks, with 8 runs in each block.
Find the other effects confounded with blocks
\[G_1G_2=BDE\]
\[G_1G_3=CDEF\]
\[G_2G_3=BCF\]
\[G_1G_2G_3=ADF\]
The seven effects confounded with blocks are: \(ABCD,\ ACE,\ ABEF,\ BDE,\ CDEF,\ BCF,\ ADF\)These effects are confounded with blocks and cannot be independently estimated.
library(GAD)
Position <- rep(c(1, 2), each = 3, times = 3)
Temperature <- rep(c(800, 825, 850), each = 6)
y <- c(
570, 565, 583, 528, 547, 521,
1063, 1080, 1043, 988, 1026, 1004,
565, 510, 590, 526, 538, 532)
Position <- as.random(Position)
Temperature <- as.fixed(Temperature)
dat <- data.frame(Temperature, Position, y)
model <- aov(y ~ Temperature + Position, data = dat)
GAD::gad(model)
## $anova
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 945342 472671 1069.257 4.924e-16 ***
## Position 1 7160 7160 16.197 0.001254 **
## Residuals 14 6189 442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At the 0.05 significance level, temperature has a significant effect on baked density. The random position effect is also significant, indicating variation among furnace positions.
MS_Position <- 7160
MS_Error <- 442
sigma_error <- MS_Error
sigma_position <- (MS_Position - MS_Error) / 9
sigma_error
## [1] 442
sigma_position
## [1] 746.4444
The estimated variance component for furnace position is 746.444, indicating substantial variability due to position
library(GAD)
Part <- rep(1:10, each = 6)
Operator <- rep(rep(1:2, each = 3), times = 10)
y <- c(50,49,50, 50,48,51,
52,52,51, 51,51,51,
53,50,50, 54,52,51,
49,51,50, 48,50,51,
48,49,48, 48,49,48,
52,50,50, 52,50,50,
51,51,51, 51,50,50,
52,50,49, 53,48,50,
50,51,50, 51,48,49,
47,46,49, 46,47,48)
Part <- as.random(Part)
Operator <- as.fixed(Operator)
model <- aov(y ~ Part * Operator)
GAD::gad(model)
## $anova
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Part 9 99.017 11.0019 7.3346 3.216e-06 ***
## Operator 1 0.417 0.4167 0.6923 0.4269
## Part:Operator 9 5.417 0.6019 0.4012 0.9270
## Residuals 40 60.000 1.5000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MS_Part <- 11.0019
MS_Interaction <- 0.6019
MS_Error <- 1.5000
sigma_part <- (MS_Part - MS_Interaction) / (2 * 3)
sigma_interaction <- (MS_Interaction - MS_Error) / 3
sigma_error <- MS_Error
sigma_part
## [1] 1.733333
sigma_interaction
## [1] -0.2993667
sigma_error
## [1] 1.5
The part variance component is 1.7333, while the interaction variance component is -0.2994 (approximately zero). The operator effect is not significant (p = 0.4269), but the part effect is significant (p < 0.001).