# heifers.R
library(lsmeans)
library(lme4)
#library(pbkrtest)
library(car)
# get the data
heifers<-read.table(file="heifers.dat",header=T,
colClasses=c("factor","factor","factor","factor","numeric"))
head(heifers)
is.factor(heifers$farm)
is.factor(heifers$consumpt)
# to fit this model with the lmer function in the lme4 package it is better
# to avoid implicitly nested factors. Instead create a heifer factor (called
# allheif below) that has 3*2=6 levels (a unique level for each heifer within
# each farm).
heifers <- within(heifers , allheif <- factor(heifer:farm))
head(heifers)
# mixed effects model with random effects of week heifer(farm) farm
# using reml in the function lmer() from package lme4
m1<-lmer(consumpt~feed+(1|week)+(1|heifer)+(1|farm),data=heifers)
summary(m1)
anova(m1,type=3)
# get the lsmeans for each level of feed
lsmeans(m1,specs=~feed)
# Notice that the results above, including the type 3 F test for feed,
# estimates of variance components, and standard errors of the estimated
# means for each feed type, do not agree with the results from PROC MIXED
# in heifers.sas. The reason for this is that to fit this model in heifers.sas
# I used method=type3, whereas the lmer function uses REML estimation.
# REML and Type3 estimation typically agree for this model, but that agreement
# breaks down when the type 3 variance component estimates are negative. In
# that case, REML constrains the corresponding variance component estimates
# to be non-negative, and that affects other results like the F tests on
# the fixed effects and standard errors on the estimated means. Type 3
# estimation of a mixed effect model is not implemented in lmer or any of
# the other functions designed for linear mixed models that are available in
# R. However, the model can be fit with lm with heifers, farms and weeks
# treated as fixed effects. That approach will give the an F test for feeds
# that agrees with the type 3 result from PROC MIXED that we obtained in
# heifers.sas. However, other results form that analysis, including standard
# errors and confidence intervals on estimated means, will not be correct.
# to get the same Type 3 Analysis of Variance as in SAS
# we need to fit the model treating week heifer farm as fixed effects using lm
# first, change the handling of unordered factors to use the sum-to-zero constraints
op<-options(contrasts=c("contr.sum","contr.poly"))
options()$contrasts
m1a<-aov(consumpt~feed+week+farm/heifer+farm,data=heifers)
# get the type III SSs
Anova(m1a,type=3)
# Alternatively, if we are not interested in separating farm and
# heifer within farm variability, we can just account for 6 heifer effects
# in the model. In that case, the Type 3 F test for feed types produced
# by the Type 3 anova estimation method (as implemented in PROC MIXED using
# using method=type3) or in a a fixed-effect model fitting routine (like
# PROC GLM or the lm function in R) will not change.
m2a<-aov(consumpt~feed+week+allheif,data=heifers)
# get the type III SSs
Anova(m2a,type=3)
# reset contrasts to their original values
options(op)