[R] low-variance warning in lmer

[Ben Bolker]

  For block effects with small variance, lmer will sometimes
estimate the variance as being very close to zero and issue
a warning.  I don't have a problem with this -- I've explored
things a bit with some simulations (see below) and conclude that
this is probably inevitable when trying to incorporate
random effects with not very much data (the means and medians
of estimates are plausibly close to the nominal values:
the fraction of runs with warnings/near-zero estimates drops
from about 50% when the between-block variance is 1% of
the total (with 2 treatments, 12 blocks nested within treatment,
3 replicates per block), to 15% when between=30% of total,
to near zero when between=50% of total.

  My question is: what should I suggest to students when this
situation comes up?  Can anyone point me to appropriate references?
(I haven't found anything relevant in Pinheiro and Bates, but
I may not have looked in the right place ...)

  thanks
    Ben Bolker


self-contained but unnecessarily complicated simulation
code/demonstration:
---------------
  library(lme4)
library(lattice)

simfun <- function(reefeff,ntreat=2,nreef=12,
                   nreefpertreat=3,
                   t.eff=10,
                   totvar=25,seed=NA) {
  if (!is.na(seed)) set.seed(seed)
  ntot = nreef*nreefpertreat
  npertreat=ntot/ntreat
  reef = gl(nreef,nreefpertreat)
  treat = gl(ntreat,npertreat)
  r.sd = sqrt(totvar*reefeff)
  e.sd = sqrt(totvar*(1-reefeff))
  y.det = ifelse(treat==1,0,t.eff)
  r.vals = rnorm(nreef,sd=r.sd)
  e.vals = rnorm(ntot,sd=e.sd)
  y <- y.det+r.vals[as.numeric(reef)]+e.vals
  data.frame(y,treat,reef)
}

getranvar <- function(x) {
  vc <- VarCorr(x)
  c(diag(vc[[1]]),attr(vc,"sc")^2)
}

estfun <- function(reefeff,...) {
  x <- simfun(reefeff=reefeff,...)
  ow <- options(warn=2)
  f1 <- try(lmer(y~treat+(1|reef),data=x))
  w <- (class(f1)=="try-error" && length(grep("effectively zero",f1))>0)
  options(ow)
  f2 <- lmer(y~treat+(1|reef),data=x)
  c(getranvar(f2),as.numeric(w))
}


rvec <- rep(c(0.01,0.05,0.1,0.15,0.2,0.3,0.5),each=100)
X <- t(sapply(rvec,estfun))
colnames(X) <- c("reefvar","resvar","warn")
rfrac <- X[,"reefvar"]/(X[,"reefvar"]+X[,"resvar"])
fracwarn <- tapply(X[,"warn"],rvec,mean)
est.mean <- tapply(rfrac,rvec,mean)


op <- par(mfrow=c(1,2))
plot(rvec,rfrac,type="n",xlim=c(-0.02,0.55),axes=FALSE)
axis(side=2)
box()
boxplot(rfrac~rvec,at=unique(rvec),add=TRUE,pars=list(boxwex=0.03),
        col="gray")
points(jitter(rvec),rfrac,col=X[,"warn"]+1)
lines(unique(rvec),fracwarn,col="blue",type="b",lwd=2)
text(0.4,0.1,"fraction\nzero",col="blue")
abline(a=0,b=1,lwd=2)
points(unique(rvec),est.mean,cex=1.5,col=5)
##
plot(rvec,rfrac,type="n",xlim=c(-0.02,0.55),axes=FALSE,log="y")
axis(side=2)
axis(side=1,at=unique(rvec))
box()
points(jitter(rvec),rfrac,col=X[,"warn"]+1)
curve(1*x,add=TRUE,lwd=2)

print(densityplot(~rfrac,groups=rvec,from=0,to=0.9,
            panel=function(...) {
              panel.densityplot(...)
              cols <- trellis.par.get("superpose.line")$col
              lpoints(unique(rvec),rep(8,length(rvec)),type="h",
                      col=cols[1:length(rvec)])
            }))


-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 254 bytes
Desc: OpenPGP digital signature
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20061124/5dbf7a26/attachment.bin