Problem 4.3
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.
obs <- c(73,68,74,71,67,73,67,75,72,70,75,68,78,73,68,73,71,75,75,69)
trt <- c(rep(1,5),rep(2,5),rep(3,5),rep(4,5))
block1 <- c(rep(seq(1,5),4))
trt<-as.fixed(trt)
block1<-as.fixed(block1)
model<-lm(obs~trt+block1)
gad(model)
## Analysis of Variance Table
##
## Response: obs
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 3 12.95 4.317 2.3761 0.1211
## block1 4 157.00 39.250 21.6055 2.059e-05 ***
## Residual 12 21.80 1.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Block1 herein denotes Bolt.
\(\overline{y_{ij}}=\mu_{i}+\beta_{j}+\epsilon_{ij}\)
\(\overline{y_{ij}}=\mu+\tau_{i}+\beta_{j}+\epsilon_{ij}\)
Corresponding to:
\(\tau_{i}=0\) null hypothesis: \(H_0: \mu_{1}=\mu_{2}=\cdots =\mu_{i}=\mu\)
\(\tau_{i}\ne 0\) alternative hypothesis: \(H_1:\) at least one \(\mu_{i}\) differs
Ans: P-value=0.1211 > 0.05, so we do not reject null hypothesis (i.e., means are equal).
————————————————————————————————
Problem 4.16
meanI1=mean(obs[1:5],)
meanI2=mean(obs[6:10],)
meanI3=mean(obs[11:15],)
meanI4=mean(obs[16:20],)
gm=mean(c(meanI1,meanI2,meanI3,meanI4))
tau1=meanI1-gm
tau2=meanI2-gm
tau3=meanI3-gm
tau4=meanI4-gm
tau1
## [1] -1.15
tau2
## [1] -0.35
tau3
## [1] 0.65
tau4
## [1] 0.85
meanJ1=mean(c(obs[1],obs[6],obs[11],obs[16]))
meanJ2=mean(c(obs[2],obs[7],obs[12],obs[17]))
meanJ3=mean(c(obs[3],obs[8],obs[13],obs[18]))
meanJ4=mean(c(obs[4],obs[9],obs[14],obs[19]))
meanJ5=mean(c(obs[5],obs[10],obs[15],obs[20]))
gm=mean(c(meanJ1,meanJ2,meanJ3,meanJ4,meanJ5))
beta1=meanJ1-gm
beta2=meanJ2-gm
beta3=meanJ3-gm
beta4=meanJ4-gm
beta5=meanJ5-gm
beta1
## [1] 1.75
beta2
## [1] -3.25
beta3
## [1] 3.75
beta4
## [1] 1
beta5
## [1] -3.25
\(\tau_{1}\)=-1.15, \(\tau_{2}\)=-0.35, \(\tau_{3}\)=0.65, \(\tau_{2}\)=0.85; \(\beta_{1}\)=1.75, \(\beta_{2}\)=-3.25, \(\beta_{3}\)=3.75, \(\beta_{1}\)=1, \(\beta_{1}\)=-3.25.
————————————————————————————————
Problem 4.22
obs <- c(8, 7, 1, 7, 3, 11, 2, 7, 3, 8, 4, 9, 10, 1, 5, 6, 8, 6, 6, 10, 4, 2, 3, 8, 8)
trt <- c(1,2,4,3,5,3,5,1,4,2,2,1,3,5,4,4,3,5,2,1,5,4,2,1,3)
block1 <- c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5))
block2 <- c(rep(seq(1,5),5))
dat1 <- cbind(obs,trt,block1,block2)
dat1 <- as.data.frame(dat1)
str(dat1)
## 'data.frame': 25 obs. of 4 variables:
## $ obs : num 8 7 1 7 3 11 2 7 3 8 ...
## $ trt : num 1 2 4 3 5 3 5 1 4 2 ...
## $ block1: num 1 1 1 1 1 2 2 2 2 2 ...
## $ block2: num 1 2 3 4 5 1 2 3 4 5 ...
dat1$trt<-as.factor(dat1$trt)
dat1$block1<-as.factor(dat1$block1)
dat1$block2<-as.factor(dat1$block2)
str(dat1)
## 'data.frame': 25 obs. of 4 variables:
## $ obs : num 8 7 1 7 3 11 2 7 3 8 ...
## $ trt : Factor w/ 5 levels "1","2","3","4",..: 1 2 4 3 5 3 5 1 4 2 ...
## $ block1: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ block2: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
\(\tau_{i}=0\) null hypothesis: \(H_0: \mu_{1}=\mu_{2}=\cdots =\mu_{i}=\mu\)
Ans: Since p-value=0.000488<0.05, we do reject the null hypothesis (i.e., at least one means differs). The factor of interest here is ingredients, batch and day are recognized as blocks, which are nuisance and not the subjects we are studing on.
we have 5 samples only, it is difficult to check normality and variance. Dr. Matis said we treat this case as only an example to practice Latin Square, Assumption checks are therefore ignored herein.
————————————————————————————————
Raw Code:
library(GAD)
obs <- c(73,68,74,71,67,73,67,75,72,70,75,68,78,73,68,73,71,75,75,69)
trt <- c(rep(1,5),rep(2,5),rep(3,5),rep(4,5))
block1 <- c(rep(seq(1,5),4))
trt<-as.fixed(trt)
block1<-as.fixed(block1)
model<-lm(obs~trt+block1)
gad(model)
meanI1=mean(obs[1:5],)
meanI2=mean(obs[6:10],)
meanI3=mean(obs[11:15],)
meanI4=mean(obs[16:20],)
gm=mean(c(meanI1,meanI2,meanI3,meanI4))
tau1=meanI1-gm
tau2=meanI2-gm
tau3=meanI3-gm
tau4=meanI4-gm
tau1
tau2
tau3
tau4
meanJ1=mean(c(obs[1],obs[6],obs[11],obs[16]))
meanJ2=mean(c(obs[2],obs[7],obs[12],obs[17]))
meanJ3=mean(c(obs[3],obs[8],obs[13],obs[18]))
meanJ4=mean(c(obs[4],obs[9],obs[14],obs[19]))
meanJ5=mean(c(obs[5],obs[10],obs[15],obs[20]))
gm=mean(c(meanJ1,meanJ2,meanJ3,meanJ4,meanJ5))
beta1=meanJ1-gm
beta2=meanJ2-gm
beta3=meanJ3-gm
beta4=meanJ4-gm
beta5=meanJ5-gm
beta1
beta2
beta3
beta4
beta5
obs <- c(8, 7, 1, 7, 3, 11, 2, 7, 3, 8, 4, 9, 10, 1, 5, 6, 8, 6, 6, 10, 4, 2, 3, 8, 8)
trt <- c(1,2,4,3,5,3,5,1,4,2,2,1,3,5,4,4,3,5,2,1,5,4,2,1,3)
block1 <- c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5))
block2 <- c(rep(seq(1,5),5))
dat1 <- cbind(obs,trt,block1,block2)
dat1 <- as.data.frame(dat1)
str(dat1)
dat1$trt<-as.factor(dat1$trt)
dat1$block1<-as.factor(dat1$block1)
dat1$block2<-as.factor(dat1$block2)
str(dat1)
