# redwing1.R
library(lsmeans)
library(car)
library(multcomp)
# get the data
# you may want to change the path to where you put the data set
redwing<-read.table(file="redwing1.dat",header=T)
head(redwing)
is.factor(redwing$block)
redwing$blockfac <- factor(redwing$block)
# fit the linear model
m1<-lm(oil~treat+blockfac,data=redwing)
summary(m1)
# get the lsmeans for treat
lsmeans(m1,specs=~treat)
# test the contrasts of interest
# Here we use the lsmeans() function in package lsmeans
c1<-c(1,1,1,1,1,-5)
c2<-c(4,-1,-1,-1,-1,0)
c3<-c(0,1,1,1,-3,0)
c4<-c(0,1,-1,0,0,0)
c5<-c(0,1,1,-2,0,0)
c61<-c(4,-1,-1,-1,-1,0)
c62<-c(0,1,1,1,-3,0)
c63<-c(0,1,-1,0,0,0)
c64<-c(0,1,1,-2,0,0)
c7<-c(0,0,0,0,0,1)
lsmeans(m1,specs=lsm~factor(treat),contr=list(lsm=list(treat.v.control=c1,
seedling.v.other=c2,
bloom.v.ripening=c3,
within.bloom.i=c4,
within.bloom.ii=c5,
#Inoculation.timing=list(c61,c62,c63,c64)
control.mean=c7)))
# Unfortunately, the lsmeans() function does not seem to be able to
# test multiple degree of freedom hypotheses on contrasts, so the
# Inoculation.timing line above is 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 inoculation timing test with glht:
K1 <- rbind( "seed.v.other" = c61,
"blm.v.rip" = c62,
"w/in blm (i)" = c63,
"w/in blm (ii)" = c64)
summary(glht(m1, linfct = mcp(treat = K1 )),test=Ftest())
# And here are the Dunnett tests of each pairwise contrast with the best
# treatment, which is treatment d.
K2 <- rbind( "a - d" = c( 1, 0, 0, -1, 0, 0),
"b - d" = c( 0, 1, 0, -1, 0, 0),
"c - d" = c( 0, 0, 1, -1, 0, 0),
"e - d" = c( 0, 0, 0, -1, 1, 0),
"f - d" = c( 0, 0, 0, -1, 0, 1))
summary(glht(m1, linfct = mcp(treat = K2 ),
alternative="less"))
# Alternatively, the linearHypothesis function in the car package
# can be used to the multiple df test on inoculation timing.
# To do this it is easier
# to refit the model with the parameterization y_ij=mu_i +beta_j +e_ij
m1a <-lm(oil~treat+blockfac-1,data=redwing)
summary(m1a)
# in the parameterization R uses, the estimated mu_i's are actually the
# treatment means for block 1. This is because the constraint that R uses
# on the beta_j's is that beta_1=0, not that the sum of the beta_j's is 0.
# However, in this additive model with no
# treatment by block interactions, any contrasts among the mu_i's in any
# given block with be the same as within any other block. Therefore, we
# can go ahead and do our contrasts based on the previous model.
c61a<-c(4,-1,-1,-1,-1,0,0,0,0)
c62a<-c(0,1,1,1,-3,0,0,0,0)
c63a<-c(0,1,-1,0,0,0,0,0,0)
c64a<-c(0,1,1,-2,0,0,0,0,0)
linearHypothesis(m1a,rbind(c61a,c62a,c63a,c64a))
# Alternatively, if we reparameterize the model y_ij=mu_i +beta_j +e_ij
# so that the sum of the beta_j's is zero, then the estimated mu_i's
# will be the estimated treatment means averaged over the 4 blocks
# That is, they will be the lsmeans of for each treatment (as given by
# the LSMEANS statement in SAS or the lsmeans() function here in R)
options()$contrasts # the default contrasts for unordered factors is
# contr.treatment, so reset it to contr.sum
op <- options(contrasts=c("contr.sum", "contr.poly"))
options()$contrasts # now you can see that they've been changed
m1b <-lm(oil~treat+blockfac-1,data=redwing)
summary(m1b)
# now the following contrast will be in the lsmeans for trt. But as mentioned
# previously, in this no interaction model, we'll get the same answer as
# when we conducted the same contrast based on the trt means in block 1
# that is, when we did this contrast based on model m1a
linearHypothesis(m1b,rbind(c61a,c62a,c63a,c64a))
options(op) # reset contrasts to their original values