lme4: constrain sigma to 0

Several millennia ago I put in a feature-request that the facility for constraining sigma to 0 be added.  So far, as far as I can tell, it hasn't been.

Apparently it is (very) difficult to add such a constraint due to the way that lmer approaches the maximization of the likelihood.  (This, rather than any intrinsic mathematical block to such a constraint.)

I wanted to be able to fit a fairly simple repeated measures model with the covariance over time being an "arbitrary" n-x-n positive definite matrix (where "n" is the number of time points).  This can't be done
using lmer() unless one can constrain the "overall" variance to be zero,
otherwise the diagonal of the covariance matrix is un-identifiable.

Since one cannot constrain sigma to be 0, one can't fit this model with lmer().  Bummer.

cheers,

Rolf Turner

On 12/02/14 11:25, Asaf Weinstein wrote:

Dear lme4 community (and Ben),

I am interested in obtaining estimates for the usual linear model but with
a known/ fixed value of residual variance.
When searching for an answer in the forum I came across this question by
Maya (machunsky at uni-mannheim.de), and it looks like there still isn't a
built-in method to obtain MLEs under constraints. Is that correct?

Thank you so much,
Asaf


On 21 November 2013 14:48, Ben Bolker <bbolker at gmail.com> wrote:

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.

	[[alternative HTML version deleted]]