# redwing3.R
library(doBy)
library(nlme)
library(lsmeans)
#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="redwing3.dat",header=T)
head(redwing)
is.factor(redwing$block)
redwing$blockfac <- factor(redwing$block)
# Method 1: compute mean response for each plot and analyze those means
# with standard model for an unreplicated RCBD
# to get the mean for each plot we can use the
# summaryBy function in the doBy package
(EUdata <- summaryBy(oil ~ treat + blockfac, data = redwing,
FUN = function(x) { c(m = mean(x)) } ) )
# produces oil.m for each
# combination of the levels of treat and blockfac
redwing2 <- groupedData(oil.m~treat|blockfac, data=EUdata)
# fit the linear mixed model
m1<-lme(oil.m~treat,data=redwing2, random= ~1|blockfac)
summary(m1)
anova(m1,Terms="treat",type="marginal")
# test the contrasts of interest
# Here we use the lsmeans() function in package lsmeans
c1<-c(1,1,1,1,1,-5)
lsmeans(m1,specs=lsm~treat,contr=list(lsm=list(treat.v.control=c1)))
# The lsmeans function call above gives a test statistic that agrees with the
# F stats from 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.
# 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 bit easier:
# refit with alternate parameterization:
m1a <-lme(oil.m~treat-1,data=redwing2, random= ~1|blockfac)
summary(m1a)
# Test average of inoculation treatment means vs control mean:
anova(m1a,type="marginal", L=c1) # agrees with SAS's result from CONTRAST
# Method 2: analyze data at subsample level.
redwing$EU <- factor(paste(redwing$treat,redwing$block, sep=""))
head(redwing)
redwing3 <- groupedData(oil~treat|blockfac/EU, data=redwing)
# fit the linear mixed model
m2 <- lme(oil~treat,data=redwing3, random=list(blockfac= ~1, EU= ~1))
summary(m2)
# Now test main effect of treat
anova(m2,Terms="treat",type="marginal")
# test the contrasts of interest
# Here we use the lsmeans() function in package lsmeans
lsmeans(m2,specs=lsm~treat,contr=list(lsm=list(treat.v.control=c1)))
# Again, lsmeans uses a large sample Z test rather than an F test.
# The F test can be obtained in the same way we did previously: refit
# the model without an intercept and use the anova() function as follows:
# refit with alternate parameterization:
m2a <- lme(oil~treat-1,data=redwing3, random=list(blockfac= ~1, EU= ~1))
summary(m2a)
# Test average of inoculation treatment means vs control mean:
anova(m2a,type="marginal", L=c1) # agrees with SAS's result from CONTRAST