Mixed Model Example — Wagner et al. (2006)

[This article was first published on fishR » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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")


Filed under: Fisheries Science, R, Statistics Tagged: Mixed model, R, Statistics

To leave a comment for the author, please follow the link and comment on their blog: fishR » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)