# R example for Two-Factor Study with Unequal Sample Sizes  

# We use the growth rate data from Chapter 23 and the example in class.  

# Save the data file into a directory and 
# use the full path name:

growth.data <- read.table(file = "z:/My Documents/teaching/stat_705/growthdata.txt", 
header=FALSE, col.names = c('growthratechange', 'gender', 'bone', 'child'))

# attaching the data frame:

attach(growth.data)

# Making "gender" and "bone" factors:

gender <- factor(gender)
bone <- factor(bone)

# Fitting the ANOVA model and producing the ANOVA table:

# The standard analysis will produce the Type I SS
# (NOT what we want here!)

# growth.fit <- lm(growthratechange ~ gender + bone + gender:bone);
# anova(growth.fit)

# To get the Type III SS, we can install the "car" package in R:

# To install the "car" package:
# Go to "Packages" menu and choose
# "Install package(s) from CRAN"
# scroll to the "car" package and choose it
# An Internet connection is needed to install it.

# Then you can proceed with the code...

# Load the "car" package:

library(car)

growth.fit <- lm(growthratechange ~ gender + bone + gender:bone);

Anova(growth.fit, type=3)

# Note that this package's definition of Type III SS is slightly 
# different for SAS's (and the book's).  In this example, the Type III
# SS for "bone" is different from what SAS reports.

# Here is a way to get the Type III SS to match SAS in this example:
# (This is the long way, following the calculations on page 955-959 
# of the book.)

x1 <- ifelse(gender==1, 1, -1)
x2 <- ifelse(bone==1, 1, (ifelse(bone==3, -1, 0)) )
x3 <- ifelse(bone==2, 1, (ifelse(bone==3, -1, 0)) )

full.model <- lm(growthratechange ~ x1 + x2 + x3 + x1:x2 + x1:x3)
reduced.model.AB <- lm(growthratechange ~ x1 + x2 + x3)
reduced.model.A <- lm(growthratechange ~ x2 + x3 + x1:x2 + x1:x3)
reduced.model.B <- lm(growthratechange ~ x1 + x1:x2 + x1:x3)

# Gives SAS-style Type III SS list for A, B, and A*B
Type.III.SS.table <- rbind(
anova(reduced.model.A,full.model)[2,],
anova(reduced.model.B,full.model)[2,],
anova(reduced.model.AB,full.model)[2,] )

row.names(Type.III.SS.table) <- c("gender", "bone", "gender*bone")

print(Type.III.SS.table)

# When the data are unbalanced (as they are here), the correct SS to  
# look at are the Type III SS, not Type I.  So   
# The correct F* for the main-effect test for factor A is 0.74 and    
# the correct F* for the main-effect test for factor B is 12.89.      
# (Compare to the values obtained at the bottom of page 958 of the    
# book.)                                                              

# Our conclusions are:  There is no significant interaction between    
# gender and bone development; there is no significant effect on mean  
# growth rate change due to gender; there is a significant effect on   
# mean growth rate change due to bone development.     

# Getting the cell sample means:

tapply(growthratechange, gender:bone, mean)

# Getting the least squares row means:
# In this example, these are for the levels of "gender":

apply(matrix(tapply(growthratechange, gender:bone, mean), nrow=2, ncol=3, byrow=T), 1, mean)

# Getting the least squares column means:
# In this example, these are for the levels of "bone":

apply(matrix(tapply(growthratechange, gender:bone, mean), nrow=2, ncol=3, byrow=T), 2, mean)

# To do the Tukey procedure with unbalanced data, we must install the "multcomp" package in R:

# To install the "multcomp" package:
# Go to "Packages" menu and choose
# "Install package(s) from CRAN"
# scroll to the "multcomp" package and choose it
# An Internet connection is needed to install it.

# Then you can proceed with the code...

# Load the "multcomp" package:

library(multcomp) 

# Tukey multiple comparisons for comparing means of the different levels of bone:
# Note the adjusted P-values show which levels are significantly different:

summary(glht(aov(growth.fit), linfct = mcp(bone = 'Tukey')) )

# Tukey simultaneous CIs (family confidence level = 0.90):

confint(glht(aov(growth.fit), linfct = mcp(bone = 'Tukey')), level=0.90 )

# Using Tukey's procedure (family significance level 0.10), we may      
# check with of the simultaneous CIs contain zero.  We thereby find     
# that levels 1 and 3 of bone development have significantly different  
# mean growth rate change; levels 1 and 3 of bone development have      
# significantly different mean growth rate change; but levels 1 and 2   
# of bone development do not have significantly different mean growth   
# rate change.                                                                     

# ##################################################################**  

# We will stick to SAS for the analysis with Empty Cells.
# The standard R methods do not seem to carry over well in the empty-cell case.