# bowling.R
library(lsmeans)
library(lme4)
library(pbkrtest)
library(car)
# get the data
bowling<-read.table(file="bowling.dat",header=T,
colClasses=c("factor","factor","factor","factor","numeric"))
head(bowling)
is.factor(bowling$day)
is.factor(bowling$score)
is.ordered(bowling$treat)
# mixed effects model with random effects of day, time, and coach
# use package lme4
# For an additive model like this one, whether we treat day, time and coach
# effects as fixed or random does not affect the F test for treatments
# when using the Type 3 ANOVA fitting method in PROC MIXED/PROC GLM
# in SAS. Furthermore, the Type 3 fitting method gives the same results for
# this model as the REML fitting method, AS LONG AS THE VARIANCE COMPONENT
# ESTIMATES ARE POSITIVE! The equivalence of type 3 and reml inference
# breaks down when type 3 variance components estimates are negative, because
# REML doesn't allow such estimates and rounds them up to 0. This affects
# the F tests on fixed effects. In this particular example, the time and day
# variance components have negative type 3 anova estimates. Therefore, REML
# and type 3 fitting and inference in this model differ. In the presence of
# negative variance component estimates, there is some controversy about
# how best to test hypotheses on the fixed effects, but there is evidence
# that the inferences (F tests for main effects, contrasts) produced by
# the type 3 anova fitting method are better than those from REML.
# Therefore, in bowling.sas we used method=type3 in PROC MIXED to do the
# analysis, which would have agreed with the analysis provided by PROC GLM
# (assuming fixed day, time and coach effects), but, because of the negative
# variance components, would not agree with what we would have obtained with
# the default method=reml results.
# In R, the functions that exist for fitting linear mixed effect models
# (lme from the nlme package and lmer from the lme4 package) do not implement
# the type 3 anova method for fitting a linear mixed effect model. They only
# implement the REML and ML (maximum likelihood) methods. Therefore, it is
# not possible to reproduce the results from the call to PROC MIXED that we
# used in bowling.sas with the lme or lmer functions in R. The only way it
# can be done is to treat the model as a purely fixed effect model and fit
# it with the lm function.
# Below we first fit the model as a mixed effect model treating coach, day
# and time as random effects using reml in the function lmer(). The type 3
# (or marginal) F test for treat can then be obtained using the anova
# function if we specify the option type=3. The results that we get are
# based on reml, and would agree with SAS if we had changed method=type3
# to method=reml in our call to PROC MIXED in bowling.sas.
m1<-lmer(score~treat+(1|day)+(1|time)+(1|coach),data=bowling)
summary(m1)
AIC(m1)
BIC(m1)
anova(m1,type=3)
# test the contrasts of interest
c1<-c(1,-1,1,-1)
c2<-c(1,1,-1,-1)
c3<-c(1,-1,-1,1)
lsmeans(m1,specs=lsm~treat,contr=list(lsm=list(motor=c1,
visual=c2,
motor_visual=c3)))
# Alternatively, we fit the model treating day, coach and time effects as
# fixed using lm.
# 1st, change the handling of unordered factors to use sum-to-zero constraints
op<-options(contrasts=c("contr.sum","contr.poly"))
options()$contrasts
m1a <- aov(score~treat+day+time+coach,data=bowling)
# and check the type III SSs (they have changed and now agree with SAS):
Anova(m1a,type=3)
# note that standard errors and confidence intervals on the treatment
# means given below will be wrong because they treat coach, day, time
# effects as fixed. But the point estimates and F tests on contrasts
# are correct. For inference on the treatment means (e.g., to get appropriate
# SEs, confidence intervals) use the reml results based on model m1 given
# above.
lsmeans(m1a,specs=lsm~treat,contr=list(lsm=list(motor=c1,
visual=c2,
motor_visual=c3)))
options(op) # reset contrasts to their original values