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

[Rolf Turner]

I can write out the specifics without much difficulty.  I want to be 
careful about being absolutely clear in what I say, but, so it will take 
me a little while.  I will get back to you fairly soon.

Thanks.

cheers,

Rolf Turner

On 12/02/14 12:48, Vincent Dorie wrote:
> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>