5.4

(a)

library(GAD)
## Loading required package: matrixStats
## Loading required package: R.methodsS3
## R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut  <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
cut <- as.fixed(cut)
FR <- as.fixed(FR)
model<-aov(finish~cut+FR+cut*FR)
GAD::gad(model)
## Analysis of Variance Table
## 
## Response: finish
##          Df  Sum Sq Mean Sq F value    Pr(>F)    
## cut       3 2125.11  708.37 24.6628 1.652e-07 ***
## FR        2 3160.50 1580.25 55.0184 1.086e-09 ***
## cut:FR    6  557.06   92.84  3.2324   0.01797 *  
## Residual 24  689.33   28.72                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Since P-value of Cut:FR=0.01797<0.05, reject null hypothesis, significant interaction between Cut and Feed Rate. So herein we stop exploration on the main factors, even if we got P-value of Cut=1.652e-07<0.05 and P-value of Feed Rate= 1.086e-09<0.05, indicating both significant impacts on means. We need to plot the interaction.

Interaction plot:

f1 <- finish[1:3]
f2 <- finish[4:6]
f3 <- finish[7:9]
f4 <- finish[10:12]
f5 <- finish[13:15]
f6 <- finish[16:18]
f7 <- finish[19:21]
f8 <- finish[22:24]
f9 <- finish[25:27]
f10 <- finish[28:30]
f11 <- finish[31:33]
f12 <- finish[34:36]
FR1 <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
avg1 <- ave(f1)
avg2 <- ave(f2)
avg3 <- ave(f3)
avg4 <- ave(f4)
avg5 <- ave(f5)
avg6 <- ave(f6)
avg7 <- ave(f7)
avg8 <- ave(f8)
avg9 <- ave(f9)
avg10 <- ave(f10)
avg11 <- ave(f11)
avg12 <- ave(f12)
?plo
## No documentation for 'plo' in specified packages and libraries:
## you could try '??plo'
plot(FR1, c(f1,f2,f3,f4),main="Scatterplot with average", xlab="Feed Rate",ylab="Finish", pch=19, col="red", ylim=range(c(min(finish), max(finish))))
lines(FR1, c(avg1,avg2,avg3,avg4), type = "b", col = "red", lwd = 3, pch = 15)

points(FR1, c(f5,f6,f7,f8),pch = 19,col="blue")
lines(FR1, c(avg5,avg6,avg7,avg8), type = "b", col = "blue", lwd = 3, pch = 15)

points(FR1, c(f9,f10,f11,f12),pch = 19,col="green")
lines(FR1, c(avg9,avg10,avg11,avg12), type = "b", col = "green", lwd = 3, pch = 15)

legend(x=0.21,y=80,c("Feed rate=0.2 in/min","Feed rate=0.25 in/min","Feed rate=0.30 in/min"),pch=c(19,19),col = c("red","blue","green"))

From interaction plot, we verify that do have significant interaction.

(b)

plot (model)

Model’s adequacy: From the residual plot and normality plot, we draw that since it passes the fat pencil test, it is roughly normally distributed. From the residual plot, variances are roughly equal.

(c)

finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))

f1 <- finish[1:12]
f2 <- finish[13:24]
f3 <- finish[25:36]

sdev <-c(rep(sd(f1),12),rep(sd(f2),12),rep(sd(f3),12)) 
avg <- c(ave(f1), ave(f2), ave(f3))

plot(FR, finish, ylim=range(c(avg-sdev, avg+sdev)), main="Scatterplot with average and error bar", xlab="Feed Rate ", ylab="Surface Finish", pch=19)

points(FR,avg,pch = 15,col="red")
arrows(FR, avg-sdev, FR, avg+sdev, length=0.05, angle=90, code=3, col="red")
lines(FR, avg, type = "b", col = "red", lwd = 3, pch = 15)

##?lines()
avg
##  [1]  81.58333  81.58333  81.58333  81.58333  81.58333  81.58333  81.58333
##  [8]  81.58333  81.58333  81.58333  81.58333  81.58333  97.58333  97.58333
## [15]  97.58333  97.58333  97.58333  97.58333  97.58333  97.58333  97.58333
## [22]  97.58333  97.58333  97.58333 103.83333 103.83333 103.83333 103.83333
## [29] 103.83333 103.83333 103.83333 103.83333 103.83333 103.83333 103.83333
## [36] 103.83333

