# Poisson GLMM vs GLM in R (lme4)

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)

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.

• Seems more like a statistics question than a programming question so it should probably be asked at Cross Validated rather than Stack Overflow. Jul 18, 2017 at 5:13
• In taking a look at the data, there doesn't appear to be a nested data structure, unless somehow count is supposed to be nested within race. Typically, there is a more obvious nesting structure. For instance, each row of data could represent homicides in a certain zip code and researchers may want to nest those observations within cities. Not knowing the specific example from the book it is hard for me to see how the nesting is supposed to work in this case. Jul 19, 2017 at 2:50

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

• Thanks for the comment. Maybe something is wrong with the question. I am almost positive that a GLMM must be involved though because this chapter is all about modeling correlated responses.
– Alex
Jul 18, 2017 at 4:14

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.