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

[Vincent Dorie]

Any chance someone could write out the specifics of such a model? If I can wrap my head around it and its not too hard, I could try to throw it into blme.

On Feb 11, 2014, at 6:18 PM, Rolf Turner wrote:

> 
> 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.
>>> 
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>> 
>> 
>> 	[[alternative HTML version deleted]]
>> 
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> 
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models