HW 10 Final
Sahil Gaikwad
11/7/2021
knitr::opts_chunk$set(echo = TRUE)Question 5.4
Hypothesis to be tested
Null Hypothesis(Ho): \(\alpha_i=0\)
Alternative Hypothesis(Ha): \(\alpha_i\neq0\)
Null Hypothesis(Ho): \(\beta_j=0\)
Alternative Hypothesis(Ha): \(\beta_j\neq0\)
Null Hypothesis(Ho): \(\alpha\beta_{ij}=0\)
Alternative Hypothesis(Ha):\(\alpha\beta_{ij}\neq0\)
Model Equation:
\(Y_{ijk} = \mu_i + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ijk}\)
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.library(agricolae)
depth<-(rep(seq(1,4),9))
rate<-c(rep(1,12),rep(2,12),rep(3,12))
obs<-c(74,79,82,99,64,68,88,104,60,73,92,96,92,98,99,104,86,104,108,110,88,88,95,99,99,104,108,114,98,99,110,111,102,95,99,107)
df<-data.frame(depth,rate,obs)
df$depth<-as.fixed(df$depth)
df$rate<-as.fixed(df$rate)
df## depth rate obs
## 1 1 1 74
## 2 2 1 79
## 3 3 1 82
## 4 4 1 99
## 5 1 1 64
## 6 2 1 68
## 7 3 1 88
## 8 4 1 104
## 9 1 1 60
## 10 2 1 73
## 11 3 1 92
## 12 4 1 96
## 13 1 2 92
## 14 2 2 98
## 15 3 2 99
## 16 4 2 104
## 17 1 2 86
## 18 2 2 104
## 19 3 2 108
## 20 4 2 110
## 21 1 2 88
## 22 2 2 88
## 23 3 2 95
## 24 4 2 99
## 25 1 3 99
## 26 2 3 104
## 27 3 3 108
## 28 4 3 114
## 29 1 3 98
## 30 2 3 99
## 31 3 3 110
## 32 4 3 111
## 33 1 3 102
## 34 2 3 95
## 35 3 3 99
## 36 4 3 107model<-aov(obs~df$depth+df$rate+df$depth*df$rate,data=df)
GAD::gad(model)## Analysis of Variance Table
##
## Response: obs
## Df Sum Sq Mean Sq F value Pr(>F)
## df$depth 3 2125.11 708.37 24.6628 1.652e-07 ***
## df$rate 2 3160.50 1580.25 55.0184 1.086e-09 ***
## df$depth:df$rate 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 ' ' 1plot(model,col="darkred")
From the Normal probability plots we can conclude that model is adequate
To find point estimates of the mean surface finish at each feed rate
?model.tables## starting httpd help server ... donemodel.tables(model,type = "means",cterms = "df$rate")## Tables of means
## Grand mean
##
## 94.33333
##
## df$rate
## df$rate
## 1 2 3
## 81.58 97.58 103.83p-value of df$depth = 1.652e-07
p-value of df$rate= 1.086e-09
p-value of df\(depth:df\)rate= 0.01797
Question 5.34
library(GAD)
library(agricolae)
depth<-(rep(seq(1,4),9))
rate<-c(rep(1,12),rep(2,12),rep(3,12))
blocks<-c(rep(1,4),rep(2,4),rep(3,4),rep(1,4),rep(2,4),rep(3,4),rep(1,4),rep(2,4),rep(3,4))
obs<-c(74,79,82,99,64,68,88,104,60,73,92,96,92,98,99,104,86,104,108,110,88,88,95,99,99,104,108,114,98,99,110,111,102,95,99,107)
df<-data.frame(depth,rate,blocks,obs)
df$depth<-as.fixed(df$depth)
df$rate<-as.fixed(df$rate)
df$blocks<-as.fixed(df$blocks)
df## depth rate blocks obs
## 1 1 1 1 74
## 2 2 1 1 79
## 3 3 1 1 82
## 4 4 1 1 99
## 5 1 1 2 64
## 6 2 1 2 68
## 7 3 1 2 88
## 8 4 1 2 104
## 9 1 1 3 60
## 10 2 1 3 73
## 11 3 1 3 92
## 12 4 1 3 96
## 13 1 2 1 92
## 14 2 2 1 98
## 15 3 2 1 99
## 16 4 2 1 104
## 17 1 2 2 86
## 18 2 2 2 104
## 19 3 2 2 108
## 20 4 2 2 110
## 21 1 2 3 88
## 22 2 2 3 88
## 23 3 2 3 95
## 24 4 2 3 99
## 25 1 3 1 99
## 26 2 3 1 104
## 27 3 3 1 108
## 28 4 3 1 114
## 29 1 3 2 98
## 30 2 3 2 99
## 31 3 3 2 110
## 32 4 3 2 111
## 33 1 3 3 102
## 34 2 3 3 95
## 35 3 3 3 99
## 36 4 3 3 107nmodel<-aov(obs~df$depth+df$rate+df$blocks+df$depth*df$rate,data=df)
GAD::gad(nmodel)## Analysis of Variance Table
##
## Response: obs
## Df Sum Sq Mean Sq F value Pr(>F)
## df$depth 3 2125.11 708.37 30.6373 4.893e-08 ***
## df$rate 2 3160.50 1580.25 68.3463 3.635e-10 ***
## df$blocks 2 180.67 90.33 3.9069 0.035322 *
## df$depth:df$rate 6 557.06 92.84 4.0155 0.007258 **
## Residual 22 508.67 23.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(nmodel,col="darkred")
From the results of the ANOVA and gad model test, the p-value of interaction of factorA and FactorB is 0.007258 As p-value is lesser than the alpha 0.05, we conclude that we reject the null hypothesis at 0.05 level of significance.
From the residual plot we can conclude that the our data is normally distributed except few out liners
From residuals vs fitted plots we can see that variances are widely spread, hence we conclude that variances are not constant
To find the variance component of blocking, according to the formula, (MSBlocking-MSE)/(IXJ)
varcomp<-((90.33-23.12)/(4*3))
varcomp## [1] 5.600833Hence variance component of blocking is 5.600833.And Blocking does not affect the experiment
Question 13.5
Model Equation for two factor interaction:
\(Y_{ijk} = \mu_i + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ijk}\)
Where
\(\alpha_i\) is Main Effects of Factor A
\(\beta_j\) is Main Effects of Factor B
\(\alpha\beta_{ij}\) is Interaction effects of Factors A and Factors B
Hypothesis to be tested
Null Hypothesis(Ho): \(\alpha_i=0\) For all i
Alternative Hypothesis(Ha): \(\alpha_i\neq0\) For some i
Null Hypothesis(Ho): \(\sigma^2\beta =0\)
Alternative Hypothesis(Ha): \(\sigma^2\beta\neq0\)
Null Hypothesis(Ho): \(\sigma^2\alpha\beta=0\)
Alternative Hypothesis(Ha): \(\alpha\beta_{ij}\neq0\)
we will test first highest order interaction effects hypothesis
Null Hypothesis(Ho): \(\sigma^2\alpha\beta=0\)
Alternative Hypothesis(Ha): \(\alpha\beta_{ij}\neq0\)
library(GAD)
library(agricolae)
temp<-rep(seq(1,3),6)
pos<-c(rep(1,9),rep(2,9))
obs<-c(570,1063,565,565,1080,510,583,1043,590,528,988,526,547,1026,538,521,1004,532)
df<-data.frame(temp,pos,obs)
df$temp<-as.fixed(df$temp)
df$pos<-as.random(df$pos)
df## temp pos obs
## 1 1 1 570
## 2 2 1 1063
## 3 3 1 565
## 4 1 1 565
## 5 2 1 1080
## 6 3 1 510
## 7 1 1 583
## 8 2 1 1043
## 9 3 1 590
## 10 1 2 528
## 11 2 2 988
## 12 3 2 526
## 13 1 2 547
## 14 2 2 1026
## 15 3 2 538
## 16 1 2 521
## 17 2 2 1004
## 18 3 2 532model<-aov(df$obs~df$pos+df$temp+df$pos*df$temp,data=df)
GAD::gad(model)## Analysis of Variance Table
##
## Response: df$obs
## Df Sum Sq Mean Sq F value Pr(>F)
## df$pos 1 7160 7160 15.998 0.0017624 **
## df$temp 2 945342 472671 1155.518 0.0008647 ***
## df$pos:df$temp 2 818 409 0.914 0.4271101
## Residual 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model,col="darkred")
?interaction.plot
interaction.plot(df$temp, df$pos,obs, type = "l", ylab = "obs", xlab = "temp", main = "Interaction Plot", col = c("darkred", "black"))
From the residual plot we can conclude that the our data is normally distributed except few out liners
From residuals vs fitted plots we can see that variances are widely spread, hence we conclude that variances are not constant
Now to Estimation the model components
MSA<-7160
MSB<-472671
MSAB<-409
MSE<-448
var_A<-(MSA-MSAB)/(3*3)
var_B<-(MSB-MSAB)/(2*3)
var_AB<-(MSAB-MSE)/(3)
var_E<-448
var<-c(var_A,var_B,var_AB,var_E)
source<-c("A","B","AB","Error")
df<-data.frame(var,source)
df## var source
## 1 750.1111 A
## 2 78710.3333 B
## 3 -13.0000 AB
## 4 448.0000 ErrorFrom the results of the ANOVA and gad model test, the p-value of Interaction between FactorA and FactorB is 0.4271101 and for factorA and factorB is 0.0017624 and 0.0008647 respectively.
Question 13.6
Model Equation for two factor interaction:
\(Y_{ijk} = \mu_i + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ijk}\)
Where
\(\alpha_i\) is Main Effects of Factor A
\(\beta_j\) is Main Effects of Factor B
\(\alpha\beta_{ij}\) is Interaction effects of Factors A and Factors B
Hypothesis to be tested
Null Hypothesis(Ho): \(\alpha_i=0\) For all i
Alternative Hypothesis(Ha): \(\alpha_i\neq0\) For some i
Null Hypothesis(Ho): \(\sigma^2\beta =0\)
Alternative Hypothesis(Ha): \(\sigma^2\beta\neq0\)
Null Hypothesis(Ho): \(\sigma^2\alpha\beta=0\)
Alternative Hypothesis(Ha): \(\alpha\beta_{ij}\neq0\)
we will test first highest order interaction effects hypothesis
Null Hypothesis(Ho): \(\sigma^2\alpha\beta=0\)
Alternative Hypothesis(Ha): \(\alpha\beta_{ij}\neq0\)
library(GAD)
library(agricolae)
operator<-rep(c(rep(1,3),rep(2,3)),10)
partnumber<-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))
measure<-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)
df<-data.frame(operator,partnumber,measure)
df## operator partnumber measure
## 1 1 1 50
## 2 1 1 49
## 3 1 1 50
## 4 2 1 50
## 5 2 1 48
## 6 2 1 51
## 7 1 2 52
## 8 1 2 52
## 9 1 2 51
## 10 2 2 51
## 11 2 2 51
## 12 2 2 51
## 13 1 3 53
## 14 1 3 50
## 15 1 3 50
## 16 2 3 54
## 17 2 3 52
## 18 2 3 51
## 19 1 4 49
## 20 1 4 51
## 21 1 4 50
## 22 2 4 48
## 23 2 4 50
## 24 2 4 51
## 25 1 5 48
## 26 1 5 49
## 27 1 5 48
## 28 2 5 48
## 29 2 5 49
## 30 2 5 48
## 31 1 6 52
## 32 1 6 50
## 33 1 6 50
## 34 2 6 52
## 35 2 6 50
## 36 2 6 50
## 37 1 7 51
## 38 1 7 51
## 39 1 7 51
## 40 2 7 51
## 41 2 7 50
## 42 2 7 50
## 43 1 8 52
## 44 1 8 50
## 45 1 8 49
## 46 2 8 53
## 47 2 8 48
## 48 2 8 50
## 49 1 9 50
## 50 1 9 51
## 51 1 9 50
## 52 2 9 51
## 53 2 9 48
## 54 2 9 49
## 55 1 10 47
## 56 1 10 46
## 57 1 10 49
## 58 2 10 46
## 59 2 10 47
## 60 2 10 48df$partnumber<-as.random(df$partnumber)
df$operator<-as.fixed(df$operator)
model<-aov(measure~partnumber+operator+partnumber*operator,data=df)
GAD::gad(model)## Analysis of Variance Table
##
## Response: measure
## Df Sum Sq Mean Sq F value Pr(>F)
## partnumber 9 99.017 11.0019 7.3346 3.216e-06 ***
## operator 1 0.417 0.4167 0.6923 0.4269
## partnumber:operator 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 ' ' 1plot(model,col="darkred")
From the residual plot we can conclude that the our data is normally distributed except few out liners
From residuals vs fitted plots we can see that variances are widely spread, hence we conclude that variances are not constant
Now to estimate the model components
MSA_partnumber<-11.0019
MSB_operator<-0.4167
MSAB_partnumberoperator<-0.6019
MSE<-1.5000
var_A<-(MSA_partnumber-MSAB_partnumberoperator)/(2*3)
var_B<-(MSB_operator-MSAB_partnumberoperator)/(10*3)
var_AB<-(MSAB_partnumberoperator-MSE)/(3)
var_E<-1.5000
var<-c(var_A,var_B,var_AB,var_E)
source<-c("A","B","AB","Error")
df<-data.frame(var,source)
df## var source
## 1 1.733333333 A
## 2 -0.006173333 B
## 3 -0.299366667 AB
## 4 1.500000000 Error