# dishwash.R
library(lsmeans)
#library(phia)
#library(doBy)
library(car)
library(effects)
library(lme4)
library(nlme)
# get the data
dishwash<-read.table(file="dishwash.dat",header=T,
colClasses=c("factor","factor","numeric"))
head(dishwash)
is.factor(dishwash$block)
is.factor(dishwash$dishes)
is.ordered(dishwash$block)
is.ordered(dishwash$dishsoap)
#### part 1
# the intra-block analysis of a BIBD using ANOVA
m1<-aov(dishes~dishsoap+block,data=dishwash)
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)
Anova(m1,type=3)
# get the contrasts of interest
c1<-c(-1,-1,-1,-1,-1,-1,-1,-1,8)
c2<-c(1,1,1,1,-1,-1,-1,-1,0)
c3<-c(-3,-1,1,3,0,0,0,0,0)
c41<-c(1,-1,-1,1,0,0,0,0,0)
c42<-c(-1,3,-3,1,0,0,0,0,0)
c5<-c(0,0,0,0,-3,-1,1,3,0)
c61<-c(0,0,0,0,1,-1,-1,1,0)
c62<-c(0,0,0,0,-1,3,-3,1,0)
lsmeans(m1,specs=lsm~dishsoap,contr=list(lsm=list(control_vs_others=c1,
detergentI_vs_detergentII=c2,
linear_in_additive_detergentI=c3,
# nonlinear_in_additive_detergentI=list(c41,c42),
linear_in_additive_detergentII=c5
# nonlinear_in_additive_detergentII=list(c61,c62)
)))
# Unfortunately, the lsmeans() function does not seem to be able to
# test multiple degree of freedom hypotheses on contrasts, so the
# nonlinear_in_additive_detergentI and nonlinear_in_additive_detergentII lines above
# are commented out (it won't work, although
# one would hope that it would). Instead we'll have to use the
# glht() function in the multcomp package. The documentation on glht
# is hard to follow and the design of glht and related functions
# is very complex, but I eventually figured out how to use it to do what
# we need to do for this example.
# Here is how to do the nonlinear_in_additive_detergentI and
# nonlinear_in_additive_detergentII tests with glht:
K1 <- rbind( "quadratic detergentI" = c41,
"cubic detergentI" = c42)
K2 <- rbind( "quadratic detergentII" = c61,
"cubic detergentII" = c62)
summary(glht(m1, linfct = mcp(dishsoap = K1 )),test=Ftest())
summary(glht(m1, linfct = mcp(dishsoap = K2 )),test=Ftest())
# use package effects to get means instead of lsmeans for each level of
# dishsoap. These are inappropriate estimates of the treatment means in an
# incomplete block design like this one; they are produced here only for
# comparison with the ls means (or fitted means from the model), which are
# appropriate here.
effect("dishsoap",m1)
# get the profile plot
# create variables "add" and "deterg"
dsoapmeans <- lsmeans(m1,specs=~dishsoap)[[1]]
dsoapmeans$add <- c(3,2,1,0,3,2,1,0,0)
dsoapmeans$deterg <- c(1,1,1,1,2,2,2,2,3)
head(dsoapmeans)
plot(dsoapmeans$add[1:4],dsoapmeans$lsmean[1:4],ylim=c(0,30),type="b",
lty=1,pch=1,col=1,
main="Mean # plates washed by amt of additive",ylab="Mean # plates washed",
xlab="Amount of additive")
lines(dsoapmeans$add[5:8],dsoapmeans$lsmean[5:8],type="b",lty=2,pch=2,col=2)
points(dsoapmeans$add[9],dsoapmeans$lsmean[9],type="b",pch=3,col=3)
legend(locator(1),legend=c("Deterg I","Deterg II", "Control"),
lty=c(1,2,3),pch=c(1,2,3),col=c(1,2,3))
#### part 2 mixed models (combined intra- and inter-block analysis)
m2<-lmer(dishes~dishsoap+(1|block),data=dishwash)
summary(m2)
# use the same contrasts
lsmeans(m2,specs=lsm~dishsoap,contr=list(lsm=list(control_vs_others=c1,
detergentI_vs_detergentII=c2,
linear_in_additive_detergentI=c3,
# nonlinear_in_additive_detergentI=list(c41,c42),
linear_in_additive_detergentII=c5
# nonlinear_in_additive_detergentII=list(c61,c62)
)))
# summary(glht(m2, linfct = mcp(dishsoap = K1 )),test=Ftest())
# summary(glht(m2, linfct = mcp(dishsoap = K2 )),test=Ftest())
# glht doesn't work for lmer, so refit the model using the lme function in
# the nlme package
# First create a groupedData data frame (same as the dishwash data frame
# but adds a grouping structure to indicate that data are grouped by block)
dish2 <- groupedData(dishes~dishsoap|block, data=dishwash)
# fit the linear mixed model (blocks random, treatments fixed)
m2a<-lme(dishes~dishsoap,data=dish2, random= ~1|block)
summary(m2a)
# Type 3 F test for dishsoap given by the following command
anova(m2a,type="marginal")
# Now test contrasts:
# the test statistics given below are labeled "z.ratio" and are
# treated as Z tests, not t tests. Nevertheless, for these 1 df
# contrasts, their squares are the F statistics given by PROC MIXED
# in SAS
lsmeans(m2a,specs=lsm~dishsoap,contr=list(lsm=list(control_vs_others=c1,
detergentI_vs_detergentII=c2,
linear_in_additive_detergentI=c3,
# nonlinear_in_additive_detergentI=list(c41,c42),
linear_in_additive_detergentII=c5
# nonlinear_in_additive_detergentII=list(c61,c62)
)))
# similarly, the tests for the 2 df contrasts given below are large sample
# chisquare tests
# not F tests. However the relationship between these chisq test statistics
# and the F test statistics given by PROC MIXED is simple
# (F test stat)*(num df) = chisq test stat.
summary(glht(m2a, linfct = mcp(dishsoap = K1 )),test=Ftest())
summary(glht(m2a, linfct = mcp(dishsoap = K2 )),test=Ftest())