HW11_EJL
Eddie J. Liu
11/14/2021
Problem 6.8
obs <- c(21,23,20,22,28,26,25,24,29,26,25,27,37,38,35,39,38,36,31,29,30,34,33,35)
time<- c(-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1,1,1,1,1)
cm <- c(-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1)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':
##
## lengthsdat1 <- data.frame(time,cm,obs)
mod <- lm(obs~time*cm,data = dat1)
coef(mod)## (Intercept) time cm time:cm
## 29.625000 4.958333 -0.625000 -1.958333halfnormal(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] time time:cm
dat1$time <- as.factor(dat1$time)
dat1$cm <- as.factor(dat1$cm)
model<-aov(obs~time+cm+time*cm,data=dat1)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## time 1 590.0 590.0 115.506 9.29e-10 ***
## cm 1 9.4 9.4 1.835 0.190617
## time:cm 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 ' ' 1plot(model)
\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\alpha\beta_{ij}+\epsilon_{ijk}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\). all i
\(H_{1}\): \(\beta_{j}\neq0\), some i
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i
P-value for time:cm, cm, time are 0.000397<0.05, 0.190617>0.05, 9.29e-10 <0.05, interaction is significant, so we do not need to investigate main effect. It confirms our halfnorm plot that time:cm is significant. We herein plot the interaction.
We do have MSE=5.1, which means we do not running out the degrees of freedom and had a normal f-test.
Model inadequacy: from residual plot and norm plot, we could roughly say they are normally distributed and variances are roughly equal.
##Analyze the residuals:
I <- c(21,25,37,31)
II <- c(22,26,39,34)
III <- c(23,24,38,29)
IV <- c(28,25,38,33)
V <- c(20,29,35,30)
VI <- c(26,27,36,35)
obs <- cbind(I,II,III,IV,V,VI)
time<- c(-1,-1,1,1)
cm<- c(-1,1,-1,1)
timecm <- time*cm
gmean <- mean(obs)
gmean## [1] 29.625mean <- mean(c(I,II,III,IV,V,VI))
mean## [1] 29.625dat <- as.data.frame(cbind(time,cm,timecm,obs))
dat## time cm timecm I II III IV V VI
## 1 -1 -1 1 21 22 23 28 20 26
## 2 -1 1 -1 25 26 24 25 29 27
## 3 1 -1 -1 37 39 38 38 35 36
## 4 1 1 1 31 34 29 33 30 35dat$sum <- c(sum(dat[1,4:9]),sum(dat[2,4:9]),sum(dat[3,4:9]),sum(dat[4,4:9]))
dat## time cm timecm I II III IV V VI sum
## 1 -1 -1 1 21 22 23 28 20 26 140
## 2 -1 1 -1 25 26 24 25 29 27 156
## 3 1 -1 -1 37 39 38 38 35 36 223
## 4 1 1 1 31 34 29 33 30 35 192effect_factor <- 2/(2^2*6)
effect_time <- (-dat[1,10]-dat[2,10]+dat[3,10]+dat[4,10])*effect_factor # effect of factor A
effect_cm <- (-dat[1,10]+dat[2,10]-dat[3,10]+dat[4,10])*effect_factor # effect of factor B
effect_timecm <- (dat[1,10]-dat[4,10]-dat[2,10]+dat[3,10])*effect_factor # effect of factor interaction AB
effect_factor## [1] 0.08333333effect_time## [1] 9.916667effect_cm## [1] -1.25effect_timecm## [1] 1.25effect_time= 9.916667, effect_cm=-1.25, effect_timecm=1.25, groundmean=29.625.
y1=29.625+9.916667/2*(-1)+(-1.25)/2*(-1)+1.25/2*(1)
y2=29.625+9.916667/2*(-1)+(-1.25)/2*(1)+1.25/2*(-1)
y3=29.625+9.916667/2*(1)+(-1.25)/2*(-1)+1.25/2*(-1)
y4=29.625+9.916667/2*(1)+(-1.25)/2*(1)+1.25/2*(1)
y1## [1] 25.91667y2## [1] 23.41667y3## [1] 34.58333y4## [1] 34.58333e1=I[1]-y1
e2=II[1]-y1
e3=II[1]-y1
e4=IV[1]-y1
e5=V[1]-y1
e6=VI[1]-y1
e7=I[2]-y2
e8=II[2]-y2
e9=III[2]-y2
e10=IV[2]-y2
e11=V[2]-y2
e12=VI[2]-y2
e13=I[3]-y3
e14=II[3]-y3
e15=III[3]-y3
e16=IV[3]-y3
e17=V[3]-y3
e18=VI[3]-y3
e19=I[4]-y4
e20=II[4]-y4
e21=III[4]-y4
e22=IV[4]-y4
e23=V[4]-y4
e24=VI[4]-y4
e1## [1] -4.916666e2## [1] -3.916666e3## [1] -3.916666e4## [1] 2.083334e5## [1] -5.916666e6## [1] 0.0833335e7## [1] 1.583334e8## [1] 2.583334e9## [1] 0.5833335e10## [1] 1.583334e11## [1] 5.583334e12## [1] 3.583334e13## [1] 2.416666e14## [1] 4.416666e15## [1] 3.416666e16## [1] 3.416666e17## [1] 0.4166665e18## [1] 1.416666e19## [1] -3.583334e20## [1] -0.5833335e21## [1] -5.583334e22## [1] -1.583334e23## [1] -4.583334e24## [1] 0.4166665residual values are liseted above.
————————————————————————————————
Problem 6.12
A <- c(-1,1,-1,1)
B <- c(-1,-1,1,1)
AB <- A*B
Factor <- cbind(A,B,AB)
Factor## A B AB
## [1,] -1 -1 1
## [2,] 1 -1 -1
## [3,] -1 1 -1
## [4,] 1 1 1I <- c(14.037,13.880,14.821,14.888)
II <- c(16.165,13.860,14.757,14.921)
III <- c(13.972,14.032,14.843,14.415)
IV <- c(13.907,13.914,14.878,14.932)
obs <- cbind(I,II,III,IV)
mean(c(I,II,III,IV))## [1] 14.51388(a):
dat <- as.data.frame(cbind(A,B,AB,obs))
dat$sum <- c(sum(dat[1,4:7]),sum(dat[2,4:7]),sum(dat[3,4:7]),sum(dat[4,4:7]))
dat## A B AB I II III IV sum
## 1 -1 -1 1 14.037 16.165 13.972 13.907 58.081
## 2 1 -1 -1 13.880 13.860 14.032 13.914 55.686
## 3 -1 1 -1 14.821 14.757 14.843 14.878 59.299
## 4 1 1 1 14.888 14.921 14.415 14.932 59.156effect_factor <- 2/(2^2*4)
effect_A <- ((-dat[1,8]-dat[3,8])+(dat[2,8]+dat[4,8]))*effect_factor # effect of factor A
effect_B <- ((-dat[1,8]-dat[2,8])+(dat[3,8]+dat[4,8]))*effect_factor # effect of factor B
effect_AB <- (dat[1,8]+dat[4,8]-dat[2,8]-dat[3,8])*effect_factor # effect of factor interaction AB
effect_factor## [1] 0.125effect_A## [1] -0.31725effect_B## [1] 0.586effect_AB## [1] 0.2815Estimate the factor effects from equation: \(effect=\frac{2(contrast)}{2^k*n}\)
Thus, effect_A=-0.31725, effect_B=0.586, effect_AB=0.2815, groundmean= 14.51388
library(DoE.base)
A <- rep(A,4)
B <- rep(B,4)
obs <- c(I,II,III,IV)
dat1 <- data.frame(A,B,obs)
dat1## A B obs
## 1 -1 -1 14.037
## 2 1 -1 13.880
## 3 -1 1 14.821
## 4 1 1 14.888
## 5 -1 -1 16.165
## 6 1 -1 13.860
## 7 -1 1 14.757
## 8 1 1 14.921
## 9 -1 -1 13.972
## 10 1 -1 14.032
## 11 -1 1 14.843
## 12 1 1 14.415
## 13 -1 -1 13.907
## 14 1 -1 13.914
## 15 -1 1 14.878
## 16 1 1 14.932mod <- lm(obs~A*B,data = dat1)
coef(mod)## (Intercept) A B A:B
## 14.513875 -0.158625 0.293000 0.140750halfnormal(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] e7 B A
dat1$A <- as.factor(dat1$A)
dat1$B <- as.factor(dat1$B)
model<-aov(obs~A+B+A*B,data=dat1)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 0.403 0.4026 1.262 0.2833
## B 1 1.374 1.3736 4.305 0.0602 .
## A:B 1 0.317 0.3170 0.994 0.3386
## Residuals 12 3.828 0.3190
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model)
Assuming critical value \(\alpha=0.05\).
\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\alpha\beta_{ij}+\epsilon_{ijk}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\). all i
\(H_{1}\): \(\beta_{j}\neq0\), some i
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i
P-value for A:B, A, B are 0.3386, 0.2833, 0.0602. Since p(interaction)>0.05, so it is insignificant, fail to reject null hypothesis. Then we do need to investigate main effect. Since p(A)>0.05 and p(B)>0.05, effect of A and effect of B are insignificant, both fail to reject null hypothesis. Repeatedly, it depends on the settings of the critical value.
(b):
Analysis of variance: Since I=2,J=2,K=4 and no interaction per se.
We can use these following equations:
\(\sigma_{\alpha}^2=\frac{MSA-MSE}{JK}=\frac{0.4026-0.3190}{2*4}=0.01045\).
\(\sigma_{\beta}^2=\frac{MSB-MSE}{IK}=\frac{1.3736-0.3190}{2*4}=0.131826\).
E[MSE]=\(\sigma^2=0.3190\).
(c):
y= groundmean+(effect_A/2)(coded factorA)+(effect_B/2)(coded factorB)+(effect_AB/2)(coded factorAB)
Since only B may affect results depending on alpha =0.05 or 0.1. equation becomes y= groundmean or y= groundmean+(effect_B/2)(coded factorB)
(d):
if assume critical value alpha=0.1,thus B does impact:
y1=14.51388+(0.586)/2*(-1)
y2=14.51388+(0.586)/2*(-1)
y3=14.51388+(0.586)/2*(1)
y4=14.51388+(0.586)/2*(1)
y1## [1] 14.22088y2## [1] 14.22088y3## [1] 14.80688y4## [1] 14.80688e1=I[1]-y1
e2=II[1]-y1
e3=II[1]-y1
e4=IV[1]-y1
e7=I[2]-y2
e8=II[2]-y2
e9=III[2]-y2
e10=IV[2]-y2
e13=I[3]-y3
e14=II[3]-y3
e15=III[3]-y3
e16=IV[3]-y3
e19=I[4]-y4
e20=II[4]-y4
e21=III[4]-y4
e22=IV[4]-y4
e1## [1] -0.18388e2## [1] 1.94412e3## [1] 1.94412e4## [1] -0.31388e7## [1] -0.34088e8## [1] -0.36088e9## [1] -0.18888e10## [1] -0.30688e13## [1] 0.01412e14## [1] -0.04988e15## [1] 0.03612e16## [1] 0.07112e19## [1] 0.08112e20## [1] 0.11412e21## [1] -0.39188e22## [1] 0.12512y1=14.51388+(0.586)/2(-1)
y2=14.51388+(0.586)/2(-1)
y3=14.51388+(0.586)/2(1)
y4=14.51388+(0.586)/2(1)
It has concern, if assume critical value alpha=0.05,thus no one impacts: y=14.51388, the resultant one should be e=variable-y. We will not calculate them all here.
(e) Discuss: It is the same we will not calculate them all. So it is hoghly depends on the critical value settings.
MSE = 0.3190, df(residuals) =12. Also from the residual plot, it is not roughly normally distributed and variances are not roughly equal.
————————————————————————————————
Problem 6.21
obs <- 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 <- c(rep(-1,7),rep(1,7),rep(-1,7),rep(1,7),rep(-1,7),rep(1,7),rep(-1,7),rep(1,7),rep(-1,7),rep(1,7),rep(-1,7),rep(1,7),rep(-1,7),rep(1,7),rep(-1,7), rep(1,7))
B <- c(rep(-1,7),rep(-1,7),rep(1,7),rep(1,7),rep(-1,7),rep(-1,7),rep(1,7),rep(1,7),rep(-1,7),rep(-1,7),rep(1,7),rep(1,7),rep(-1,7),rep(-1,7),rep(1,7),rep(1,7))
C <- c(rep(-1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(1,7),rep(1,7),rep(1,7),rep(1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(1,7),rep(1,7),rep(1,7),rep(1,7))
D <- c(rep(-1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(-1,7),rep(1,7),rep(1,7),rep(1,7),rep(1,7),rep(1,7),rep(1,7),rep(1,7),rep(1,7))
dat1 <- data.frame(A,B,C,D,obs)str(dat1)## 'data.frame': 112 obs. of 5 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 ...
## $ obs: num 10 18 14 12.5 19 16 18.5 0 16.5 4.5 ...library(DoE.base)
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
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 e95 e28 e44 e49 B e84 e32 e78
(a):
From halfnorm plot, A and B are significant.
dat1$A <- as.factor(dat1$A)
dat1$B <- as.factor(dat1$B)
dat1$C <- as.factor(dat1$C)
dat1$D <- as.factor(dat1$D)
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 ' ' 1plot(model)
(b):
\(\overline{y_{ijklm}}=\mu+\alpha_{i}+\beta_{j}+\gamma_{k}+\delta_{l}+\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_{ijklm}\)
We only retain A as well as B and exclude the others.
\(\overline{y_{iklm}}=\mu+\alpha_{i}+\beta_{k}+\epsilon_{iklm}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{k}=0\), all j
\(H_{1}\): \(\beta_{k}\neq0\), some j
P-value for A and B are 0.00136, 0.03469. They are all below 0.05. We reject null hypothesis for A and B. Their effects on means are significant, which confirms the halfnorm plot.
For residuals, Df=109,SSE=9249, MSE=84.9
Since A B are the onlt significant.
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 <- cbind(I,II,III,IV,V,VI,VII)
mean(c(I,II,III,IV,V,VI,VII))## [1] 12.29911effect and RESIDUAL analysis similar to questions above y= groundmean+(effect_A/2)(coded factorA)+(effect_B/2)(coded factorB)+(effect_AB/2)(coded factorAB), will not repeat bulky calculations.
model inadequacy: from the halfnorm plot and norm plot, we could say most of the data follows normal distribution, while only some of them violates. But overall, they could pass fat pencil test. So it is ok. From residual plot, also variances are roughly equal.
————————————————————————————————
Problem 6.36
obs <- 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)
A <- c(-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1)
B <- c(-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1)
C <- c(-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1)
D <- c(-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1)
dat1 <- data.frame(A,B,C,D,obs)str(dat1)## 'data.frame': 16 obs. of 5 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 ...
## $ obs: num 1.92 11.28 1.09 5.75 2.13 ...library(DoE.base)
mod <- lm(obs~A*B*C*D,data = dat1)
coef(mod)## (Intercept) A B C D A:B
## 4.680625 3.160625 -1.501875 -0.220625 -0.079375 -1.069375
## A:C B:C A:D B:D C:D A:B:C
## -0.298125 0.229375 -0.056875 -0.046875 0.029375 0.344375
## A:B:D A:C:D B:C:D A:B:C:D
## -0.096875 -0.010625 0.094375 0.141875## ?halfnormal
halfnormal(mod)##
## Significant effects (alpha=0.05, Lenth method):## [1] A B A:B A:B:C
(a):
From halfnorm plot, we have A, B, A:B, A:B:C, which is significant.
dat1$A <- as.factor(dat1$A)
dat1$B <- as.factor(dat1$B)
dat1$C <- as.factor(dat1$C)
dat1$D <- as.factor(dat1$D)
model<-aov(obs~A+B+C+A*B*C,data=dat1)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 159.83 159.83 1563.061 1.84e-10 ***
## B 1 36.09 36.09 352.937 6.66e-08 ***
## C 1 0.78 0.78 7.616 0.02468 *
## A:B 1 18.30 18.30 178.933 9.33e-07 ***
## A:C 1 1.42 1.42 13.907 0.00579 **
## B:C 1 0.84 0.84 8.232 0.02085 *
## A:B:C 1 1.90 1.90 18.556 0.00259 **
## Residuals 8 0.82 0.10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model)
\(\overline{y_{ijklm}}=\mu+\alpha_{i}+\beta_{j}+\gamma_{k}+\delta_{l}+\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_{ijklm}\)
(b):
For residuals, df=8, SSE=0.82, MSE=0.10
model inadequacy: from residual plot and norm plot, we could roughly say they are normally distributed, but variances are not roughly equal.
(d):
\(\overline{y_{iklm}}=\mu+\alpha_{i}+\beta_{j}\gamma_{k}+\alpha\beta_{ij}+\alpha\gamma_{ik}+\beta\gamma_{jk}+\alpha\beta\gamma_{ijk}+\epsilon_{ijkm}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{k}=0\), all j
\(H_{1}\): \(\beta_{k}\neq0\), some j
\(H_{0}\): \(\gamma_{k}=0\), all k
\(H_{1}\): \(\gamma_{k}\neq0\), some k
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i,j
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i,j
\(H_{0}\): \(\alpha\gamma_{ik}=0\), all i,k
\(H_{1}\): \(\alpha\gamma_{ik}\neq0\), some i,k
\(H_{0}\): \(\beta\gamma_{jk}=0\), all j,k
\(H_{1}\): \(\beta\gamma_{jk}\neq0\), some j,k
\(H_{0}\): \(\alpha\beta\gamma_{ijk}=0\), all j,k
\(H_{1}\): \(\alpha\beta\gamma_{ijk}\neq0\), some j,k
P(A:B:C)=0.00259<0.05, so we verify three order interaction is significant.
(c): Transformation
dat1 <- data.frame(A,B,C,D,log(obs))
str(dat1)## 'data.frame': 16 obs. of 5 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 ...
## $ log.obs.: num 0.6523 2.423 0.0862 1.7492 0.7561 ...library(DoE.base)
mod <- lm(log(obs)~A*B*C*D,data = dat1)
coef(mod)## (Intercept) A B C D A:B
## 1.185417116 0.812870345 -0.314277554 -0.006408558 -0.018077390 -0.024684570
## A:C B:C A:D B:D C:D A:B:C
## -0.039723700 -0.004225796 -0.009578245 0.003708723 0.017780432 0.063434408
## A:B:D A:C:D B:C:D A:B:C:D
## -0.029875960 -0.003740235 0.003765760 0.031322043## ?halfnormal
halfnormal(mod)##
## Significant effects (alpha=0.05, Lenth method):## [1] A B A:B:C
From halfnorm plot, we have A, B, A:B:C, which is significant.
dat1$A <- as.factor(dat1$A)
dat1$B <- as.factor(dat1$B)
dat1$C <- as.factor(dat1$C)
dat1$D <- as.factor(dat1$D)
model<-aov(log(obs)~A+B+C+A*B*C,data=dat1)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 10.572 10.572 1994.556 6.98e-11 ***
## B 1 1.580 1.580 298.147 1.29e-07 ***
## C 1 0.001 0.001 0.124 0.73386
## A:B 1 0.010 0.010 1.839 0.21207
## A:C 1 0.025 0.025 4.763 0.06063 .
## B:C 1 0.000 0.000 0.054 0.82223
## A:B:C 1 0.064 0.064 12.147 0.00826 **
## Residuals 8 0.042 0.005
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model)
linear equation remains the same and P(A:B:C)=0.00826<0.05, so we verify three order interaction is significant, we stop explorations on two order interactions and main effects.
Transformation has been useful. Since from residual plot, we can see that variances’ differneces have been narrowed down. Now roughly equal.
model adequacy: after transformation, from residual plot and norm plot, we could roughly say they are normally distributed and variances are roughly equal.
#(d): ### same liner model and same hypothesis as above.
effect and RESIDUAL analysis similar to questions above y= groundmean+(effect_A/2)(coded factorA)+(effect_B/2)(coded factorB)+(effect_AB/2)(coded factorAB), will not repeat bulky calculations.
————————————————————————————————
Problem 6.39
obs <- c(8.11,5.56,5.77,5.82,9.17,7.8,3.23,5.69,8.82,14.23,9.2,8.94,8.68,11.49,6.25,9.12,7.93,5,7.47,12,9.86,3.65,6.4,11.61,12.43,17.55,8.87,25.38,13.06,18.85,11.78,26.05)
A <- c(-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1)
B <- c(-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1)
C <- c(-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1)
D <- c(-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1)
E <- c(rep(-1,16),rep(1,16))
dat1 <- data.frame(A,B,C,D,E,obs)str(dat1)## 'data.frame': 32 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 ...
## $ E : num -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ obs: num 8.11 5.56 5.77 5.82 9.17 ...library(DoE.base)
mod <- lm(obs~A*B*C*D*E,data = dat1)
coef(mod)## (Intercept) A B C D E
## 10.1803125 1.6159375 0.0434375 -0.0121875 2.9884375 2.1878125
## A:B A:C B:C A:D B:D C:D
## 1.2365625 -0.0015625 -0.1953125 1.6665625 -0.0134375 0.0034375
## A:E B:E C:E D:E A:B:C A:B:D
## 1.0271875 1.2834375 0.3015625 1.3896875 0.2503125 -0.3453125
## A:C:D B:C:D A:B:E A:C:E B:C:E A:D:E
## -0.0634375 0.3053125 1.1853125 -0.2590625 0.1709375 0.9015625
## B:D:E C:D:E A:B:C:D A:B:C:E A:B:D:E A:C:D:E
## -0.0396875 0.3959375 -0.0740625 -0.1846875 0.4071875 0.1278125
## B:C:D:E A:B:C:D:E
## -0.0746875 -0.3553125## ?halfnormal
halfnormal(mod)##
## Significant effects (alpha=0.05, Lenth method):## [1] D E A:D A D:E B:E A:B A:B:E A:E A:D:E
# (a): ### From halfnorm plot, we have A, E, D, A:B, A:D, A:E, B:E, D:E, A:B:E A:D:E, which is significant.
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$E <- as.factor(dat1$E)
model<-aov(obs~A+B+D+E+A*B+A*D+A*E+B*E+D*E+A*B*E+A*D*E,data=dat1)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 83.56 83.56 51.362 6.10e-07 ***
## B 1 0.06 0.06 0.037 0.849178
## D 1 285.78 285.78 175.664 2.30e-11 ***
## E 1 153.17 153.17 94.149 5.24e-09 ***
## A:B 1 48.93 48.93 30.076 2.28e-05 ***
## A:D 1 88.88 88.88 54.631 3.87e-07 ***
## A:E 1 33.76 33.76 20.754 0.000192 ***
## B:E 1 52.71 52.71 32.400 1.43e-05 ***
## D:E 1 61.80 61.80 37.986 5.07e-06 ***
## A:B:E 1 44.96 44.96 27.635 3.82e-05 ***
## A:D:E 1 26.01 26.01 15.988 0.000706 ***
## Residuals 20 32.54 1.63
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model)
\(\overline{y_{ijlmn}}=\mu+\alpha_{i}+\beta_{j}+\delta_{l}+\zeta_{m}+\alpha\beta_{ij}+\alpha\delta_{il}+\alpha\zeta_{im}+\beta\zeta_{jm}+\delta\zeta_{lm}+\alpha\beta\zeta_{ijm}+\alpha\delta\zeta_{ilm}+\epsilon_{ijlmn}\)
\(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
\(H_{0}\): \(\delta_{l}=0\), all l
\(H_{1}\): \(\delta_{l}\neq0\), some l
\(H_{0}\): \(\zeta_{m}=0\), all m
\(H_{1}\): \(\zeta_{m}\neq0\), some m
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i,j
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i,j
\(H_{0}\): \(\alpha\delta_{il}=0\), all i,l
\(H_{1}\): \(\alpha\delta_{il}\neq0\), some i,l
\(H_{0}\): \(\alpha\zeta_{im}=0\), all i,m
\(H_{1}\): \(\alpha\zeta_{im}\neq0\), some j,m
\(H_{0}\): \(\beta\zeta_{jm}=0\), all j,m
\(H_{1}\): \(\beta\zeta_{jm}\neq0\), some j,m
\(H_{0}\): \(\delta\zeta_{lm}=0\), all l,m
\(H_{1}\): \(\delta\zeta_{lm}\neq0\), some l,m
\(H_{0}\): \(\alpha\beta\zeta_{ijm}=0\), all i,j,m
\(H_{1}\): \(\alpha\beta\zeta_{ijm}\neq0\), some i,j,m
\(H_{0}\): \(\alpha\delta\zeta_{ilm}=0\), all i,l,m
\(H_{1}\): \(\alpha\delta\zeta_{ilm}\neq0\), some i,l,m
P(A:B:E)=3.82e-05<0.05 and P(A:D:E)=0.000706<0.05, so we verify three order interactions of these two are significant. We stop explorations on two order interactions and main effects.
(b): We do have MSE=1.63, which means we do not running out the degrees of freedom and had a normal f-test. model inadequacy: from residual plot and norm plot, we could roughly say they are normally distributed and variances are roughly equal.
(c): Since p(B)=0.849178>0.05, B factor is not significant. If we drop this factor, We only retain A, E, D, A:D, A:E, D:E, A:D:E and their related components and exclude the others.
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$E <- as.factor(dat1$E)
model<-aov(obs~A+D+E+A*D+A*E+D*E+A*D*E,data=dat1)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 83.56 83.56 11.191 0.002696 **
## D 1 285.78 285.78 38.275 2.16e-06 ***
## E 1 153.17 153.17 20.514 0.000137 ***
## A:D 1 88.88 88.88 11.903 0.002084 **
## A:E 1 33.76 33.76 4.522 0.043943 *
## D:E 1 61.80 61.80 8.277 0.008297 **
## A:D:E 1 26.01 26.01 3.484 0.074249 .
## Residuals 24 179.20 7.47
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model)
