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

[Rolf Turner]

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
>