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@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@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.
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]