For feed rate 0.20, average is 81.58333; For feed rate 0.25, average is 97.58333; For feed rate 0.30, average is 103.83333.

(d)

Please see part (a)

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

5.34

All fixed effect, we can use ANOVA or GAD

finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut  <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
block1 <- c(rep(seq(1,3),12))


dat1 <- cbind(finish,cut,FR,block1)
dat1 <- as.data.frame(dat1)

dat1$cut <- as.factor(dat1$cut)
dat1$FR <- as.factor(dat1$FR)
dat1$block1 <- as.factor(dat1$block1)
str(dat1)
## 'data.frame':    36 obs. of  4 variables:
##  $ finish: num  74 64 60 79 68 73 82 88 92 99 ...
##  $ cut   : Factor w/ 4 levels "0.15","0.18",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ FR    : Factor w/ 3 levels "0.2","0.25","0.3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ block1: Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
model<-aov(finish~cut+FR+block1+cut*FR,data=dat1)
summary(model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## cut          3 2125.1   708.4  30.637 4.89e-08 ***
## FR           2 3160.5  1580.2  68.346 3.64e-10 ***
## block1       2  180.7    90.3   3.907  0.03532 *  
## cut:FR       6  557.1    92.8   4.015  0.00726 ** 
## Residuals   22  508.7    23.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\gamma_{k}+\alpha\beta_{ij}+\epsilon_{ijkl}\)

\(\gamma_{k}\) is block effect.

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

Here is an important concept: block do not interact with other factors, so we have the linear equation and hypothesis above.

Since P-value of Cut: FR=0.00726<0.05, reject null hypothesis, significant interaction between Cut and Feed Rate. We stop exploration of main effect, even if P-value of Cut=4.89e-08<0.05 and P-value of Feed Rate= 3.64e-10<0.05, indicating significant impact of both on means. We need to plot interaction.

Interaction plot is the same as I plotted in Question 5.4.

Since we introduce the block effect, nuisance variability will introduce the sum square block (SSB) and decrease the sum square error(SSE). MSE is likely to drop, f increases and p-value decreases, all the p-value for Cut:FR, Cut and FR have been reduced. Blocking was useful in this experiment.

variance component: Since I=4,J=3,K=3,L=1

We can use these following equations:

\(\sigma_{\alpha\beta}^2=\frac{MSAB-MSE}{KL}=\frac{92.8-23.1}{3*1}=23.2333\)

\(\sigma_{\alpha}^2=\frac{MSA-MSAB}{JKL}=\frac{708.4-92.8}{3*3X1}=68.4\).

\(\sigma_{\beta}^2=\frac{MSB-MSAB}{IKL}=\frac{1580.2-92.8}{4*3X1}=123.95\).

\(\sigma_{\gamma}^2=\frac{MSC-MSE}{IJL}=\frac{90.3-23.1}{4*3X1}=5.6\).

E[MSE]=\(\sigma^2=23.1\).

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

Problem 13.5

Question assumes no interaction. If we count interaction in:

density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position+temp*position)
GAD::gad(model)
## Analysis of Variance Table
## 
## Response: density
##               Df Sum Sq Mean Sq  F value    Pr(>F)    
## temp           2 945342  472671 1155.518 0.0008647 ***
## position       1   7160    7160   15.998 0.0017624 ** 
## temp:position  2    818     409    0.914 0.4271101    
## Residual      12   5371     448                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (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

Question assumes no interaction. Since P-value of temp:position=0.4271101>0.05, we confirm no interaction.

P-value of Temp=0.0008647<0.05, P-value of Position=0.0017624<0.05, we reject both null hypothesis, both of them significantly affect the means.

Model’s adequacy: From norm plot and residual plot, we know that it is normally distributed and variances are roughly equal.

Analysis of variance: Since I=3,J=2,K=3, since no interaction,\(\sigma_{\alpha\beta}^2=0\)

We can use these following equations:

\(\sigma_{\alpha\beta}^2=\frac{MSAB-MSE}{K}=\frac{409-448}{3}<0\), so \(\sigma_{\alpha\beta}^2=0\).

