Homework week 11
Dowthyaksai Lagadapati
2022-11-14
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(GAD)## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following object is masked from 'package:dplyr':
##
## count## Loading required package: R.methodsS3## R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.library(DoE.base)## Warning: package 'DoE.base' was built under R version 4.2.2## 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':
##
## lengths6.8)
cm <- rep(c(-1,-1,1,1),6)
time <- c(rep(-1,12),rep(1,12))
growth <- c(21,22,25,26,23,28,24,25,20,26,29,27,37,39,31,34,38,38,29,33,35,36,30,35)Model Equation :
\(y_{ijk}\) = \(\alpha_{i}\) + \(\beta_{j}\) + \(\alpha\beta_{ij}\) + \(\epsilon_{ijk}\)
Hypothesis :
Medium Effect :
\(H_{0}\) : \(\alpha_{i}\) = 0 \(\forall\) i ,
\(H_{a}\) : \(\alpha_{i}\) \(\neq\) 0 \(\exists\) i ,
Time Effect :
\(H_{0}\) : \(\beta_{j}\) = 0 \(\forall\) j ,
\(H_{a}\) : \(\beta_{j}\) \(\neq\) 0 \(\exists\) j ,
Interaction Effect :
\(H_{0}\) : \(\alpha\beta_{ij}\) = 0 \(\forall\) i, j ,
\(H_{a}\) : \(\alpha\beta_{ij}\) \(\neq\) \(\exists\) i, j ,
Manipulating Data and Running AOV
cm <- as.fixed(cm)
time <- as.fixed(time)
dat6.8 <- data.frame(growth,cm,time)
model6.8 <- aov(growth ~ time + cm + time*cm)
GAD::gad(model6.8)## Analysis of Variance Table
##
## Response: growth
## Df Sum Sq Mean Sq F value Pr(>F)
## time 1 590.04 590.04 115.5057 9.291e-10 ***
## cm 1 9.38 9.38 1.8352 0.1906172
## time:cm 1 92.04 92.04 18.0179 0.0003969 ***
## Residual 20 102.17 5.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From the variance analysis table we can observe that the p value of the interaction is less than the threshold \(\alpha\) = 0.05. So, we reject the null hypothesis and the interaction is significant.
interaction.plot(time,cm,growth)
plot(model6.8)
From the plots we can observe that the interaction is normally distributed and the variance also looks similar. So, we can conclude that the model is adequate.
6.12)
library(DoE.base)
af <- rep(c(-1,1,-1,1),4)
dt <- rep(c(-1,-1,1,1),4)
thickness <- c(14.037,13.880,14.821,14.888,16.165,13.860,14.757,14.921,13.972,14.032,14.843,14.415,13.907,13.914,14.878,14.932)
af <- as.factor(af)
dt <- as.factor(dt)
data6.12 <- data.frame(af,dt,thickness)a)
One <- c(14.037,16.165,13.972,13.907)
af <- c(13.88,13.86,14.032,13.914)
dt <- c(14.821,14.757,14.843,14.878)
inter <- c(14.888,14.921,14.415,14.932)
S1 <- sum(One)
SA <- sum(af)
SB <- sum(dt)
SAB <- sum(inter)
effectaf <- (2*(SA+SAB-S1-SB)/(4*4))
effectdt <- (2*(SB+SAB-S1-SA)/(4*4))
effectinter <- (2*(SA+SB-S1-SAB)/(4*4))
effectaf## [1] -0.31725effectdt## [1] 0.586effectinter## [1] -0.2815b)
model6.12 <- aov(thickness ~ af*dt, data = data6.12)
summary(model6.12)## Df Sum Sq Mean Sq F value Pr(>F)
## af 1 0.403 0.4026 1.262 0.2833
## dt 1 1.374 1.3736 4.305 0.0602 .
## af:dt 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 ' ' 1From the results we can conclude that the interaction between factors is insignificant. So, we remove the interaction and test the main effects.
model6.12 <- aov(thickness ~ af + dt, data = data6.12)
summary(model6.12)## Df Sum Sq Mean Sq F value Pr(>F)
## af 1 0.403 0.4026 1.263 0.2815
## dt 1 1.374 1.3736 4.308 0.0584 .
## Residuals 13 4.145 0.3189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From the anova analysis, we can conclude that the main effects are significant.
c)
model <- lm(thickness ~ af*dt, data = data6.12)
coef(model6.12)## (Intercept) af1 dt1
## 14.37950 -0.31725 0.58600summary(model6.12)## Df Sum Sq Mean Sq F value Pr(>F)
## af 1 0.403 0.4026 1.263 0.2815
## dt 1 1.374 1.3736 4.308 0.0584 .
## Residuals 13 4.145 0.3189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1d)
plot(model6.12)
By plotting the residuals, we can see the that the model is inadequate.
e)
We can perform a boxcox transformation on the data and then perform anova analysis on transformed data.
6.21)
FA <- rep(rep(c(-1,1),8),7)
FB <- rep(rep(c(1,1,-1,-1),4),7)
FC <- rep(rep(c(rep(1,4),rep(-1,4)),2),7)
FD <- rep(c(rep(1,8),rep(-1,8)),7)
R1 <- 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)
R2 <- 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)
R3 <- 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)
R4 <- 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)
R5 <- 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)
R6 <- 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)
R7 <- 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)
Obs6.21 <- c(R1,R2,R3,R4,R5,R6,R7)a)
dat6.21 <- data.frame(FA,FB,FC,FD,Obs6.21)
model6.21 <- lm(Obs6.21 ~ FA*FB*FC*FD,data=dat6.21)
halfnormal(model6.21)## Warning in halfnormal.lm(model6.21): halfnormal not recommended for models with
## more residual df than model df##
## Significant effects (alpha=0.05, Lenth method):## [1] FA e74 e92 FB e53
From the above plot clearly factors A and B are significant.
b)
model6.21 <- aov(Obs6.21 ~ FA*FB,data=dat6.21)
plot(model6.21)
## hat values (leverages) are all = 0.03571429
## and there are no factor predictors; no plot no. 5
6.36)
FacA <- c(-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1)
FacB <- c(-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1)
FacC <- c(-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1)
FacD <- c(-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1)
resistivity <- 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)
dat6.36 <- data.frame(FacA,FacB,FacC,FacD,resistivity)a)
model6.36 <- lm(resistivity ~ FacA*FacB*FacC*FacD,data = dat6.36)
coef(model6.36)## (Intercept) FacA FacB FacC
## 4.680625 3.160625 -1.501875 -0.220625
## FacD FacA:FacB FacA:FacC FacB:FacC
## -0.079375 -1.069375 -0.298125 0.229375
## FacA:FacD FacB:FacD FacC:FacD FacA:FacB:FacC
## -0.056875 -0.046875 0.029375 0.344375
## FacA:FacB:FacD FacA:FacC:FacD FacB:FacC:FacD FacA:FacB:FacC:FacD
## -0.096875 -0.010625 0.094375 0.141875halfnormal(model6.36)##
## Significant effects (alpha=0.05, Lenth method):## [1] FacA FacB FacA:FacB FacA:FacB:FacC
summary(model6.36)##
## Call:
## lm.default(formula = resistivity ~ FacA * FacB * FacC * FacD,
## data = dat6.36)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.68062 NaN NaN NaN
## FacA 3.16062 NaN NaN NaN
## FacB -1.50187 NaN NaN NaN
## FacC -0.22062 NaN NaN NaN
## FacD -0.07937 NaN NaN NaN
## FacA:FacB -1.06938 NaN NaN NaN
## FacA:FacC -0.29812 NaN NaN NaN
## FacB:FacC 0.22937 NaN NaN NaN
## FacA:FacD -0.05687 NaN NaN NaN
## FacB:FacD -0.04688 NaN NaN NaN
## FacC:FacD 0.02937 NaN NaN NaN
## FacA:FacB:FacC 0.34437 NaN NaN NaN
## FacA:FacB:FacD -0.09688 NaN NaN NaN
## FacA:FacC:FacD -0.01063 NaN NaN NaN
## FacB:FacC:FacD 0.09438 NaN NaN NaN
## FacA:FacB:FacC:FacD 0.14188 NaN NaN NaN
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NAFrom the half normal plot we can find that factors A, B, AB are significant. we will consider the tenative model.
\(Y_{ijk}\) = 4.680625 + 3.160625\(\alpha_{i}\) - 1.501875\(\beta_{j}\) - 1.069375\(\alpha\beta_{ij}\) + \(\epsilon_{ijkl}\)
b)
model6.36.2 <- aov(resistivity ~ FacA + FacB + FacC + FacA*FacB + FacA*FacB*FacC,data = dat6.36)
summary(model6.36.2)## Df Sum Sq Mean Sq F value Pr(>F)
## FacA 1 159.83 159.83 1563.061 1.84e-10 ***
## FacB 1 36.09 36.09 352.937 6.66e-08 ***
## FacC 1 0.78 0.78 7.616 0.02468 *
## FacA:FacB 1 18.30 18.30 178.933 9.33e-07 ***
## FacA:FacC 1 1.42 1.42 13.907 0.00579 **
## FacB:FacC 1 0.84 0.84 8.232 0.02085 *
## FacA:FacB:FacC 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(model6.36.2)
## hat values (leverages) are all = 0.5
## and there are no factor predictors; no plot no. 5
By observing the plots above we can observe that the model is inadequate.
The anova done later shows that all factors and interactions are significant at \(\alpha\) = 0.05.
c)
logresistivity <- log(resistivity)
dat6.36log <- data.frame(FacA,FacB,FacC,FacD,logresistivity)
model6.36log <- lm(logresistivity ~ FacA*FacB*FacC*FacD,data = dat6.36log)
coef(model6.36log)## (Intercept) FacA FacB FacC
## 1.185417116 0.812870345 -0.314277554 -0.006408558
## FacD FacA:FacB FacA:FacC FacB:FacC
## -0.018077390 -0.024684570 -0.039723700 -0.004225796
## FacA:FacD FacB:FacD FacC:FacD FacA:FacB:FacC
## -0.009578245 0.003708723 0.017780432 0.063434408
## FacA:FacB:FacD FacA:FacC:FacD FacB:FacC:FacD FacA:FacB:FacC:FacD
## -0.029875960 -0.003740235 0.003765760 0.031322043halfnormal(model6.36log)##
## Significant effects (alpha=0.05, Lenth method):## [1] FacA FacB FacA:FacB:FacC
summary(model6.36log)##
## Call:
## lm.default(formula = logresistivity ~ FacA * FacB * FacC * FacD,
## data = dat6.36log)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.185417 NaN NaN NaN
## FacA 0.812870 NaN NaN NaN
## FacB -0.314278 NaN NaN NaN
## FacC -0.006409 NaN NaN NaN
## FacD -0.018077 NaN NaN NaN
## FacA:FacB -0.024685 NaN NaN NaN
## FacA:FacC -0.039724 NaN NaN NaN
## FacB:FacC -0.004226 NaN NaN NaN
## FacA:FacD -0.009578 NaN NaN NaN
## FacB:FacD 0.003709 NaN NaN NaN
## FacC:FacD 0.017780 NaN NaN NaN
## FacA:FacB:FacC 0.063434 NaN NaN NaN
## FacA:FacB:FacD -0.029876 NaN NaN NaN
## FacA:FacC:FacD -0.003740 NaN NaN NaN
## FacB:FacC:FacD 0.003766 NaN NaN NaN
## FacA:FacB:FacC:FacD 0.031322 NaN NaN NaN
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NAAfter the log transformation we can observe that the half normal plot only characterizes factors A, B and ABC as significant. ANOVA analysis also presents factors A, B and ABC as significant.
model6.36log1 <- aov(logresistivity ~ FacA + FacB + FacC + FacA*FacB*FacC,data = dat6.36log)
summary(model6.36log1)## Df Sum Sq Mean Sq F value Pr(>F)
## FacA 1 10.572 10.572 1994.556 6.98e-11 ***
## FacB 1 1.580 1.580 298.147 1.29e-07 ***
## FacC 1 0.001 0.001 0.124 0.73386
## FacA:FacB 1 0.010 0.010 1.839 0.21207
## FacA:FacC 1 0.025 0.025 4.763 0.06063 .
## FacB:FacC 1 0.000 0.000 0.054 0.82223
## FacA:FacB:FacC 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(model6.36log1)
## hat values (leverages) are all = 0.5
## and there are no factor predictors; no plot no. 5
After the log transformation, we can observe that data is almost falling in a straight line .It hasn't improved considerable but variance has been stabilized to some extent and thus we can conclude that transformation has been useful.
d)
\(Y_{ijkl}\) = 1.185417 + 0.812870345\(\alpha_{i}\) - 0.314277554\(\beta_{j}\) + \(\epsilon_{ijkl}\)
6.39)
A <- rep(c(-1,1),16)
B <- rep(c(-1,-1,1,1),8)
C <- rep(c(-1,-1,-1,-1,1,1,1,1),4)
D <- rep(c(-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1),2)
E <- c(rep(-1,16),rep(1,16))
y <- 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)
dat6.39 <- data.frame(A,B,C,D,E,y)a)
model6.39 <- lm(y ~ A*B*C*D*E,data = dat6.39)
halfnormal(model6.39)##
## 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
coef(model6.39)## (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.3553125dat6.39 <- data.frame(A,B,D,E,y)
model6.39.2 <- aov(y ~ A + B + D + E + A*B + B*E + D*E + A*D + A*E + A*B*E + A*D*E,data = dat6.39)From half normal plot the factors A, D, E, A:D, D:E, B:E, A:B, A:E, A:B:E, A:D:E as significant. ANOVA analysis also presents the factors A, D, E, AB, AD, AE, BE, DE, ABE, ADE as significant.
b)
plot(model6.39.2)
## hat values (leverages) are all = 0.375
## and there are no factor predictors; no plot no. 5
From the Normal and Residual plot we can conclude that the model is inadequate. The data though satisfies the normality it fails as variance is spread wide.
c)
model6.39.3 <- aov(y ~ A*B*D*E,data = dat6.39)
halfnormal(model6.39.3)##
## 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 e10
coef(model6.39.3)## (Intercept) A B D E A:B
## 10.1803125 1.6159375 0.0434375 2.9884375 2.1878125 1.2365625
## A:D B:D A:E B:E D:E A:B:D
## 1.6665625 -0.0134375 1.0271875 1.2834375 1.3896875 -0.3453125
## A:B:E A:D:E B:D:E A:B:D:E
## 1.1853125 0.9015625 -0.0396875 0.4071875summary(model6.39.3)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 83.56 83.56 57.233 1.14e-06 ***
## B 1 0.06 0.06 0.041 0.841418
## D 1 285.78 285.78 195.742 2.16e-10 ***
## E 1 153.17 153.17 104.910 1.97e-08 ***
## A:B 1 48.93 48.93 33.514 2.77e-05 ***
## A:D 1 88.88 88.88 60.875 7.66e-07 ***
## B:D 1 0.01 0.01 0.004 0.950618
## A:E 1 33.76 33.76 23.126 0.000193 ***
## B:E 1 52.71 52.71 36.103 1.82e-05 ***
## D:E 1 61.80 61.80 42.328 7.24e-06 ***
## A:B:D 1 3.82 3.82 2.613 0.125501
## A:B:E 1 44.96 44.96 30.794 4.40e-05 ***
## A:D:E 1 26.01 26.01 17.815 0.000650 ***
## B:D:E 1 0.05 0.05 0.035 0.854935
## A:B:D:E 1 5.31 5.31 3.634 0.074735 .
## Residuals 16 23.36 1.46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model6.39log5 <- lm(y ~ A + B + D + E, data = dat6.39)
coef(model6.39log5)## (Intercept) A B D E
## 10.1803125 1.6159375 0.0434375 2.9884375 2.1878125Factor C was insignificant hence it is droped and did analysis again but got similar results in part a.
d)
\(Y_{ijkl}\) = 10.1803125 + 1.6159375\(\alpha_{i}\) + 0.0434375\(\beta_{j}\) + 2.9884375\(\gamma_{k}\) + 2.1878125\(\delta_{l}\) + \(\epsilon_{ijkl}\)
From the above model equation we can maximize the predicted response.