Significance of Repeatablity Estimates

[cc'd to r-sig-mixed-models]

Hi Ben,

I was reading through Nakagawa and Schielzeth (2010) about how to
calculate the significance of repeatability estimates and they mention
your 2009 paper.  I can't seem to find this paper and was wondering if
you could clear up the procedure for me.  They state that you compare
models with and without a random intercept parameter, which would mean
comparing a linear model with a linear mixed-effect model, which can't
be done in R and returns 'Error: $ operator not defined for this S4
class'.  I also thought this wasn't a good idea to do in general.  Can
you briefly describe how this can be done in R?  Thank you for the help.

Jacob

The 2009 paper is available from
http://ms.mcmaster.ca/~bolker/bb-pubs.html (username: bbpapers,
password: research).

There are a variety ways of doing this:

library(lme4)

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy, REML=FALSE)
fm0 <- lm(Reaction ~ Days, sleepstudy)

## likelihood ratio test -- subject to the various caveats listed
##  at http://glmm.wikidot.com/faq#random-sig
## The log-likelihoods reported by lm() and lmer()
##  *are* commensurate; see test below ...

ddiff <- -2*(logLik(fm0)-logLik(fm1))
pchisq(ddiff,3,lower.tail=FALSE)

## use RLRsim (redo model with only a single random effect -- RLRsim
##  is limited in this way
fm2 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy, REML=FALSE)
library(RLRsim)
exactLRT(fm2,fm0)

## check commensurateness
dd <- update(fm2,devFunOnly=TRUE)
all.equal(dd(0),c(-2*logLik(fm0)))  ## TRUE

## parametric bootstrapping
bootSim <- function() {
   s <- simulate(fm0)[[1]]
   fm1B <- refit(fm1,s)
   fm0B <- update(fm0,data=transform(sleepstudy,Reaction=s))
   -2*(logLik(fm0B)-logLik(fm1B))
}
set.seed(101)
rr <- replicate(500,bootSim())
hist(rr,breaks=50,col="gray")
mean(ddiff<=c(rr,ddiff))
##  =1/501; essentially limited by size of bootstrap set