Problem 7.12
A <- rep(rep(c(-1,1),8),7) # Factor A
B <- rep(rep(c(1,1,-1,-1),4),7) # Factor B
C <- rep(rep(c(rep(1,4),rep(-1,4)),2),7) # Factor C
D <- rep(c(rep(1,8),rep(-1,8)),7) # Factor D
I <- c(10.0,0.0,4.0,0.0,0.0,5.0,6.5,16.5,4.5,19.5,15.0,41.5,8.0,21.5,0.0,18.0)
II <- c(18.0,16.5,6.0,10.0,0.0,20.5,18.5,4.5,18.0,18.0,16.0,39.0,4.5,10.5,0.0,5.0)
III <- c(14.0,4.5,1.0,34.0,18.5,18.0,7.5,0.0,14.5,16.0,8.5,6.5,6.5,6.5,0.0,7.0)
IV <- c(12.5,17.5,14.5,11.0,19.5,20.0,6.0,23.5,10.0,5.5,0.0,3.5,10.0,0.0,4.5,10.0)
V <- c(19.0,20.5,12.0,25.5,16.0,29.5,0.0,8.0,0.0,10.0,0.5,7.0,13.0,15.5,1.0,32.5)
VI <- c(16.0,17.5,14.0,21.5,15.0,19.0,10.0,8.0,17.5,7.0,9.0,8.5,41.0,24.0,4.0,18.5)
VII <- c(18.5,33.0,5.0,0.0,11.0,10.0,0.0,8.0,6.0,36.0,3.0,36.0,14.0,16.0,6.5,8.0)
obs <- c(I,II,III,IV,V,VI,VII)
block <- c(rep(1,16),rep(2,16),rep(3,16),rep(4,16),rep(5,16),rep(6,16),rep(7,16))
dat1 <- cbind(A,B,C,D,block,obs)
dat1 <- as.data.frame(dat1)
str(dat1)
## 'data.frame': 112 obs. of 6 variables:
## $ A : num -1 1 -1 1 -1 1 -1 1 -1 1 ...
## $ B : num 1 1 -1 -1 1 1 -1 -1 1 1 ...
## $ C : num 1 1 1 1 -1 -1 -1 -1 1 1 ...
## $ D : num 1 1 1 1 1 1 1 1 -1 -1 ...
## $ block: num 1 1 1 1 1 1 1 1 1 1 ...
## $ obs : num 10 0 4 0 0 5 6.5 16.5 4.5 19.5 ...
library(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
mod <- lm(obs~A*B*C*D,data = dat1)
coef(mod)
## (Intercept) A B C D A:B
## 12.2991071 2.8616071 1.8616071 1.1383929 0.1116071 -1.3973214
## A:C B:C A:D B:D C:D A:B:C
## 0.3258929 -1.0133929 -0.9151786 0.7098214 -0.1205357 -0.2544643
## A:B:D A:C:D B:C:D A:B:C:D
## 1.0044643 -0.5937500 0.5491071 -0.9241071
halfnormal(mod)
## Warning in halfnormal.lm(mod): halfnormal not recommended for models with more
## residual df than model df
##
## Significant effects (alpha=0.05, Lenth method):
## [1] A e74 e92 B e53

A, B are found to be significant.
Supplymentary: We also try to take block into the linear model to make conclusion more safe (although this is not necessary). But noticeably, block here is a virtue factor, and normally should not be taken into account.
mod <- lm(obs~A*B*C*D*block,data = dat1)
coef(mod)
## (Intercept) A B C D
## 9.968750000 2.852678571 -2.147321429 2.334821429 -1.441964286
## block A:B A:C B:C A:D
## 0.582589286 -1.531250000 -1.245535714 -1.370535714 -2.183035714
## B:D C:D A:block B:block C:block
## 0.477678571 -4.147321429 0.002232143 1.002232143 -0.299107143
## D:block A:B:C A:B:D A:C:D B:C:D
## 0.388392857 -2.933035714 0.504464286 -0.424107143 1.611607143
## A:B:block A:C:block B:C:block A:D:block B:D:block
## 0.033482143 0.392857143 0.089285714 0.316964286 0.058035714
## C:D:block A:B:C:D A:B:C:block A:B:D:block A:C:D:block
## 1.006696429 -1.397321429 0.669642857 0.125000000 -0.042410714
## B:C:D:block A:B:C:D:block
## -0.265625000 0.118303571
halfnormal(mod)
## Warning in halfnormal.lm(mod): halfnormal not recommended for models with more
## residual df than model df
##
## Creation of Ablock
## Projected out: A
##
## Creation of Bblock
## Projected out: B
##
## Creation of Cblock
## Projected out: C
##
## Creation of Dblock
## Projected out: D
##
## Creation of ABblock
## Projected out: AB
##
## Creation of ACblock
## Projected out: AC
##
## Creation of BCblock
## Projected out: BC
##
## Creation of ADblock
## Projected out: AD
##
## Creation of BDblock
## Projected out: BD
##
## Creation of CDblock
## Projected out: CD
##
## Creation of ABCblock
## Projected out: ABC
##
## Creation of ABDblock
## Projected out: ABD
##
## Creation of ACDblock
## Projected out: ACD
##
## Creation of BCDblock
## Projected out: BCD
##
## Creation of ABCDblock
## Projected out: ABCD
##
## Significant effects (alpha=0.05, Lenth method):
## [1] A block lof6 lof72 C:D:block B:block lof45
##
## [8] B lof66

\(\overline{y_{ijklmn}}=\mu+\alpha_{i}+\beta_{j}+\gamma_{k}+\delta_{l}+\zeta_{m}+\alpha\beta_{ij}+\alpha\gamma_{ik}+\alpha\delta_{il}+\beta\gamma_{jk}+\beta\delta_{jl}+\gamma\delta_{kl}+\alpha\beta\gamma_{ijk}+\alpha\beta\delta_{ijl}+\alpha\gamma\delta_{ikl}+\beta\gamma\delta_{jkl}+\alpha\beta\gamma\delta_{ijkl}+\epsilon_{ijklmn}\)
Linear equation then becomes as follow:
\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\epsilon_{ijk}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\), all j
\(H_{1}\): \(\beta_{j}\neq0\), some j
We do not need to actually investigate the block, which is a virtue concept. Since P(A)=0.00136<0.05, P(B)=0.03469<0.05, we confirmed the halfnorm plot that A and B do significantly affect means and we reject the null hypothesis.
Model adequency: From residual and qqnorm plot, It has rough normal distribution and equal variance, even if there is a slight violence on equal variance.
\(\overline{y_{ijklmn}}=\mu+\alpha_{i}+\beta_{j}+\gamma_{k}+\delta_{l}+\zeta_{m}+\alpha\beta_{ij}+\alpha\gamma_{ik}+\alpha\delta_{il}+\beta\gamma_{jk}+\beta\delta_{jl}+\gamma\delta_{kl}+\alpha\beta\gamma_{ijk}+\alpha\beta\delta_{ijl}+\alpha\gamma\delta_{ikl}+\beta\gamma\delta_{jkl}+\alpha\beta\gamma\delta_{ijkl}+\epsilon_{ijklmn}\)
P(A)=0.00176<0.05, P(B)= 0.03875<0.05, and rest of those P-values are >0.05. So we confirmed the abovementioned conclusion that A and B do significantly affect means and reject null hypothesis.
Model adequency: From residual and qqnorm plot, It has rough normal distribution and equal variance.
————————————————————————————————
Raw code:
A <- rep(rep(c(-1,1),8),7) # Factor A
B <- rep(rep(c(1,1,-1,-1),4),7) # Factor B
C <- rep(rep(c(rep(1,4),rep(-1,4)),2),7) # Factor C
D <- rep(c(rep(1,8),rep(-1,8)),7) # Factor D
I <- c(10.0,0.0,4.0,0.0,0.0,5.0,6.5,16.5,4.5,19.5,15.0,41.5,8.0,21.5,0.0,18.0)
II <- c(18.0,16.5,6.0,10.0,0.0,20.5,18.5,4.5,18.0,18.0,16.0,39.0,4.5,10.5,0.0,5.0)
III <- c(14.0,4.5,1.0,34.0,18.5,18.0,7.5,0.0,14.5,16.0,8.5,6.5,6.5,6.5,0.0,7.0)
IV <- c(12.5,17.5,14.5,11.0,19.5,20.0,6.0,23.5,10.0,5.5,0.0,3.5,10.0,0.0,4.5,10.0)
V <- c(19.0,20.5,12.0,25.5,16.0,29.5,0.0,8.0,0.0,10.0,0.5,7.0,13.0,15.5,1.0,32.5)
VI <- c(16.0,17.5,14.0,21.5,15.0,19.0,10.0,8.0,17.5,7.0,9.0,8.5,41.0,24.0,4.0,18.5)
VII <- c(18.5,33.0,5.0,0.0,11.0,10.0,0.0,8.0,6.0,36.0,3.0,36.0,14.0,16.0,6.5,8.0)
obs <- c(I,II,III,IV,V,VI,VII)
block <- c(rep(1,16),rep(2,16),rep(3,16),rep(4,16),rep(5,16),rep(6,16),rep(7,16))
dat1 <- cbind(A,B,C,D,block,obs)
dat1 <- as.data.frame(dat1)
str(dat1)
library(DoE.base)
mod <- lm(obs~A*B*C*D*block,data = dat1)
coef(mod)
halfnormal(mod)
dat1$A <- as.factor(dat1$A)
dat1$B <- as.factor(dat1$B)
dat1$C <- as.factor(dat1$C)
dat1$D <- as.factor(dat1$D)
dat1$block <- as.factor(dat1$block)
model<-aov(obs~A+B,data=dat1)
summary(model)
plot(model)
dat1$A <- as.factor(dat1$A)
dat1$B <- as.factor(dat1$B)
dat1$C <- as.factor(dat1$C)
dat1$D <- as.factor(dat1$D)
dat1$block <- as.factor(dat1$block)
model<-aov(obs~A+B+C+D+block+A*B*C*D,data=dat1)
summary(model)
plot(model)