homework5
Oghenetega Dibie
4/5/2020
Chapter 6: 1, 2, 3, 4, 5, 6, 7, 9, 10, 12, 13, 21
#1
#df of model = 2^k - 1 = 15
# 2 levels each with 4 replicates.
# f- none of above#2
# 2^3 factorial is replicated twice
# 3 factors- with two levels each
# replicated twice
# df_total = 2 * 2^3 - 1= 15
#dfmodel= maineffect + 2way_effect + 3way_effect
#dfmodel = 3 + 3 + 1
#dferror = dftotal - dfmodel= 15-7 -8#3
# 2^3 factorial is replicated twice
# 3 factors- with two levels each
# replicated twice
# df_total = 2 * 2^3 - 1= 15
#dfmodel= maineffect
#dfmodel = 3 + 3 + 1
#dferror = dftotal - dfmodel= 15-3=12#4
# 2^3 factorial is replicated thrice
# 3 factors- with two levels each
# replicated three times
# df_total = 3 * 2^3= 23
#dfmodel= maineffect
#dfmodel = 3 + 1 _ 1
#dferror = dftotal - dfmodel= 23- 5 = 18#5
#a
y <- c(22,31,25,
32,43,29,
35,34,50,
55,47,46,
44,45,38,
40,37,36,
60,50,54,
39,41,47)
A <- rep(c(rep("-",3),rep("+",3)),4)
B <- rep(c(rep("-",6),rep("+",6)),2)
C <- c(rep("-",12),rep("+",12))
ym <- matrix(y,ncol=3,byrow=T)
tot <- rowSums(ym)
a1 <- rep(c(-1,1),4)
b1 <- rep(c(-1,-1,1,1),2)
c1 <- c(rep(-1,4),rep(1,4))
sum(a1*tot)## [1] 4sum(b1*tot)## [1] 136sum(c1*tot)## [1] 82bin <- cbind(a1,b1,c1)
contr <- tot%*%bin
factor_effect <- contr/(2^(3-1)*3)
df <- data.frame(y,A,B,C)
m <- aov(y~A+B+C)
anova(m)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 0.67 0.67 0.0128 0.911142
## B 1 770.67 770.67 14.7661 0.001016 **
## C 1 280.17 280.17 5.3680 0.031232 *
## Residuals 20 1043.83 52.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A <- rep(c(-1,1),4)
B <- rep(c(-1,-1,1,1),2)
C <- c(rep(-1,4),rep(1,4))
bin <- cbind(A,B,C)
colnames(bin) <- c("A","B","C")
# (b)
SS <- contr^2/(2^3*3)
SS## a1 b1 c1
## [1,] 0.6666667 770.6667 280.1667SST <- sum(SS)
SStot <- sum((y-mean(y))^2)
SSE <- SStot-SST
## deg freedom
dftot <- 2^3*3-1
dftr <- 3
dferr <- dftot-dftr
MSE <- SSE/dferr
F <- (SS/1)/MSE
1-pf(F,1,dferr)## a1 b1 c1
## [1,] 0.9111418 0.001016122 0.03123236#FE
ybar <- mean(y)
#C
a <- rep(c(rep(-1,3),rep(1,3)),4)
b <- rep(c(rep(-1,6),rep(1,6)),2)
c <- c(rep(-1,12),rep(1,12))
x1 <- 2*(a-(max(a)+min(a))/2)/(max(a)-min(a))
x2 <- 2*(b-(max(b)+min(b))/2)/(max(b)-min(b))
x3 <- 2*(c-(max(c)+min(c))/2)/(max(c)-min(c))
summary(lm(y~x1*x2*x3))##
## Call:
## lm(formula = y ~ x1 * x2 * x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.667 -3.500 -1.167 3.167 10.333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.8333 1.1211 36.421 < 2e-16 ***
## x1 0.1667 1.1211 0.149 0.883680
## x2 5.6667 1.1211 5.054 0.000117 ***
## x3 3.4167 1.1211 3.048 0.007679 **
## x1:x2 -0.8333 1.1211 -0.743 0.468078
## x1:x3 -4.4167 1.1211 -3.939 0.001172 **
## x2:x3 -1.4167 1.1211 -1.264 0.224475
## x1:x2:x3 -1.0833 1.1211 -0.966 0.348282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.492 on 16 degrees of freedom
## Multiple R-squared: 0.7696, Adjusted R-squared: 0.6689
## F-statistic: 7.637 on 7 and 16 DF, p-value: 0.0003977## The regression model is
# y=40.8333+5.6667x2+3.4167x3-4.4167x13
m <- lm(y~x1*x2*x3)
#d
plot(m$res,m$pred)
abline(h=0)
qqnorm(m$res)
qqline(m$res)
shapiro.test(m$res)##
## Shapiro-Wilk normality test
##
## data: m$res
## W = 0.92106, p-value = 0.06166# The normality test passes but the fitted vs residuals plot has a pattern which means equal variance test fails
#e
#Since B has a positive effect, set B at the high level to increase life. The AC interaction plot reveals that life would be maximized with C at the high level and A at the low level.
#6
#Result provided in book
#7
ngrid <- 25
x <- seq(-1,1,length.out=ngrid)
z1 <- matrix(0,ngrid,ngrid)
for(i in 1:ngrid){
for(j in 1:ngrid){
z1[i,j] <- predict(m,newdata=data.frame(x1=1,x2=x[i],x3=x[j],x13=x[j]))
}
}
persp(x,x,z1,theta=-30,phi=15)
z2 <- matrix(0,ngrid,ngrid)
for(i in 1:ngrid){
for(j in 1:ngrid){
z2[i,j] <- predict(m,newdata=data.frame(x1=-1,x2=x[i],x3=x[j],x13=-x[j]))
}
}
persp(x,x,z2,theta=-60,phi=15)
#9
y <- c(18.2,18.9,12.9,14.4,
27.2,24.0,22.4,22.5,
15.9,14.5,15.1,14.2,
41.0,43.9,36.3,39.9)
A <- rep(c(rep("-",4),rep("+",4)),2)
B <- c(rep("-",8),rep("+",8))
df <- data.frame(y,A,B)
#fit the anova model
model <- aov(y~A*B)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 1107.23 1107.23 185.252 1.175e-08 ***
## B 1 227.26 227.26 38.023 4.826e-05 ***
## A:B 1 303.63 303.63 50.801 1.201e-05 ***
## Residuals 12 71.72 5.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Main effect and interaction effect are significant
qqnorm(model$residuals)
qqline(model$residuals)
plot(model$fitted.values,model$residuals)
abline(h=0)
#Normality and equal variance test holds
## interaction plot
interaction.plot(x.factor=A,trace.factor=B,response=y,
fun=mean,type="b",col=c(2,3),
pch=c(19,15),fixed=TRUE,leg.bty="o",main="Interaction graph")
#To reduce the vibration, use the smaller bit. Once the small bit is specified, either speed will work equally well, because the slope of the curve relating vibration to speed for the small tip is approximately zero. The process is robust to speed changes if the small bit is used -- Answer provided in book.10
A <- rep(c(-1,1),4)
B <- rep(c(-1,-1,1,1),2)
C <- c(rep(-1,4),rep(1,4))
y <- c(22,32,35,55,44,40,60,39)
bin <- cbind(A,B,C)
y <- c(y,36,40,43,45)
k <- 3
nc <- 4
bin <- rbind(bin,matrix(0,nc,k))
df <- data.frame(y,bin)
model <- aov(y~I (A^2)+A*B*C,data=df)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## I(A^2) 1 0.04 0.04 0.0027 0.96170
## A 1 3.12 3.12 0.2038 0.68231
## B 1 325.12 325.12 21.2038 0.01926 *
## C 1 190.12 190.12 12.3995 0.03888 *
## A:B 1 6.12 6.12 0.3995 0.57225
## A:C 1 378.12 378.12 24.6603 0.01568 *
## B:C 1 55.12 55.12 3.5951 0.15423
## A:B:C 1 91.13 91.13 5.9429 0.09268 .
## Residuals 3 46.00 15.33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1qqnorm(model$residuals)
qqline(model$residuals)
plot(model$fitted.values,model$residuals)
#Effects B, C and AC are significant at 5%. There is no effect of curvature.
#The model does not differ substantially from the model in Problem 6.5, part (c).
#(c) # y = 40.84+(0.33/2)*A+(11.33/2)*B+(6.83/2)*C
#d Residuals plot seem to indicate normality.
#e. Solution provided in book12
# response varaible
y <- c(21,23,20,22,28,26,
37,38,35,39,38,36,
25,24,29,26,25,27,
31,29,30,34,33,35)
A <- rep(c(rep(-1,6),rep(1,6)),2)
B <- c(rep(-1,12),rep(1,12))
df <- data.frame(y,A,B)
model <- aov(y~A*B)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 590.04 590.04 115.5057 9.291e-10 ***
## B 1 9.37 9.37 1.8352 0.1906172
## A:B 1 92.04 92.04 18.0179 0.0003969 ***
## Residuals 20 102.17 5.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# factor A and AB are significant.
# residuals
qqnorm(model$residuals)
qqline(model$residuals)
plot(model$fitted.values,model$residuals)
abline(h=0)
#Assumption holds: normality and equality of variance
# interaction plot
interaction.plot(x.factor=A,trace.factor=B,response=y,
fun=mean,type="b",col=c(2,3),
pch=c(19,15),fixed=TRUE,leg.bty="o",main="Interaction graph")
13
y <- c(5.12,4.98,4.89,5.00,
4.95,4.27,4.43,4.25,
6.65,5.49,6.24,5.55,
5.28,4.75,4.91,4.71)
A <- rep(c(rep("Glass",4),rep("Plastic",4)),2)
B <- c(rep("1",8),rep("2",8))
df <- data.frame(y,A,B)
model <- aov(y~A*B)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 2.53606 2.53606 20.4092 0.0007047 ***
## B 1 2.02351 2.02351 16.2844 0.0016531 **
## A:B 1 0.29976 0.29976 2.4123 0.1463480
## Residuals 12 1.49112 0.12426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Main effect is significant
# residuals
qqnorm(model$residuals)
qqline(model$residuals)
plot(model$fitted.values,model$residuals)
abline(h=0)
# Assumptions holds21
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, 0.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)
A <- rep(c(rep(10,7),rep(30,7)),8)
B <- rep(c(rep("Mallet",14),rep("Cavity back",14)),4)
C <- rep(c(rep("Straight",28),rep("Breaking",28)),2)
D <- c(rep("Level",56),rep("Downhill",56))
df<- data.frame(y,A,B,C,D)
# fit the anova model
model<- aov(y~A*B*C*D)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 917.1 917.15 10.5878 0.001572 **
## B 1 388.1 388.15 4.4809 0.036862 *
## C 1 145.1 145.15 1.6756 0.198615
## D 1 1.4 1.40 0.0161 0.899280
## A:B 1 218.7 218.68 2.5245 0.115377
## A:C 1 11.9 11.90 0.1373 0.711776
## B:C 1 115.0 115.02 1.3278 0.252054
## A:D 1 93.8 93.81 1.0829 0.300658
## B:D 1 56.4 56.43 0.6515 0.421588
## C:D 1 1.6 1.63 0.0188 0.891271
## A:B:C 1 7.3 7.25 0.0837 0.772939
## A:B:D 1 113.0 113.00 1.3045 0.256228
## A:C:D 1 39.5 39.48 0.4558 0.501207
## B:C:D 1 33.8 33.77 0.3899 0.533858
## A:B:C:D 1 95.6 95.65 1.1042 0.295994
## Residuals 96 8315.8 86.62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# factor A and B are significant
# Rerun the mode with significant terms.
model2 <- aov(y~A+B)
anova(model2)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 917.1 917.15 10.8087 0.00136 **
## B 1 388.1 388.15 4.5743 0.03469 *
## Residuals 109 9248.9 84.85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1qqnorm(model$residuals)
qqline(model$residuals)
plot(model$fitted.values,model$residuals)
abline(h=0)
qqnorm(model2$residuals)
qqline(model2$residuals)
plot(model2$fitted.values,model2$residuals)
abline(h=0)
#Assumption holds for both model (Normality and Equality of Variance)
interaction.plot(x.factor=A,trace.factor=B,response=y,
fun=mean,type="b",col=c(2,3),
pch=c(19,15),fixed=TRUE,leg.bty="o",main="Interaction graph")