3
$\begingroup$

The context is about the use of a given model deviance (often referred to as “Residual deviance” in R) and that of its “Null deviance” to calculate D2, the deviance explained for models with non-normal error.

This metric (D2) is somehow an approximation of what a coefficient of determination (R2) is for a linear model, providing an approximate idea of the variation explained by your model, although it is probably better described as a value reflecting how close its fit is from being perfect when compared to a saturated model, see:

https://bookdown.org/egarpor/PM-UC3M/glm-deviance.html

See also p. 166-167 of Guisan & Zimmermann (2000) for the calculation of D2 and its related adjusted version:

https://www.wsl.ch/staff/niklaus.zimmermann/papers/ecomod135_147.pdf

With most packages, such as “stats” and its glm() function, “MASS” and glm.nb(), and “mgcv” and gam(), those two deviance-related values can easily be obtained with "deviance" and "null.deviance" (preceeded by a dollar symbol - see further below). The gam() function actually also provides D2 directly in its output.

For the glmmTMB package and its glmmTMB() function, I haven't been able to extract the residual deviance nor the null deviance to calculate D2.

Does anyone has an idea of how to achieve this?

Here's an example of what I'm looking for, using real count data from a sample of walleyes through a monitoring program. The idea is to extract the deviance-related information for a same model using each function/package when fitting these with the same error structure (NB2, link=log) and see how they compare, especially when contrasted to glmmTMB (not being able to extract such information from this package).

Walleye catch curve example:

From the "descending limb of a catch curve" in fisheries, one can model the rate at which counts decrease with age to estimate the instantaneous mortality (Z: the absolute value of the age coefficient) on the log scale from a sample of randomly-captured fish. The age-frequencies data are as follow:

age<-seq(1,15,by=1)
count<-c(151,56,117,10,12,21,8,2,2,1,2,0,1,1,2)
walleye<-data.frame(age,count)

These data are over-dispersed (variance > mean) and as such, the Poisson family distribution is inadequate. Using the glm.nb() function of the "MASS" package allows to model the variance in "extra" as follow:

summary(m.walleye.nb2<-glm.nb(count~age,data=walleye))

Call:
glm.nb(formula = count ~ age, data = walleye, init.theta = 3.114212171, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6522  -0.8761  -0.2121   0.5188   1.8561  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.21804    0.36336  14.361  < 2e-16 ***
age         -0.42792    0.05395  -7.931 2.17e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(3.1142) family taken to be 1)

    Null deviance: 109.841  on 14  degrees of freedom
Residual deviance:  15.903  on 13  degrees of freedom
AIC: 95.102

Number of Fisher Scoring iterations: 1


              Theta:  3.11 
          Std. Err.:  1.58 

 2 x log-likelihood:  -89.102

Note that both the Null deviance (109.841) and the Residual deviance (15.903) are provided directly in the output and already tell me that the predicted values of the considered model adjust pretty well to the observed (count) data given the large difference between the Null and Residual deviances.

These two deviance-related values can also be directly extracted as follow:

m.walleye.nb2$null.deviance
[1] 109.8413

m.walleye.nb2$deviance
[1] 15.90327

Using first the hnp() function of the hnp package (Moral et al. 2017) I can tell that the Pearson residuals of my model are behaving as expected given the distributional assumptions of this Poisson extension, i.e. the negative binomial type II (nb2). This model is thus adequate (i.e., goodness-of-fit), but I'm interested too in obtaining a "calibration" metric to get an idea of how well the predictions adjust to the observed data. That's where D2 is useful, and more so from an "explanatory power" than "adequacy" perspective:

D2<-100*(1-m.walleye.nb2$deviance/m.walleye.nb2$null.deviance)
D2
[1] 85.52159

This quite high value is not unexpected, as counts are deacreasing with age (walleyes die as they age) and despite the high variation observed in these age-frequencies data, the .nb2 model is capturing most of the signal for the central tendency and its associated variance.

