[R-sig-ME] Constraining error variance to 0 in a mixed model.
Mollie Brooks
mollieebrook@ @ending from gm@il@com
Fri Jan 11 13:40:25 CET 2019
Hi Rolf,
Thanks for the enthusiasm. I believe the ability to force the variance into the RE in glmmTMB was Kasper Kristensen’s idea.
I had a look at your example and it seems that it could be a false positive coming from nlminb. The maximum gradient component is 2.4e-6, so it looks to me like the model converged.
In the future, the glmmTMB developers will reconsider the convergence check that relies on nlminb’s convergence codes.
if (!is.null(fit$convergence) && fit$convergence != 0) warning()
cheers,
Mollie
> On 11Jan 2019, at 8:53, D. Rizopoulos <d.rizopoulos using erasmusmc.nl> 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)
>
> Best,
> Dimitris
>
>
> -----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.
>
> --
> Honorary Research Fellow
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list