On 11/24/06, Ben Bolker <[EMAIL PROTECTED]> wrote: > For block effects with small variance, lmer will sometimes > estimate the variance as being very close to zero and issue > a warning.
The "very close to zero" is merely a computational convenience. The theoretical lower bound on the variance components is zero but evaluating the profiled log-likelihood or log-restricted likelihood when the variance-covariance matrix of the random effects is singular would require a completely different algorithm. Evaluating the gradient and the ECME update in this situation would be even more complicated. Instead of doing that I set the lower bounds to a value close to zero but still large enough that the evaluation algorithm used elsewhere works here too. The value I use for the relative variances is 5e-10. The corresponding estimate of the variance component is 5e-10 * s^2. So the algorithm hasn't really converged to 5e-10 * s^2. The estimate of the variance component should be zero but I haven't allowed it to go all the way to zero. > I don't have a problem with this I hope not. The maximum likelihood or REML estimates of variance components can legitimately be zero. > -- 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? If it is a variance component that is estimated as zero then just drop that term from the model. Having the MLE or the REML estimate of a variance component be zero is pretty strong evidence that it is not significantly greater than zero (in the usual sense of failing to reject the null hypothesis - it still could be the case that it is greater than zero but we don't have strong enough evidence to reject the case that it is zero). Hence we go with the simpler model that eliminates this random effect. The more difficult case is when the estimate of the variance-covariance matrix of vector-valued random effects is singular. Say you have a random effect for both the intercept and the slope w.r.t. time for each subject. It is possible to converge to a variance-covariance matrix for these random effects with a non-zero variance for the intercept and a non-zero variance for the slope but perfect correlation between these random effects. I'm not sure what to suggest in that case as the reduced model is, as far as I can see, not in the class of linear mixed models. I think I might just not bother mentioning this in a classroom. > 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)]) > })) > > > > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > > > > ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.