If I use the gam() function of the mgcv package, I can run the exact same model when not recoursing to a smoothing function s() and specifying the argument method="ML" as REML would otherwise be used by default:

library(mgcv)
summary(m.walleye.nb2.GAM<-gam(
count~age,family=nb(theta=NULL),method="ML",data=walleye))

Family: Negative Binomial(3.114) 
Link function: log 

Formula:
count ~ age

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.21804    0.36336  14.361  < 2e-16 ***
age         -0.42792    0.05395  -7.931 2.17e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


R-sq.(adj) =  0.774   Deviance explained = 85.5%
-ML = 44.551  Scale est. = 1         n = 15

Note that the parameter estimates are identical to those obtained with glm.nb() of the "MASS" package. The Deviance explained provided in the output corresponds to the one that was calculated for the previous .nb2 model. If we calculate D2 by hand for this nb2.GAM model:

D2<-100*(1-m.walleye.nb2.GAM$deviance/m.walleye.nb2.GAM$null.deviance)
D2
[1] 85.52159

Identical.

Now, fitting the same model in glmmTMB using .NB2 instead of .nb2 to differentiate it from the one obtained with glm.nb(), we get:

library(glmmTMB)
summary(m.walleye.NB2<-glmmTMB(count~age,family=nbinom2,data=walleye))  
 Family: nbinom2  ( log )
Formula:          count ~ age
Data: walleye

     AIC      BIC   logLik deviance df.resid 
    95.1     97.2    -44.6     89.1       12 


Dispersion parameter for nbinom2 family (): 3.11 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.21804    0.34947  14.931  < 2e-16 ***
age         -0.42792    0.05485  -7.802 6.09e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The parameter estimates are again identical to those of .nb2 and .nb2.GAM and the deviance provided is of 89.1.

All models have the same log-likelihood:

logLik(m.walleye.nb2)
'log Lik.' -44.55119 (df=3)

logLik(m.walleye.nb2.GAM)
'log Lik.' -44.55119 (df=3)

logLik(m.walleye.NB2)
'log Lik.' -44.55119 (df=3)

As indicated by @Ben Bolker, the deviance reported in glmmTMB is:

deviance<-(-2*-44.55119)
deviance
[1] 89.10238

The "deviance" here is obviously the same for all three models (89.10238) and is used for instance in the calculation of the AIC = -2(log-likelihood) + 2K where K is the number of parameters (df=3 from above):

89.10238+(2*3)  
[1] 95.10238

In the end, the more specific question would be:

How can I extract the deviance-related information (null and residual deviances) from a glmmTMB object for the computation of D2 in such a simple model (one predictor, fixed effect)?

Although not of prime importance, this metric nonetheless helps to provide an idea of how "good" is your model at explaining the variation in your data, knowing that a top-ranking, adequate model may have a low explanatory power while still being useful. As such, knowing this information is desirable to sometimes tone down a statement related to the model predictions. I guess that Simon Wood has not provided D2 in the gam() function output of his "mgcv" package for no reason.

$\endgroup$

3 Answers 3

5
$\begingroup$

This would be a comment, but I need to provide some code, so I'll include it as an answer.

This may be in essence a programming question, but there might be a statistics question as well.

  • With glmmTMB objects, you may need to fit a null model to determine the null deviance. What constitutes a null model may depend on your purposes.

  • You should be able to derive the deviance from the log likelihood. The latter can be extracted with logLik(model).

  • The anova() function displays the deviance.

library(glmmTMB)
     
model = glmmTMB(count ~ mined + (1|site), family=poisson, data=Salamanders)
     
model.null = glmmTMB(count ~ 1+ (1|site), family=poisson, data=Salamanders)
      
logLik(model)
     
   ### 'log Lik.' -1104.849 (df=3)
     
logLik(model.null)
    
   ### 'log Lik.' -1120.773 (df=2)
    
anova(model, model.null)
    
   ###            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
   ### model.null  2 2245.6 2254.5 -1120.8   2241.6                             
   ### model       3 2215.7 2229.1 -1104.8   2209.7 31.847      1  1.668e-08 ***
