####################################Chapter 23, TA01, unbalanced design###############################
growth<-read.table("~/Desktop/jenn/teaching/stat445545/data/CH23TA01.txt",col.names=c("y","A","B","obs"))
n<-nrow(growth)
n
## [1] 14
growth
## y A B obs
## 1 1.4 1 1 1
## 2 2.4 1 1 2
## 3 2.2 1 1 3
## 4 2.1 1 2 1
## 5 1.7 1 2 2
## 6 0.7 1 3 1
## 7 1.1 1 3 2
## 8 2.4 2 1 1
## 9 2.5 2 2 1
## 10 1.8 2 2 2
## 11 2.0 2 2 3
## 12 0.5 2 3 1
## 13 0.9 2 3 2
## 14 1.3 2 3 3
growth$A<-factor(growth$A)
growth$B<-factor(growth$B)
attach(growth)
#summary statistics
tapply(y,A,mean) #row means
## 1 2
## 1.657143 1.628571
tapply(y,B, mean) #column means
## 1 2 3
## 2.10 2.02 0.90
aggregate(y~A*B, data=growth,mean) #cell means
## A B y
## 1 1 1 2.0
## 2 2 1 2.4
## 3 1 2 1.9
## 4 2 2 2.1
## 5 1 3 0.9
## 6 2 3 0.9
aggregate(y~A*B, data=growth,length) #cell size, unequal
## A B y
## 1 1 1 3
## 2 2 1 1
## 3 1 2 2
## 4 2 2 3
## 5 1 3 2
## 6 2 3 3
##View interaction plots, also called profile plots, page 955
interaction.plot(B,A,y,type='b',
col=1:2, pch=1:2)
#full model
#adding restrictions, A=contr.sum,B=contr.sum
myfit<-lm(y~A+B+A:B,contrasts = c(A=contr.sum,B=contr.sum))
## The contrast statement above must be included identifying
## each main effect with "contr.sum" in order for the correct
## Type III SS to be computed.
#without restriction
myfitwot<-lm(y~A+B+A:B)
#type I
anova(myfit) #type I SS
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 0.0029 0.00286 0.0176 0.897785
## B 2 4.3960 2.19800 13.5262 0.002713 **
## A:B 2 0.0754 0.03771 0.2321 0.798034
## Residuals 8 1.3000 0.16250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(myfitwot) #type I SS
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 0.0029 0.00286 0.0176 0.897785
## B 2 4.3960 2.19800 13.5262 0.002713 **
## A:B 2 0.0754 0.03771 0.2321 0.798034
## Residuals 8 1.3000 0.16250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(car)
## Loading required package: carData
#type III
Anova(myfit,type=3) #type III SS
## Anova Table (Type III tests)
##
## Response: y
## Sum Sq Df F value Pr(>F)
## (Intercept) 34.680 1 213.4154 4.729e-07 ***
## A 0.120 1 0.7385 0.415160
## B 4.190 2 12.8914 0.003145 **
## A:B 0.075 2 0.2321 0.798034
## Residuals 1.300 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(myfitwot,type=3) #type III SS
## Anova Table (Type III tests)
##
## Response: y
## Sum Sq Df F value Pr(>F)
## (Intercept) 12.0000 1 73.8462 2.6e-05 ***
## A 0.1200 1 0.7385 0.41516
## B 1.6171 2 4.9758 0.03944 *
## A:B 0.0754 2 0.2321 0.79803
## Residuals 1.3000 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#reduced model
myfit2<-lm(y~A+B,contrasts = c(A=contr.sum,B=contr.sum))
anova(myfit2) #type I SS
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 0.0029 0.00286 0.0208 0.8882630
## B 2 4.3960 2.19800 15.9805 0.0007687 ***
## Residuals 10 1.3754 0.13754
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(myfit2,type=3) #type III SS
## Anova Table (Type III tests)
##
## Response: y
## Sum Sq Df F value Pr(>F)
## (Intercept) 38.855 1 282.493 1.166e-08 ***
## A 0.093 1 0.673 0.4311159
## B 4.396 2 15.980 0.0007687 ***
## Residuals 1.375 10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
####################################Bakery data, balanced design###############################
bakery<- read.table(file="~/Desktop/jenn/teaching/stat445545/data/CH19TA07.txt",col.names=c("y","height","width","obs"))
nt<-nrow(bakery)
nt
## [1] 12
attach(bakery)
## The following objects are masked from growth:
##
## obs, y
bakery
## y height width obs
## 1 47 1 1 1
## 2 43 1 1 2
## 3 46 1 2 1
## 4 40 1 2 2
## 5 62 2 1 1
## 6 68 2 1 2
## 7 67 2 2 1
## 8 71 2 2 2
## 9 41 3 1 1
## 10 39 3 1 2
## 11 42 3 2 1
## 12 46 3 2 2
height<-factor(height)
width<-factor(width)
aggregate(y~height*width, data=bakery,length) #cell size, balanced design
## height width y
## 1 1 1 2
## 2 2 1 2
## 3 3 1 2
## 4 1 2 2
## 5 2 2 2
## 6 3 2 2
myfitbakery<-lm(y~height+width+height*width,contrasts = c(height=contr.sum,width=contr.sum))
anova(myfitbakery) #Type I SS
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## height 2 1544 772.00 74.7097 5.754e-05 ***
## width 1 12 12.00 1.1613 0.3226
## height:width 2 24 12.00 1.1613 0.3747
## Residuals 6 62 10.33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(myfitbakery,type=2) #Type II SS
## Anova Table (Type II tests)
##
## Response: y
## Sum Sq Df F value Pr(>F)
## height 1544 2 74.7097 5.754e-05 ***
## width 12 1 1.1613 0.3226
## height:width 24 2 1.1613 0.3747
## Residuals 62 6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(myfitbakery,type=3) #Type III SS
## Anova Table (Type III tests)
##
## Response: y
## Sum Sq Df F value Pr(>F)
## (Intercept) 31212 1 3020.5161 2.437e-09 ***
## height 1544 2 74.7097 5.754e-05 ***
## width 12 1 1.1613 0.3226
## height:width 24 2 1.1613 0.3747
## Residuals 62 6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#type I, II, and III SS are all the same