lme4: constrain sigma to 0

 Dear Dr. Bolker,

I am using lmer to estimate fixed effects for exerimental data (in
Psychology). At the moment I am trying to run a model in which the
residual variance is fixed to 0 and I am wondering whether it is
possible to run such a model with lmer. I found messages from 2011 (a
discussion between you and Rolf Turner) according to which it was not
possible to constrain sigma to 0. Now I am asking myself whether it is
now with a new version of the lme4 package possible.

Thanks a lot in advance and
best regards,
Maya Machunsky

No, still not possible (or at least not easy), because the
residual variance is profiled out, so it doesn't enter directly
into the calculations -- it's not a parameter you can set.

  You can sort of achieve this using blmer to set a 'point' prior
on sigma (i.e. fix it at a specific value) -- but you can't set the
residual variance to *exactly* zero.

library(lme4)
## fit model
fm0 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
library(blme)
fm1 <- blmer(Reaction ~ Days + (1 | Subject), sleepstudy,
       cov.prior=NULL,
       resid.prior=point(0.0001))

sumfun <- function(fm) {
    c(sd_Subject=unname(sqrt(unlist(VarCorr(fm)))),
      sd_Residual=sigma(fm),
      REMLdev=unname(deviance(fm)))

}
cfun <- function(sigma=1) {
    fm <- blmer(Reaction ~ Days + (1 | Subject), sleepstudy,
                cov.prior=NULL,
                resid.prior=point(sigma))
    sumfun(fm)
}
svec <- c(10^seq(-5,2,by=0.25),20:60)
sigmaprof <- rbind(sumfun(fm0),t(sapply(svec,cfun)))

sigmaprof <- as.data.frame(sigmaprof[order(sigmaprof[,"sd_Residual"]),])
par(las=1,bty="l")

## REML deviance vs. residual value -- LOG scale
plot(REMLdev~sd_Residual,data=sigmaprof,type="b",log="xy")
abline(v=sigma(fm0),col=2)
## leave out last point to improve scaling
plot(sd_Subject~sd_Residual,data=sigmaprof[-30,],type="b",log="x")
abline(v=sigma(fm0),col=2)

plot(REMLdev~sd_Residual,data=sigmaprof,
    subset=REMLdev-min(REMLdev)<30,type="b",log="xy")

Compare with the likelihood profile:

pp <- profile(fm1)
xyplot(pp)


There *might* also be a way to do this by digging into the
profiling machinery a bit more.

  As always, corrections and improvements are welcome.