\(\sigma_{\alpha}^2=\frac{MSA-MSAB}{JK}=\frac{472671-409}{2*3}=78710.3333\).

\(\sigma_{\beta}^2=\frac{MSB-MSAB}{IK}=\frac{7160-409}{3*3}=750.1111\).

E[MSE]=\(\sigma^2=448\).

Question assumes no interaction.If we do not count interaction in:

density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position)
GAD::gad(model)
## Analysis of Variance Table
## 
## Response: density
##          Df Sum Sq Mean Sq  F value    Pr(>F)    
## temp      2 945342  472671 1069.257 4.924e-16 ***
## position  1   7160    7160   16.197  0.001254 ** 
## Residual 14   6189     442                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (model)

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

\(H_{1}\): \(\beta_{j}\neq0\), some i

Question assumes no interaction.

P-value of Temp=4.924e-16<0.05, P-value of Position=0.001254<0.05, we reject both null hypothesis, both of them significantly affect the means.

Model’s adequacy: From norm plot and residual plot, we know that it is normally distributed and variances are roughly equal.

Analysis of variance: Since I=3,J=2,K=3 and no interaction per se.

We can use these following equations:

\(\sigma_{\alpha}^2=\frac{MSA-MSE}{JK}=\frac{472671-442}{2*3}=78704.8333\).

\(\sigma_{\beta}^2=\frac{MSB-MSE}{IK}=\frac{7160-442}{3*3}=746.4444\).

E[MSE]=\(\sigma^2=442\).

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

Problem 13.6

dat1 <- c(50,49,50,50,48,51,52,52,51,51,51,51,53,50,50,54,52,51,49,51,50,48,50,51,48,49,48,48,49,48,52,50,50,52,50,50,51,51,51,51,50,50,52,50,49,53,48,50,50,51,50,51,48,49,47,46,49,46,47,48)
operator <- c(rep(c(rep(1,3),rep(2,3)),10))
operator
##  [1] 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1
## [39] 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2
part <- c(rep(1,6),rep(2,6),rep(3,6),rep(4,6),rep(5,6),rep(6,6),rep(7,6),rep(8,6),rep(9,6),rep(10,6))
part
##  [1]  1  1  1  1  1  1  2  2  2  2  2  2  3  3  3  3  3  3  4  4  4  4  4  4  5
## [26]  5  5  5  5  5  6  6  6  6  6  6  7  7  7  7  7  7  8  8  8  8  8  8  9  9
## [51]  9  9  9  9 10 10 10 10 10 10
operator <- as.fixed(operator)
part <- as.random(part)
model<-aov(dat1~operator+part+operator*part)
GAD::gad(model)
## Analysis of Variance Table
## 
## Response: dat1
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## operator       1  0.417  0.4167  0.6923    0.4269    
## part           9 99.017 11.0019  7.3346 3.216e-06 ***
## operator:part  9  5.417  0.6019  0.4012    0.9270    
## Residual      40 60.000  1.5000                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (model)

Since test were made in random order indicated by question, there is not block effect, herein the measurements 1 2 3 do not mean blocking. This a mixed effect.

