How to interpret a two-dimensional contingency table?

Question

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).

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

Yielding:

        citizen
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)

Outputting:

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

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

Coefficients:
(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...)

Answer 1

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:

chisq.test(tab.sex.citizen)
#  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):

summary(phd.mod.sat)
# ...
# 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.

Answer 2

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)

An association plot:

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

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).

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

Yielding:

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

Then the plots:

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

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.