Solved – How to interpret a two-dimensional contingency table

categorical datacontingency tablesinterpretationlog-linearr

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

Best Answer

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.

Related Question