\(\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 operator:part=0.9270>0.05, fail to reject null hypothesis, no significant interaction. So analyze main effect. P-value for operator=0.4269>0.05, no impact by operator on means. P-value for part=3.216e-06<0.05, there is impact by part on means.

Model’s adequacy: From norm plot and residual plot, we can draw that it is roughly noramlly distributed, and variances are roughly equal. Assumptions are verified.

Analysis of variance: Since I=2,J=10,K=3, E[MSE]=\(\sigma^2\)=1.5000, E[MSAB]=\(\sigma^2+K\sigma_{\alpha\beta}^2\), \(0.6019=0.1500+3*\sigma_{\alpha\beta}^2\), so \(\sigma_{\alpha\beta}^2\)<0, assume \(\sigma_{\alpha\beta}^2\)=0.

We can ALSO use these following equations:

\(\sigma_{\alpha\beta}^2=\frac{MSAB-MSE}{K}=\frac{0.6019-1.5000}{3}<0\), so \(\sigma_{\alpha\beta}^2=0\).

\(\sigma_{\alpha}^2=\frac{MSA-MSAB}{JK}=\frac{0.4167-0.6019}{10*3}<0\), so \(\sigma_{\alpha}^2\)=0.

\(\sigma_{\beta}^2=\frac{MSB-MSAB}{IK}=\frac{11.0019-0.6019}{2*3}=1.7333\),

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

Raw Code:

library(GAD)
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut  <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
cut <- as.fixed(cut)
FR <- as.fixed(FR)
model<-aov(finish~cut+FR+cut*FR)
GAD::gad(model)
f1 <- finish[1:3]
f2 <- finish[4:6]
f3 <- finish[7:9]
f4 <- finish[10:12]
f5 <- finish[13:15]
f6 <- finish[16:18]
f7 <- finish[19:21]
f8 <- finish[22:24]
f9 <- finish[25:27]
f10 <- finish[28:30]
f11 <- finish[31:33]
f12 <- finish[34:36]
FR1 <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
avg1 <- ave(f1)
avg2 <- ave(f2)
avg3 <- ave(f3)
avg4 <- ave(f4)
avg5 <- ave(f5)
avg6 <- ave(f6)
avg7 <- ave(f7)
avg8 <- ave(f8)
avg9 <- ave(f9)
avg10 <- ave(f10)
avg11 <- ave(f11)
avg12 <- ave(f12)
?plo
plot(FR1, c(f1,f2,f3,f4),main="Scatterplot with average", xlab="Feed Rate",ylab="Finish", pch=19, col="red", ylim=range(c(min(finish), max(finish))))
lines(FR1, c(avg1,avg2,avg3,avg4), type = "b", col = "red", lwd = 3, pch = 15)

points(FR1, c(f5,f6,f7,f8),pch = 19,col="blue")
lines(FR1, c(avg5,avg6,avg7,avg8), type = "b", col = "blue", lwd = 3, pch = 15)

points(FR1, c(f9,f10,f11,f12),pch = 19,col="green")
lines(FR1, c(avg9,avg10,avg11,avg12), type = "b", col = "green", lwd = 3, pch = 15)

legend(x=0.21,y=80,c("Feed rate=0.2 in/min","Feed rate=0.25 in/min","Feed rate=0.30 in/min"),pch=c(19,19),col = c("red","blue","green"))
plot (model)
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))

f1 <- finish[1:12]
f2 <- finish[13:24]
f3 <- finish[25:36]

sdev <-c(rep(sd(f1),12),rep(sd(f2),12),rep(sd(f3),12)) 
avg <- c(ave(f1), ave(f2), ave(f3))

plot(FR, finish, ylim=range(c(avg-sdev, avg+sdev)), main="Scatterplot with average and error bar", xlab="Feed Rate ", ylab="Surface Finish", pch=19)

points(FR,avg,pch = 15,col="red")
arrows(FR, avg-sdev, FR, avg+sdev, length=0.05, angle=90, code=3, col="red")
lines(FR, avg, type = "b", col = "red", lwd = 3, pch = 15)
##?lines()
avg
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut  <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
block1 <- c(rep(seq(1,3),12))


dat1 <- cbind(finish,cut,FR,block1)
dat1 <- as.data.frame(dat1)

dat1$cut <- as.factor(dat1$cut)
dat1$FR <- as.factor(dat1$FR)
dat1$block1 <- as.factor(dat1$block1)
str(dat1)
model<-aov(finish~cut+FR+block1+cut*FR,data=dat1)
summary(model)

density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position+temp*position)
GAD::gad(model)
plot 

density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position)
GAD::gad(model)
plot (model)

dat1 <- c(50,49,50,50,48,51,52,52,51,51,51,51,53,50,50,54,52,51,49,51,50,48,50,51,48,49,48,48,49,48,52,50,50,52,50,50,51,51,51,51,50,50,52,50,49,53,48,50,50,51,50,51,48,49,47,46,49,46,47,48)
operator <- c(rep(c(rep(1,3),rep(2,3)),10))
operator
part <- c(rep(1,6),rep(2,6),rep(3,6),rep(4,6),rep(5,6),rep(6,6),rep(7,6),rep(8,6),rep(9,6),rep(10,6))
part
operator <- as.fixed(operator)
part <- as.random(part)
model<-aov(dat1~operator+part+operator*part)
GAD::gad(model)
plot (model)