$\endgroup$
2
  • $\begingroup$ Thanks for the feedback @Sal Mangiafico $\endgroup$ Apr 11, 2023 at 18:33
  • $\begingroup$ The deviances obtained with the anova() function in the example you have provided are the ones reported too in a glmmTMB output (see Ben Bolker's answer) and can be used for instance to calculate the AIC of the related model. Here, the AIC of "model.null" would be -2logLik + 2 K = -2*(-1120.8)+2*2 = 2245.6. But is there a way to calculate the Residual deviance from the deviance of "model" and the Null deviance from the deviance of "model.null"? Otherwise, if one could extract the deviance residuals of each models, that could resolve the issue. See my comment to Ben Bolker. $\endgroup$ Apr 12, 2023 at 13:15
4
$\begingroup$

Another general comment, from the details section of ?lme4::deviance.merMod:

Deviance and log-likelihood of GLMMs:

One must be careful when defining the deviance of a GLM. For example, should the deviance be defined as minus twice the log-likelihood or does it involve subtracting the deviance for a saturated model? To distinguish these two possibilities we refer to absolute deviance (minus twice the log-likelihood) and relative deviance (relative to a saturated model, e.g. Section 2.3.1 in McCullagh and Nelder 1989). With GLMMs however, there is an additional complication involving the distinction between the likelihood and the conditional likelihood. The latter is the likelihood obtained by conditioning on the estimates of the conditional modes of the spherical random effects coefficients, whereas the likelihood itself (i.e. the unconditional likelihood) involves integrating out these coefficients. The following table summarizes how to extract the various types of deviance for a ‘glmerMod’ object:

                     conditional        unconditional 
   relative   ‘deviance(object)’         NA in ‘lme4’ 
   absolute  ‘object@resp$aic()’  ‘-2*logLik(object)’ 
library(lme4)
library(glmmTMB)
m1 <- glmmTMB(count~spp * mined + (1|site), data = Salamanders, family="nbinom2")
m2 <- glmer.nb(count~spp * mined + (1|site), data = Salamanders)
  • -2*logLik() is the same (up to a small numerical difference) for glmmTMB and lme4 (approx 1631.3)
  • deviance(m2) is 501.8 (NULL for m1)
  • m2@resp$aic() is 1584.3 (undefined for m1)

For what it's worth, the canonical definition of deviance is what's called "relative" deviance above: people (including me) are sometimes sloppy and call $-2 \log {\cal L}$ (i.e. the "absolute deviance" from above) a deviance, but I believe that's technically incorrect. If we were only interested in differences between deviances rather than ratios this distinction wouldn't matter ...

Since you're not dealing with GLMMs (no random effect in your example) this gets a little bit less messy, but you still have a problem. We can fit the null model via update(my_model, . ~ 1), but that only gets us the values of $-2 \log{\cal L}$ ("absolute") for the full and null models, not the deviances.

$\endgroup$
2
  • $\begingroup$ Tangentially, I am currently having a brain fart: I can't see why the sum of the squared deviance residuals [which is how the (relative) deviance is usually computed] is actually different from $-2 \log{\cal L}$ ...? $\endgroup$
    – Ben Bolker
    Apr 11, 2023 at 14:50
  • $\begingroup$ Thanks for the feedback @Ben Bolker. I've read that we can indeed use the sum of the squared deviance residuals to estimate the residual deviance, but unless I'm wrong, glmmTMB only allows to extract either the “response”, “pearson”, or “working” residuals. If I use: sum(residuals(m.walleye.nb2,type="deviance")^2) for the model obtained with glm.nb() of the MASS package, I get the exact same value regarding the residual deviance: [1] 15.90327. The same hold true for the .nb2.GAM from gam() of the mgcv package. Identical (correct) value. Being able to do so with glmmTMB would offer a solution. $\endgroup$ Apr 11, 2023 at 18:52
1
$\begingroup$

Until a more optimal answer is provided, here's a solution based on the information that I've found in a "GLM blog" by Matthias Dörig:

The glmmTMB() function allows to extract the "response", "working" and "pearson" residuals of a considered model, but (to my knowledge) not the "deviance" residuals. This is unfortunate, as the "Residual deviance" of a GLM is simply the sum of the squared "deviance" residuals.

Using the same walleye data, I'll start by using the glm() function with family=poisson as the solution relies on transforming the "response" residuals into "deviance" residuals based on an error-specific function that was provided by Dr. Dörig for the Poisson case.

summary(m.walleye.p<-glm(count~age,family=poisson,data=walleye))

Call:
glm(formula = count ~ age, family = poisson, data = walleye)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.0251  -1.4546  -0.1779   0.8031   7.0199  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.45450    0.08303   65.70   <2e-16 ***
age         -0.47326    0.02485  -19.05   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 851.04  on 14  degrees of freedom
Residual deviance: 111.96  on 13  degrees of freedom
AIC: 169.06

Number of Fisher Scoring iterations: 5

When using glmmTMB() to model the same data in an identical manner, we aim to obtain a Residual deviance of 111.96. Not displaying the output, this is coded as follow, using .P instead of .p to differentiate it from the first GLM:

summary(m.walleye.P<-glmmTMB(count~age,family=poisson,data=walleye))

The "response" residuals correspond to the simplest form, i.e. the difference between the observed and predicted values:

observed<-walleye$count
predicted<-fitted(m.walleye.P)

resid<-observed-predicted

The "response" residuals obained (i.e., resid object) would be identical as those:

residuals_response<-residuals(m.walleye.P,type="response")
resid == residuals_response
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

The function to transform the "response" to "deviance" residuals for family=poisson is based on the "observed" and "predicted" objects. Note that the function needs to be modified according to the considered error structure, the Poisson family representing (I believe) the simplest case when analyzing count data:

Poisson.dev<-function(y, mu)
2*(y*log(ifelse(y == 0, 1, y/mu))-(y-mu))

residuals.deviance<- sqrt(Poisson.dev(observed,predicted))*ifelse(observed>predicted,1,-1)

The Residual deviance is finally obtained as:

sum(residuals.deviance^2)
[1] 111.9586

It's the correct value.

Here's a last check to make sure that we would get the same result for m.walleye.p from the "deviance" residuals (see this link) that can be directly extracted from a glm() object:

sum(residuals(m.walleye.p,type="deviance")^2)
[1] 111.9586

It's the same value as:

m.walleye.p$deviance
[1] 111.9586

And the Null deviance would be:

m.walleye.p$null.deviance
[1] 851.0357

So the deviance explained (D2) would be:

100*(1-m.walleye.p$deviance/m.walleye.p$null.deviance)
[1] 86.84443

Although this value is pretty high, this model is completely inadequate. The count data are over-dispersed and thus, the required equi-dispersion (variance = mean) for a Poisson model is not respected. To be adequate, about 95% of the residuals should be found within the hnp simulated envelope below if the distributional assumptions of the Poisson family were to be respected, which is clearly not the case here, as 80% are found outside:

library(hnp)
hnp(m.walleye.p,resid.type="pearson",how.many.out=TRUE,paint=TRUE)
Poisson model 
Total points: 15 
Points out of envelope: 12 ( 80 %)

enter image description here

So, in the end, contrasting model adequacy (hnp) and model approximated explanatory power (D2) is valuable for more reliable statistical inferences.

For glmmTMB, one need to also run the null counterpart of the model (i.e., ~ 1), as indicated by @Sal Mangiafico and @Ben Bolker, to obtain the Null deviance too and then, calculate D2 from these values.

Maybe at some point the glmmTMB package will allow to obtain these deviance-related metrics, but the solution here is a bit too long to apply and I'm afraid that the required functions for either the Generalized Poisson (family=genpois), mean-parameterized Conway-Maxwell-Poisson (family=compois) or esle such as family=nbinom1 may be too difficult to code under this proposed approach.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.