[R-sig-ME] lme4: constrain sigma to 0

[Ben Bolker]

On 13-11-21 12:42 PM, Maya Machunsky wrote:
>  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.