# bleach.R
library(lsmeans)
library(car)
library(phia)
library(doBy)
# get the data
bleach<-read.table(file="bleach.dat",header=T)
head(bleach)
bleach$amtfac <- factor(bleach$amtblch)
is.ordered(bleach$amtfac)
is.ordered(bleach$staintyp)
######PART 1
# fit the two-way anova model
m1 <- aov(time~amtfac*staintyp,data=bleach)
summary(m1)
# note that the summary function gives F tests based on type I
# or sequential sums of squares. To get F tests based on the type II and III
# SS, we need the Anova() function from the car package.
Anova(m1,type=2)
# the Anova() function from the car package claims to produce Type III
# SSs and F tests, but if you run the line below you'll see that you get
# results quite different from those provided by SAS.
Anova(m1,type=3)
# The reason for this discrepancy is realted to the method of computing
# type III SSs in the Anova() function. That method depends upon the
# parameterization of the model. R and SAS parameterize factors differently.
# SAS uses a "non-full rank" parameterization (aka an overparameterization),
# whereas R removes the overparameterization. The user has some control over
# how that overparameterization is removed. That is, the user has some control
# over the parameterization that R uses. For the Anova(,type=3) function to
# produce the "right" type III SSs (the ones that agree with SAS and test
# hypotheses on marginal means as described in our lecture notes), one
# has to choose a sum-to-zero parameterization of factor effects. This is
# not the default in R. Therefore it has to be changed and the model refit.
# Change the handling of unordered factors (like amtfac and staintyp) to
# use the sum-to-zero constraints:
op <- options(contrasts=c("contr.sum", "contr.poly"))
options()$contrasts # now you can see that they've been changed
# Now refit the model with the new parameterization:
m1a <- aov(time~amtfac*staintyp,data=bleach)
# now check the type 1 SSs (they haven't changed):
summary(m1a)
# and check the type II SSs (they haven't changed):
Anova(m1a,type=2)
# and check the type III SSs (they have changed and now agree with SAS):
Anova(m1a,type=3)
######PART 2
# Now get the marginal and joint lsmeans:
lsmeans(m1, specs= list(~staintyp,~amtfac,~amtfac:staintyp))
# get profile plots (uses package phia)
( trtmns <- interactionMeans(m1) )
plot(trtmns)
# alternatively can use function lsmip from lsmeans package
lsmip(m1,amtfac~staintyp)
lsmip(m1,staintyp~amtfac)
######PARTS 3 & 4
# check resids versus fitted plot and normal QQ plot of resids
dev.off(which=dev.cur())
plot(m1,which=1)
plot(m1,which=2)
# looks like non-constant variance. Use regression method to estimate
# transformation:
# to get summary statistics (mean, sd) by treatment we can use the
# summaryBy function in the doBy package
(sumStats <- summaryBy(time ~ amtfac + staintyp, data = bleach,
FUN = function(x) { c(m = mean(x), s = sd(x)) } ) )
# produces time.m and time.s for each
# combination of the levels of amtfac and staintyp
#regress log(SD) on log(mean) to estimate power transformation
lm(I(log(time.s)) ~ I(log(time.m)),data=sumStats)
# slope is .53 which is approximately 0.5, so this suggests
# taking y to the 1-.5=0.5 power (sqrt(y))
# creat sqrttime=sqrt(time) in data bleach
bleach$sqrttime<-sqrt(bleach$tim)
# fit Model 2: 2-way ANOVA model after taking sqrt transformation for time
m2 <- aov(sqrttime~amtfac*staintyp,data=bleach)
summary(m2)
# Type III anova table:
Anova(m2,type=3)
# Now get the marginal and joint lsmeans:
lsmeans(m2, specs= list(~staintyp,~amtfac,~amtfac:staintyp))
# get profile plots (uses package phia)
( trtmns2 <- interactionMeans(m2) )
plot(trtmns2)
# alternatively can use function lsmip from lsmeans package
lsmip(m2,amtfac~staintyp)
lsmip(m2,staintyp~amtfac)
# check resids versus fitted plot and normal QQ plot of resids
dev.off(which=dev.cur())
plot(m2,which=1)
plot(m2,which=2)
# These look better.
######PART 5
# below are the linear and quadratic contrasts for amt of bleach within each
# stain type:
c1<-c(-1,0,1,0,0,0,0,0,0)
c2<-c(1,-2,1,0,0,0,0,0,0)
c3<-c(0,0,0,-1,0,1,0,0,0)
c4<-c(0,0,0,1,-2,1,0,0,0)
c5<-c(0,0,0,0,0,0,-1,0,1)
c6<-c(0,0,0,0,0,0,-1,2,1)
lsmeans(m2,specs=lsm~amtfac:staintyp,
contr=list(lsm=list(linear.bleach.stain1=c1,
nonlinear.bleach.stain1=c2,
linear.bleach.stain2=c3,
nonlinear.bleach.stain2=c4,
linear.bleach.stain3=c5,
nonlinear.bleach.stain3=c6)))
options(op) # reset contrasts to their original values