Poisson GLMM vs GLM in R (lme4)

Question

I have been trying to sharpen my GLMM knowledge by working through some problems in Foundations of Linear and Generalized Linear Models. I am stuck on problem 9.36 which gives some homicide data then states "fit a Poisson GLMM. Interpret estimates. Show that the deviance decreases by 116.6 compared with the Poisson GLM, and intercept"

This is what I did

library(lme4)
homi <- read.table("http://www.stat.ufl.edu/~aa/glm/data/Homicides.dat",
                   header = TRUE)

fit1 <- glm(count~race, family=poisson(link = log),data=homi)

fit2 <- glmer(count~1+(1|race), family=poisson(link = log),data=homi)
summary(fit1)
summary(fit2)

For fit1

Call:
glm(formula = count ~ race, family = poisson(link = log), data = homi)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0218  -0.4295  -0.4295  -0.4295   6.1874  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.38321    0.09713  -24.54   <2e-16 ***
race         1.73314    0.14657   11.82   <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: 962.80  on 1307  degrees of freedom
Residual deviance: 844.71  on 1306  degrees of freedom
AIC: 1122

Number of Fisher Scoring iterations: 6

and for fit2

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: count ~ 1 + (1 | race)
   Data: homi

     AIC      BIC   logLik deviance df.resid 
  1132.5   1142.8   -564.2   1128.5     1306 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.7174 -0.3054 -0.3054 -0.3054 19.3426 

Random effects:
 Groups Name        Variance Std.Dev.
 race   (Intercept) 0.739    0.8596  
Number of obs: 1308, groups:  race, 2

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -1.5235     0.6122  -2.489   0.0128 *

It seems like the deviance for the GLMM is 1128.5 while the deviance for the GLM is 844.71. This is shows the GLM is fitting better than the GLMM which I think is the exact opposite solution from what question implied. I am not sure if I am looking at the correct output or if I setup the problem wrong.

Answer 1

There is probably something wrong with the question. The deviance for fit1 can be computed with

deviance(fit1)

# same as 
sum(resid(fit1)^2)

But for the GLMM, the lme4 packages uses another method documented in ?llikAIC, which gives 1128 - higher than for the GLM.

Perhaps the question actually wants to test for the decrease in deviance when adding the single fixed effect, because it comes out as 118, close to your 116.6:

fit1 <- glm(count~race, family=poisson(link = log),data=homi)
fit0 <- glm(count~1, family=poisson(link = log),data=homi)

> anova(fit0, fit1)
Analysis of Deviance Table

Model 1: count ~ 1
Model 2: count ~ race
Resid. Df Resid. Dev Df Deviance
1      1307     962.80            
2      1306     844.71  1   118.09

Answer 2

This is a tricky exercise. Let's break it down in two parts.

1. Specify a Poisson GLMM for the 1990 General Society Survey data. 1308 subjects responded to the question: Within the past 12 months, how many people have you known personally that were victims of homicide? The responses are broken down by race ("white" and "black").
2. Fit the Poisson GLMM.

This is survey data; it seems natural to let the effects vary across participants. So let's make the GLMM a random-intercepts model.

In mathematical notation we will compare:

• A Poisson GLM with fixed race effects. The intercept for "white" subjects is $\beta_0$ and the intercept for "black" subjects is $\beta_0 + \beta_1$. $$ \begin{aligned} \operatorname{count}_{i} &\sim \operatorname{Poisson}(\lambda_i) \\ \log(\lambda_i) &= \beta_{0} + \beta_{1}(\operatorname{black}) \end{aligned} $$

• A Poisson GLMM with random race effects. The intercepts for "white" subjects are $\operatorname{Normal}(\gamma_0, \sigma^2)$ and the intercepts for "black" subjects are $\operatorname{Normal}(\gamma_0 + \gamma_1, \sigma^2)$. $$ \begin{aligned} \operatorname{count}_{i} &\sim \operatorname{Poisson}(\lambda_i) \\ \log(\lambda_i) &=\alpha_{j[i]} \\ \alpha_{j} &\sim N \left(\gamma_{0}+ \gamma_{1}(\operatorname{black}), \sigma^2\right) \end{aligned} $$ We have only observation from each participant, so we don't expect to get reliable estimates of the random intercepts. This is not a very sophisticated mixed-effects model. In fact, a Negative Binomial GLM fits the data better .

Let's fit the two models: the GLM with stats::glm and the GLMM with lme4::glmer.

fit.glm <- glm(
  count ~ race,
  family = poisson,
  data = homicides
)
fit.glmer <- glmer(
  count ~ race + (1 | Obs),
  family = poisson,
  data = homicides,
  nAGQ = 20
)

And now the tricky part. Let's calculate the deviance two ways: a) two times the negative log likelihood and b) the sum of the deviance residuals squared.

# Compute deviance as -2 times the log likelihood.
-2 * logLik(fit.glm)
#> 1117.99 (df=2)
-2 * logLik(fit.glmer)
#> 728.0926 (df=3)

# Compute deviance as the sum of the deviance residuals squared.
sum(resid(fit.glm)^2)
#> 844.7073
sum(resid(fit.glmer)^2)
#> 214.0758

Notice that 844.7073 - 728.0926 = 116.6147. This gives the "right answer" though the computation is not meaningful. The answer @RemkoDuursma points it out as well.

To choose between the models we can use the deviance = -2 × log-likelihood or the AIC = -2 × log-likelihood + 2 × #parameters. See also Residual deviance, residuals, and log-likelihood in [weighted] logistic regression

PS. I came across this error when I used a different R package to fit the Poisson GLMM, glmmML.

fit.glmmML <- glmmML(
  count ~ race,
  cluster = Obs,
  family = poisson,
  method = "ghq",
  data = homicides
)

deviance(fit.glm)
#> 844.7073
deviance(fit.glmer)
#> 214.0758
deviance(fit.glmmML)
#> 728.107

lme4::lmer and glmmML::glmmML report very different "deviance" for the same Poisson GLMM even though both use maximum likelihood and their parameter estimates are almost the same. It took me a while to realize they have a different definition of deviance.