I am trying to understand how to interpret log-linear models for contingency tables, fitted by way of Poisson GLMs.

Consider this example from CAR (Fox and Weisberg, 2011, p. 252).

(tab.sex.citizen <- xtabs(count ~ sex + citizen, data=AMSsurvey))


sex      Non-US  US
  Female    260 202
  Male      501 467

Then we fit the model of (mutual) independence:

AMS2 <- as.data.frame(tab.sex.citizen)
(phd.mod.indep <- glm(Freq ~ sex + citizen, family=poisson, data=AMS2))
pchisq(2.57, df=1, lower.tail=FALSE)


> (phd.mod.indep <- glm(Freq ~ sex + citizen, family=poisson, data=AMS2))

Call:  glm(formula = Freq ~ sex + citizen, family = poisson, data = AMS2)

(Intercept)      sexMale    citizenUS  
     5.5048       0.7397      -0.1288  

Degrees of Freedom: 3 Total (i.e. Null);  1 Residual
Null Deviance:      191.5 
Residual Deviance: 2.572    AIC: 39.16
> pchisq(2.57, df=1, lower.tail=FALSE)
[1] 0.1089077

The p value is close to 0.1 indicating weak evidence to reject independence. However, let us assume that we have sufficient evidence to reject the NULL (i.e. for our purposes, the 0.10 p value is indicative of an association between the two variables).

Question: How, then, do we interpret this loglinear model?

(Do we fit the saturated model (i.e. update(phd.mod.indep, . ~ . + sex:citizen))? Do we interpret the estimated regression coefficients? In CAR they stop at this point, because of weak evidence for rejecting the NULL, but I'm interested in understanding the mechanics of the interpretation of this simple log-linear model as if the "interaction" were significant...)


In general, there isn't much to a 2-way contingency table, but you are trying to unpack this at such a level of detail that some confusions are arising. Typically, with a simple contingency table like this, people just want to know if the variables (sex and citizen) are independent. To assess that, you can run a chi-squared test:

#  Pearson's Chi-squared test with Yates' continuity correction
# data:  tab.sex.citizen
# X-squared = 2.389, df = 1, p-value = 0.1222

You can perform a likelihood ratio version of this test (the chi-squared test above is a score test), by performing a nested model test of the Poisson GLM you ran against the full (saturated) model:

phd.mod.indep <- glm(Freq ~ sex + citizen, family=poisson, data=AMS2)
phd.mod.sat   <- glm(Freq ~ sex * citizen, family=poisson, data=AMS2)
anova(phd.mod.indep, phd.mod.sat, test="LRT")
# Analysis of Deviance Table
# Model 1: Freq ~ sex + citizen
# Model 2: Freq ~ sex * citizen
#   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
# 1         1     2.5721                     
# 2         0     0.0000  1   2.5721   0.1088
deviance(phd.mod.indep)                                        # [1] 2.572123
deviance(phd.mod.sat)                                          # [1] 3.308465e-14
1-pchisq(deviance(phd.mod.indep)-deviance(phd.mod.sat), df=1)  # [1] 0.1087617

You can also get the Wald test of the interaction term (which is the test of independence):

# ...
# Coefficients:
#                   Estimate Std. Error z value Pr(>|z|)    
# (Intercept)        5.56068    0.06202  89.663  < 2e-16 ***
# sexMale            0.65592    0.07643   8.582  < 2e-16 ***
# citizenUS         -0.25241    0.09379  -2.691  0.00712 ** 
# sexMale:citizenUS  0.18214    0.11373   1.602  0.10926    
# ...
#     Null deviance: 1.9148e+02  on 3  degrees of freedom
# Residual deviance: 3.3085e-14  on 0  degrees of freedom
# AIC: 38.586
# ...

(To read about score vs Wald vs likelihood ratio tests, see my answer here: Why do my p-values differ between logistic regression output, chi-squared test, and the confidence interval for the OR?)

Note however, that the way you conducted your test of phd.mod.indep is incorrect (see here: Test GLM model using null and model deviances). The test of that model against the null model is the test of whether all cells have the same probability in the population. It would be implemented as follows:

1-pchisq(191.5-2.57, df=3-1)  # [1] 0

Setting aside the test of whether all cell probabilities are equal (which I doubt is of much interest to you), if you believed that there was an association between sex and citizen, you would not interpret the model phd.mod.indep. That would be a misspecified model.

  • $\begingroup$ So I'm looking too much into it, and it's mostly about detecting association... Thank you for the thorough answer. I'll do some reading and some checks, and later get back on the subject of the incorrect test... $\endgroup$ – landroni Jun 19 '15 at 8:20
  • 1
    $\begingroup$ "Note however, that the way you conducted your test of phd.mod.indep is incorrect" I suspect that the difference comes from testing different hypotheses. The 1 - pchisq(191.5-2.57, df=3-1) test will look at whether "the [fitted] model as a whole is better than the null model". The 1 - pchisq(2.57, df=1) test is intended as a lack-of-fit test (see CAR p. 252); in my understanding, it compares the fitted model to the saturated model (i.e. update(phd.mod.indep, . ~ . + sex:citizen), which has 0 residual deviance and 0 df). $\endgroup$ – landroni Jun 19 '15 at 13:59
  • $\begingroup$ @landroni, fair enough. I suppose you can do it that way. $\endgroup$ – gung Jun 19 '15 at 14:02
  • $\begingroup$ I added an answer covering the use of mosaic plots to aid interpretation of contingency tables. $\endgroup$ – landroni Jun 21 '15 at 8:20

In addition to @gung's answer, I discovered that mosaic plots are very useful in assessing the relative relationships with categorical data. Once we establish (actually, we assume here) departure from independence, one can use package vcd to graphically display such contingency tables.

A mosaic plot:

mosaic(tab.sex.citizen, shade=TRUE, legend=TRUE)

enter image description here

An association plot:

assoc(tab.sex.citizen, shade=TRUE, legend=TRUE)

enter image description here

From the second plot, it's clear that the proportion of women is bigger for non-citizens than for US citizens. Equivalently, the proportion of non-citizens is bigger for females than for males.

The graphs are clearer when the departure from independence is clear-cut, so here's an example on the Titanic data (see this tutorial).

sex.survived <- margin.table(Titanic, c(2,4))


> ftable(sex.survived)
       Survived   No  Yes
Male            1364  367
Female           126  344

Then the plots:

mosaic(sex.survived, shade=TRUE, legend=TRUE)
assoc(sex.survived, shade=TRUE, legend=TRUE)

enter image description here

enter image description here

From the Pearson residuals (a measure of the departure of the Observed Frequencies from the Expected Frequencies, or the bits in the data not explained by the loglinear model), it's clear that for males there is a disproportionately high number that haven't survived and a disproportionately low number that have survived. Equivalently, out of those that survived the number of females is unexpectedly high.

  • 2
    $\begingroup$ +1, looking at your data is always important. Mosaic plots & association plots are a nice basic option; there are others as well, see here. Friendly's Visualizing Categorical Data (the book) has lots of nice ideas, this pdf is a nice short introduction. If your contingency tables are larger than 2x2, correspondence analysis can be helpful (see here). $\endgroup$ – gung Jun 21 '15 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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