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

Ben Bolker bbo|ker @end|ng |rom gm@||@com
Fri May 21 15:27:42 CEST 2021


   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>, 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')
> 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>
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



More information about the R-sig-mixed-models mailing list