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

Block here is a virtue factor, even if it shows that C:D:block, B:block, block, A, B are significant. So we only keep A and B.

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)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## A             1    917   917.1  10.809 0.00136 **
## B             1    388   388.1   4.574 0.03469 * 
## Residuals   109   9249    84.9                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(model)

\(\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.

Supplymentary: In order to make our conclusion more safe, we run a full-scale-factor ANOVA to further confirm our conclusion.

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)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## A            1    917   917.1  10.396 0.00176 **
## B            1    388   388.1   4.400 0.03875 * 
## C            1    145   145.1   1.645 0.20290   
## D            1      1     1.4   0.016 0.90021   
## block        6    376    62.7   0.710 0.64202   
## A:B          1    219   218.7   2.479 0.11890   
## A:C          1     12    11.9   0.135 0.71433   
## B:C          1    115   115.0   1.304 0.25655   
## A:D          1     94    93.8   1.063 0.30522   
## B:D          1     56    56.4   0.640 0.42594   
## C:D          1      2     1.6   0.018 0.89227   
## A:B:C        1      7     7.3   0.082 0.77499   
## A:B:D        1    113   113.0   1.281 0.26073   
## A:C:D        1     39    39.5   0.448 0.50520   
## B:C:D        1     34    33.8   0.383 0.53767   
## A:B:C:D      1     96    95.6   1.084 0.30055   
## Residuals   90   7940    88.2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(model)

\(\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.

————————————————————————————————

Problem 7.20

Different from Table 7.8 and Table 7.9, we use ABCD, CDEF as the block generator, so Interactions Confounded with Blocks are, ABCD, CDEF and ABEF.

I did this in excel file, this will be attached in the corresponding pdf profile at the end. 7.20 part is marked in GREEN in the EXCEL file.

Final four blocks list (we assuming we put them in random order):

block 1 block 2 block 3 block 4

“1” a d e

c b ab f

bc bd ac ad

bf be af ae

ef df ce cd

abc abd cf de

abf abe ade acf

acd ace adf aef

bef bdf bce bcd

cde cdf bcf bde

def cef abde abcf

abcd abce abdf abef

acde acdf bcde bcdf

adef acef bdef bcef

abdef acdef cdef abcde

abcdef bcdef abcef abcdf

————————————————————————————————

Problem 7.21

Since we consider the 2^6 design in eight blocks of eight runs each with ABCD, ACE, and ABEF as the independent effects chosen to be confounded with blocks.

Interactions Confounded with Blocks are, ABCD, ACE, ABEF as well as BDE, CDEF, BCF and ADF. 7.21 part is marked in ORANGE in the EXCEL file.

Final four blocks list (we assuming we put them in random order):

block 1 block 2 block 3 block 4 block 5 block 6 block 7 block 8

“1” a b c d e f ab

bf be bd bc af ae ad ac

abc abd df ef cf de cd ce

abf abe ace acd ade acf aef adf

bef bdf cdf cde bce bcd bde bcf

def abce acdf def abdf abef abcf abde

adef acef bcdef acde bdef bcef bcdf bcde

abdef acdef cef abcdef cdef abcde abcdf abcef

————————————————————————————————

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)