Mixed Model Example — Wagner et al. (2006)

I am preparing for a workshop on mixed models and looked at the paper “Accounting for multilevel data structures in fisheries data using mixed models” by Wagner et al. (2006) (PDF available here).  Wagner et al. (2006) used two examples, with the first being simulated fish density data from eight streams, four of which had a barrier and four of which did not, and each stream was sampled above and below a barrier (i.e., a reference location in the non-barrier stream).  The authors analyzed the data (incorrectly, but as an illustration) as a traditional linear model without acknowledging nesting within the streams and as a mixed model that acknowledged this dependency.  They also noted, but did not include, that another traditional analysis would be to average among sites in the stream and analyze the averages.  My recreation of their example (including the mentioned but included incorrect analysis), with some comments, is shown below.

#=================================================================
# Create the data
#=================================================================
dens <- c(15,14, 7,38,13,24,
          49,47,39,45,50,50,
          14,15,25,21,31,33,
          12, 9,18,20,24,30,
          46,41,42,49,50,55,
          67,70,82,76,76,74,
          45,45,57,39,37,35,
          74,67,72,66,64,72)
stream <- rep(c("A","B","C","D","E","F","G","H"),each=6)
treat <- rep(c("barrier","control"),each=6*4)
pos <- rep(rep(c("above","below"),each=3),times=8)
site <- rep(1:6,times=8)
df <- data.frame(stream,treat,pos,site,dens)
str(df)

#=================================================================
# Load packages
#=================================================================
library(lsmeans)   # for lsmeans()
library(plyr)      # for ddply()
library(nlme)      # for lme()

#=================================================================
# WRONG ANALYSIS -- Disaggregated analysis (ignores stream)
#   used type-I SS in anova() as they matched the paper (which
#     said type-III)
#   lsmeans match paper, Tukey results match paper
#=================================================================
lm1 <- lm(dens~treat*pos,data=df)
anova(lm1)
lsmeans(lm1,pairwise~pos:treat,adjust="tukey")

#=================================================================
# WRONG ANALYSIS -- Aggregated analysis (avg within a stream)
#   NOT in the paper
#=================================================================
# First, create a data.frame that contains the mean density at
#   each stream, barrier, position combination.
df2 <- ddply(df,.(stream,pos,treat),function(df) mean(df$dens))
names(df2)[4] <- "meanDens"
# Second, fit the two-way model to the means.  Note that the
#   lsmeans are the same in both methods.  However, the SE are
#   very different as are the multiple comparisons (above the
#   barrier is different than both control locations)
lm2 <- lm(meanDens~treat*pos,data=df2)
anova(lm2)
lsmeans(lm2,pairwise~pos:treat,adjust="tukey")

#=================================================================
# PROPER ANALYSIS - Mixed-Model ... random intercepts for streams
#   lsmeans results generally match paper, though a strict
#     interpretation of the Tukey's results would be slightly
#     different
#   Note the significant interaction term here.
#=================================================================
lme1 <- lme(dens~treat*pos,random=~1|stream,data=df)
anova(lme1)
summary(lme1)
lsmeans(lme1,pairwise~pos:treat,adjust="tukey")