[R-sig-ME] glmmTMB predicted values and AIC, vs glmer
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Sat May 22 02:16:15 CEST 2021
What allFit indicates to me is that all of the different optimizers
reach very similar log-likelihoods (max difference of 1-4), and if I
recall correctly all of the fixed effects are also very similar, so I
would feel free to disregard the convergence warning as a false positive.
The AIC values are still a lot more different than I would like, and
digging into the gory details (available on request) it seems that they
are actually giving different likelihood calculations -- e.g. for the
parameters reached by glmer, glmmTMB returns a deviance of 607.037 while
glmer returns 609.3362. I suspect something deep in the guts of e.g.
the way the Laplace approximation is computed (not technically a bug in
either code, but some kind of numerical instability).
This is good news if what you're interested in is making sure you get
the correct parameter values. It's bad news if you're interested in
mixing and matching glmmTMB and lme4 fits in a model comparison.
I'm not horribly surprised that things are a bit numerically unstable
- 12 parameters for 95 observations, and the extremely rapid tail
behaviour of cloglog often makes things worse.
I do intend to track down where the infinite prediction is coming
from too, but haven't gotten there yet.
I'm going to include some of my code in case it's useful, but it's
poorly commented, messy *and may be in the wrong order* - I was bouncing
around developing it interactively.
====
Obsglmer1.cll <- glmer(form1, nAGQ=1,
family=binomial(link="cloglog"), data=ffly)
aa <- allFit(Obsglmer1.cll)
summary(aa)$fixef
ll <- -summary(aa)$llik
ll-min(ll)
dotwhisker::dwplot(modList[-2])
library(dotwhisker)
modList <- list(glmmTMB=ObsTMB.cll, nAGQ0=Obsglmer.cll,
nAGQ1=Obsglmer1.cll)
f1 <- ObsTMB.cll$obj$fn
## TMB params, TMB format
p1 <- with(environment(f1),last.par.best[-random])
## TMB params, glmer format
p1B <- c(exp(p1[names(p1)=="theta"]), p1[names(p1)=="beta"])
f3 <- getME(Obsglmer1.cll,"devfun")
## glmer params, glmer format
p3 <- unlist(getME(Obsglmer1.cll,c("theta","beta")))
## glmer params, TMB format
p3B <- c(p3[-(1:2)],log(p3[1:2]))
f3(p1B) ## 609.3559
f3(p3) ## 609.3362
2*f1(p3B) ## 608.037
2*f1(p1) ## 608.0166
Obsglmer2.cll <- update(Obsglmer1.cll,
start=list(theta=p1B[1:2],
fixef=p1B[-(1:2)]))
dwplot(modList)
library(bbmle)
AICtab(modList, logLik=TRUE)
modList2 <- c(modList, list(nAGQ1_newstart=Obsglmer2.cll))
AICtab(modList2, logLik=TRUE)
On 5/21/21 6:24 PM, John Maindonald wrote:
> Setting nAGQ=1 does bring the AIC values closer together, albeit with
> convergence failure when other parameters are left at their defaults.
> Running allFit() does not provide (to me) any obvious clues on what
> might work better. Or, should one just accept the less strict tolerance?
> Nor does centering and scale scTime help in this case.
>
> llik <- summary(allFit(Obsglmer1.cll))$llik
> max(lik)-lik
> ## bobyqa Nelder_Mead
> ## 0.00000000000e+00 1.31201379190e-04
> ## nlminbwrap nmkbw
> ## 1.76015646502e-09 3.48215053236e-06
> ## optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
> ## 1.09246872171e-06 1.13990722639e-07
> ## nloptwrap.NLOPT_LN_BOBYQA
> ## 5.81495783081e-07
>
> John Maindonaldemail: john.maindonald using anu.edu.au
> <mailto:john.maindonald using anu.edu.au>
>
>
>> On 22/05/2021, at 01:27, Ben Bolker <bbolker using gmail.com
>> <mailto:bbolker using gmail.com>> wrote:
>>
>> I don't know yet, I will dig in and try to see what's going on. An
>> infinite predicted value certainly seems like an issue.
>>
>> The most obvious difference is that nAGQ=0 is actually doing
>> something slightly different (it's fitting the fixed-effect parameters
>> as part of the PIRLS loop rather than maximizing over them
>> explicitly); I would rather compare the nAGQ=1 case, just to minimize
>> the number of differences.
>>
>> On 5/21/21 6:30 AM, John Maindonald wrote:
>>> The code that follows below demonstrates what I find a very odd
>>> issue with r-sig-mixed-models using r-project.org
>>> <mailto:r-sig-mixed-models using r-project.org><mailto:r-sig-mixed-models using r-project.org
>>> <mailto:r-sig-mixed-models using r-project.org>>, and with the comparison
>>> between
>>> glmmTMB::glmmTMB and lme4::glmer, for models that have the
>>> same model formulae.
>>> 1) In the glmmTMB model as fitted to `ffly`, one predicted value is
>>> infinite. This, even though all coefficients and standard errors
>>> have finite values. This surely indicates that there is something
>>> odd in the way that predictions on the scale of the linear
>>> predictor are calculated.
>>> 2) Initially, I'd thought that the `Inf` value explained the
>>> difference between the two AIC values. But the difference
>>> pretty much stays the same when the "offending" point is
>>> removed.
>>> I've not checked what happens when there are no observation level
>>> random effects. Is there some difference in the way that the
>>> model formulae are interpreted in the two cases?
>>> ffly <-
>>> read.csv('https://github.com/jhmaindonald/dataR/files/6521483/ffly.csv <https://github.com/jhmaindonald/dataR/files/6521483/ffly.csv>')
>>> ffly$obs <- factor(ffly$obs)
>>> form1 <- cbind(Dead,Live)~0+trtGp/TrtTime+(1|obs)+(1|trtGpRep)
>>> library(lme4); library(glmmTMB)
>>> ObsTMB.cll <- glmmTMB(form1,
>>> family=binomial(link="cloglog"), data=ffly)
>>> Obsglmer.cll <- glmer(form1, nAGQ=0,
>>> family=binomial(link="cloglog"), data=ffly)
>>> round(AIC(Obsglmer.cll, ObsTMB.cll), 2)
>>> ## df AIC
>>> ## Obsglmer.cll 14 639.44
>>> ## ObsTMB.cll 14 636.02
>>> round(c(max(predict(ObsTMB.cll)),max(predict(Obsglmer.cll))),3)
>>> ## [1] Inf 3.833
>>> range(fixef(ObsTMB.cll))
>>> ## [1] -5.28780064553 1.15952249272
>>> range(vcov(ObsTMB.cll))
>>> ## [1] -0.0546986905338 0.2887049942215
>>> ## Try also; there are small deviations from the line
>>> plot(predict(Obsglmer.cll)~predict(ObsTMB.cll))
>>> abline(0,1)
>>> ## Remove offending point
>>> ffly1 <- ffly[-which.max(predict(ObsTMB.cll)),]
>>> ObsTMB1.cll <- glmmTMB(form1,
>>> family=binomial(link="cloglog"), data=ffly1)
>>> Obsglmer1.cll <- glmer(form1, nAGQ=0,
>>> family=binomial(link="cloglog"), data=ffly1)
>>> cbind(AIC(Obsglmer.cll, ObsTMB.cll), AIC1=AIC(Obsglmer1.cll,
>>> ObsTMB1.cll)[,2])
>>> ## df AIC AIC1
>>> ## Obsglmer.cll 14 639.441888969 639.441889196
>>> ## ObsTMB.cll 14 636.016597730 636.016597723
>>> ## Observe that changes are in the final four decimal places.
>>> round(rbind(glmer=fixef(Obsglmer1.cll),
>>> TMB=unclass(fixef(ObsTMB1.cll)$cond)),3)
>>> trtGpspAEgg trtGpspAL2 trtGpspAL3 trtGpspBEgg trtGpspBL2 trtGpspBL3
>>> glmer 0.772 -2.255 -3.299 -1.803 -2.632 -5.114
>>> TMB 0.790 -2.367 -3.441 -1.889 -2.757 -5.288
>>> trtGpspAEgg:TrtTime trtGpspAL2:TrtTime trtGpspAL3:TrtTime
>>> glmer 0.231 0.372 0.563
>>> TMB 0.288 0.398 0.602
>>> trtGpspBEgg:TrtTime trtGpspBL2:TrtTime trtGpspBL3:TrtTime
>>> glmer 0.278 0.517 1.112
>>> TMB 0.299 0.554 1.160
>>> John Maindonald email: john.maindonald using anu.edu.au
>>> <mailto:john.maindonald using anu.edu.au><mailto:john.maindonald using anu.edu.au
>>> <mailto:john.maindonald using anu.edu.au>>
>>> [[alternative HTML version deleted]]
>>> _______________________________________________
>>> R-sig-mixed-models using r-project.org
>>> <mailto:R-sig-mixed-models using r-project.org> mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>>>
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org
>> <mailto:R-sig-mixed-models using r-project.org> mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>
More information about the R-sig-mixed-models
mailing list