[R-sig-ME] glmmTMB predicted values and AIC, vs glmer

[Ben Bolker]

   I figured out the infinite-prediction thing (didn't solve it, but 
figured it out):

https://github.com/glmmTMB/glmmTMB/issues/696

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>
>