[R-sig-ME] Constraining error variance to 0 in a mixed model.

[Rolf Turner]

On 1/11/19 8:53 PM, D. Rizopoulos wrote:

> For what it's worth, the same model can be fitted with nlme::gls without any warning messages:
> 
> library(nlme)
> fit4 <- gls(test ~ 0 + occasion, data = Dat,
>              weights = varIdent(form = ~ 1 | occasion),
>              correlation = corSymm(form = ~ 1 | person))
> 
> vcov(fit4)

Thanks.  That's interesting, but a bit obscure from the point of view of 
a naive young ( :-) ) thing like myself.  The syntax of the call to 
gls() is mind-boggling --- at least to my feeble mind.

What I really wanted (as I have remarked in similar posts in the past) 
is to be able to fit simple mixed models in a simple-minded way, to 
check that my call to the more sophisticated software is correct.  This 
desideratum kind of falls apart if the sophisticated software 
cannot/won't fit the simple model (at least not without throwing a 
hissy-fit).

I think that what I *can* take away from your fitting procedure is that 
there are/should be no insurmountable intrinsic numerical barriers to 
fitting the simple model that I propose by means of sophisticated software.

cheers,

Rolf

-- 
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276


> -----Original Message-----
> From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> On Behalf Of Rolf Turner
> Sent: Friday, January 11, 2019 3:30 AM
> To: Help Mixed Models <r-sig-mixed-models using r-project.org>
> Subject: [R-sig-ME] Constraining error variance to 0 in a mixed model.
> 
> 
> Let me start off by apologising for beating this issue to death on this mailing list.  I have raised it (without getting satisfactory answers) on several occasions.  If you are fed up to the teeth with my maundering s, please just press the delete button.
> 
> This latest occasion was triggered by my trying to learn to use the glmmTMB function (from the package of the same name).  I noticed (to my initial delight) that glmmTMB() has an argument "dispformula" and the help says (in particular) "In Gaussian mixed models, dispformula=~0 fixes the [dispersion] parameter to be 0, forcing variance into the random effects."
> 
> I said "Ah-ha!  The very thing."  So I tried it out, using the following data set (supplied by Ben Pelzer in a posting to this list, and previously used by my very good self to illustrate the problem that I am trying to solve):
> 
> Dat <- structure(list(person = c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 9L, 9L, 9L, 10L, 10L, 10L, 11L, 11L, 11L, 12L, 12L, 12L, 13L, 13L, 13L, 14L, 14L, 14L, 15L, 15L, 15L, 16L, 16L, 16L, 17L, 17L, 17L, 18L, 18L, 18L, 19L, 19L, 19L, 20L, 20L, 20L), occasion = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1", "2", "3"), class = "factor"), test = c(25L, 21L, 27L, 27L, 19L, 23L, 20L, 18L, 18L, 26L, 33L, 46L, 25L, 36L, 47L, 36L, 35L, 41L, 30L, 30L, 37L, 23L, 21L, 19L, 39L, 37L, 33L, 29L, 36L, 49L, 29L, 34L, 44L, 39L, 32L, 30L, 26L, 32L, 36L, 33L, 21L, 8L, 30L, 34L, 34L, 27L, 34L, 39L, 41L, 40L, 44L, 37L, 34L, 34L, 23L, 26L, 30L, 31L, 28L, 27L)), row.names = c(NA, -60L), class = "data.frame")
> 
> I tried three ways of "analysing" these data:
> 
> (1) Simple-minded multivariate analysis:
> 
> X   <- matrix(Dat$test,byrow=TRUE,ncol=3)
> colnames(X) <- paste0("occasion",1:3)
> mu  <- apply(X,2,mean)
> Sig <- var(X)/20 # X has 20 rows
> 
> (2) Using lmer():
> library(lme4)
> fit2 <- lmer(test ~ 0 + occasion + (0 + occasion | person),data=Dat,
>                control=lmerControl(check.nobs.vs.nRE = "ignore"))
> 
> (The "control" argument was gathered from a posting by Maarten Jung, on which I have previously commented.)
> 
> Note that this "works" but throws (as I have previously noted) a disconcerting warning message:
> 
>> Warning message:
>> In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>>    Model is nearly unidentifiable: large eigenvalue ratio
>>   - Rescale variables?
> 
> (3) Using glmmTMB():
> library(glmmTMB)
> fit3 <- glmmTMB(test ~ 0 + occasion + (0 + occasion | person),data=Dat,
>                   dispformula = ~0)
> 
> (No warning message; initially I said "Wheeeee!!!")
> 
> The "parameter estimates" mu, fixef(fit2) and fixef(fit3) all agree:
> 
> occasion1 occasion2 occasion3
>       29.80     30.05     33.30
> 
> The covariance Sig and vcov(fit2) agree modulo numerical noise:
> 
>   > Sig
>             occasion1 occasion2 occasion3
> occasion1 1.7821053  1.150526 0.5768421
> occasion2 1.1505263  2.249868 3.0623684
> occasion3 0.5768421  3.062368 5.8531579
> 
>   > vcov(fit2)
> 3 x 3 Matrix of class "dpoMatrix"
>             occasion1 occasion2 occasion3
> occasion1 1.7821054  1.150526 0.5768422
> occasion2 1.1505265  2.249868 3.0623683
> occasion3 0.5768422  3.062368 5.8531577
> 
> However vcov(fit3) differs; after some head scratching I realised that it is equal to 19/20 times the others.  After a little more head scratching I said ah-ha!  The glmmTMB() function sets REML=FALSE by default.  So I'll just call it with REML=TRUE:
> 
> fit3 <- glmmTMB(test ~ 0+occasion + (0+occasion | person),data=Dat,
>                   dispformula = ~0,REML=TRUE)
> 
> And now vcov(fit3) more-or-less agrees with the other results:
> 
>   > vcov(fit3)
> Conditional model:
>             occasion1 occasion2 occasion3
> occasion1 1.7819337  1.150276 0.5765386
> occasion2 1.1502758  2.249755 3.0623895
> occasion3 0.5765386  3.062390 5.8534164
> 
> But, to my extreme irritation, glmmTMB() now also throws a warning:
> 
>> Warning message:
>> In fitTMB(TMBStruc) :
>> Model convergence problem; false convergence (8). See
>> vignette('troubleshooting')
> 
> (I've looked at the vignette, and I can sort of see a connection, but not really.)
> 
> *Why* does the bloody universe always do this sort of thing to me?  Why on earth should dividing by 19 rather than 20, effectively, cause such a warning to be thrown?
> 
> Isn't something a little out of whack here?
> 
> I reiterate at this point a remark that I made in a previous post on the issue of constraining the error variance to 0:
> 
>> What we have here is a "boundary" or "edge" or "corner" case, and I
>> have heard it asserted by someone knowledgeable and respected (can't
>> remember who, it was a long while ago) that in testing your software
>> it is crucial to see how it behaves with such "boundary" cases.
> 
> Am I being unreasonable?  (*Me*?  *Unreasonable*?  Perish the thought!) :-)
> 
> cheers,
> 
> Rolf Turner
> 
> P.S.  Note that setting REML=FALSE in the call to lmer() causes it to produce a covariance matrix that agrees with that produced by glmmTMB() with default setting.  However this has no impact upon the error message issued by lmer().
> 
> R. T.
>