# redwing2.R
library(lsmeans)
library(nlme)
#library(lme4) # not used
#library(lmerTest) # not used
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)
redwing2 <- groupedData(oil~treat|blockfac, data=redwing)
# fit the linear model
m1<-lme(oil~treat,data=redwing2, random= ~1|blockfac)
summary(m1)
anova(m1,Terms="treat",type="marginal")
# Notice that the F test for treatments given above is different
# from the one given in redwing1.R when block effects are fixed.
# This contradicts what I said in class that the F test for treatments
# is unaffected by whether block effects are fixed or random. Actually,
# when Type 3 estimation of the mixed effect model is used, by statement
# is correct. For Type 3 estimation, the F test for treatments is the same
# regardless of whether block effects are fixed or random. However, the lme()
# function does not use Type 3 estimation. It uses REML estimation. REML
# and Type 3 estiamtion give the same results in the RCB model unless
# the Type 3 estimate of the block-to-block variance component is <0.
# If it is >=0, then reml and type 3 coincide and it remains true that
# treating block effects as fixed or random does not alter the F test on
# the treatments. However, in the unusual case that there is no
# evidence of block to block variability, then the type 3 estimate of the
# block variance component can be negative, whereas the REML estimator
# corrects this estimate and makes it 0. In that case (only) the F test for
# treatments given by REML estimation of the mixed effect model
# (block effects random) differs from that of the fixed effect model (block
# effects fixed). The truth is that this example (the redwing flaxseed
# example) is anomolous, it is unusual for there to be no block to block
# variance. In such a situation (and in general), you are really better off
# using REML estimation.
# Note that if you change the method of estimation used in PROC MIXED in
# redwing2.sas, you will see that it gives the same results as provided
# here in this program where we use lme to fit our model.
# Now get the lsmeans for treat. Note the SEs are different than when
# block effects are fixed (see redwing1.R).
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)))
# The lsmeans function call above gives test statistics that agree with those
# from SAS's PROC MIXED, but it treats these test statistics as Z tests, not
# t-tests. They are contrasts tests of the form C/s.e.(C) where C is the
# estimated contrast, but the lsmeans function does not implement these as
# t tests. Their squares are still the F statistics produced by the CONTRAST
# statement in SAS, but lsmeans in R treats these as Z tests and uses their
# large sample distribution, which is standard normal under the null
# hypothesis, rather than t (or F if squared). This is emphasized by the
# fact that the denominator DF given by lsmeans is listed as NA.
# Because the lsmeans function refers the statistics to their large sample
# std normal distribution rather than to a t distribution, it gives
# different p-values than does PROC MIXED.
# Similarly, if we use glht to do the 4 df test on the inoculation timing
# contrasts, it uses a large sample chisquare test rather than an F test.
# The chisq test statistic it uses is numDF*F where numDF=4 for this
# hypothesis (there are four contrasts being tested simultaneously) and
# F is the usual F statistic for the contrast (e.g., the F statistic produced
# by SAS' PROC MIXED.
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())
# The anova function can be used to test hypotheses on the treatment means
# with F tests. It is easier, though, if we first refit the model without
# an intercept (that is, in the parameterization y_ij=mu_i+b_j+e_ij)
# Then the fitted mu_i's are the estimated treatment means and specifying
# the contrasts become a but easier:
# refit with alternate parameterization:
m1a <-lme(oil~treat-1,data=redwing2, random= ~1|blockfac)
summary(m1a)
# Test inoculation timing:
anova(m1a,type="marginal", L=K1)
# Test average of inoculation treatment means vs control mean:
anova(m1a,type="marginal", L=c1)
# Again, note that the tests given above agree with PROC MIXED when
# method = REML is used. Given the negative variance component for blocks
# here, REML is the better metod. If blocks had not had a negative variance
# component, the REML and TYPE 3 results (including the above F tests)
# would have agreed.