2.4 K-Way Tables

These notes rely on PSU STATS 504 course notes.

A k-way table has k independent variables. Each cell in the k-way table has joint probability \(\pi_{ij..k}\). You can collapse a k-way table along a dimension to create a marginal table and the resulting joint probabilities in the marginal table are referred to as marginal distributions. Alternatively, you can consider a single level of a dimension, creating a conditional distribution.

Here are two case studies to illustrate the concepts.

Study 1: Death Penalty. Data is collected for n = 326 murder cases along three dimensions: defendant’s color, victim’s color, and whether or not the defendant received the death penalty.

deathp <- array(c(19, 132, 11, 52, 0, 9, 6, 97), dim=c(2,2,2))
dimnames(deathp) <- list(
  DeathPen = c("yes","no"),
  Defendant=c("white","black"),
  Victim=c("white","black")
)
ftable(deathp, row.vars = c("Defendant", "Victim"), col.vars = "DeathPen")
##                  DeathPen yes  no
## Defendant Victim                 
## white     white            19 132
##           black             0   9
## black     white            11  52
##           black             6  97

Study 2: Boy Scouts. Data is collected for n = 800 Boy Scouts along three dimensions: socioeconomic status, scout status, and delinquency status.

scout <- expand.grid(
  delinquent = c("no","yes"),
  scout = c("no", "yes"),
  SES = c("low", "med","high")
)
scout <- cbind(scout, count = c(169,42,43,11,132,20,104,14,59,2,196,8))

scout_ct <- xtabs(count ~ ., scout)

ftable(scout_ct, row.vars = c("SES", "scout"), col.vars = "delinquent")
##            delinquent  no yes
## SES  scout                   
## low  no               169  42
##      yes               43  11
## med  no               132  20
##      yes              104  14
## high no                59   2
##      yes              196   8

2.4.1 Odds Ratio

In the DeathP study, the marginal odds ratio is 1.18, meaning the odds of death penalty for a white defendant are 1.18 times as high as they are for a black defendant.

deathp_m <- margin.table(deathp, margin = c(1, 2))
deathp_p <- prop.table(deathp_m, margin = 2)
deathp_odds <- deathp_p[1, ] / deathp_p[2, ]
deathp_or <- deathp_odds[1] / deathp_odds[2]
print(deathp_or)
## white 
##   1.2

The conditional odds ratio given the victim is white is 0.68, and 0.79 given the victim is black, meaning the odds of death penalty for a white defendant are 0.68 times as high as they are for a black defendent if the victim is white, and 0.79 times as high as they are for a black defendent if the victim is black.

deathp_m <- margin.table(deathp[, , 1], margin = c(2, 1))
deathp_p <- prop.table(deathp_m, margin = 2)
deathp_odds <- deathp_p[1, ] / deathp_p[2, ]
deathp_or <- deathp_odds[1] / deathp_odds[2]
print(deathp_or)
##  yes 
## 0.68
deathp_m <- margin.table(deathp[, , 2], margin = c(2, 1)) + 0.5
deathp_p <- prop.table(deathp_m, margin = 2)
deathp_odds <- deathp_p[1, ] / deathp_p[2, ]
deathp_or <- deathp_odds[1] / deathp_odds[2]
print(deathp_or)
##  yes 
## 0.79

Notice above that the second margin table had a zero in one cell. To get an odds ratio in this case, the convention is to add 0.5 to all cells.

Interesting that the marginal odds of a white defendent receiving the death penalty are >1, but the conditional odds of a white defendent receiving the death penalty given the race of the victim are <1. Simpson’s paradox is the phenomenon that a pair of variables can have marginal association and partial (conditional) associations in opposite direction. Another way to think about this is that the nature and direction of association changes due to presence or absence of a third (possibly confounding) variable.

2.4.2 Chi-Square Independence Test

There are several ways to think about “independence” when dealing with a k-way table.

  • Mutual independence. All variables are independent from each other, (A, B, C).
  • Joint independence. Two variables are jointly independent of the third, (AB, C).
  • Marginal independence. Two variables are independent if you ignore the third, (A, B).
  • Conditional independence Two variables are independent given the third, (AC, BC).
  • Homogeneous associations Conditional (partial) odds-ratios are not related on the value of the third, (AB, AC, BC).

Under the assumption that the model of independence is true, once you know the marginal probability values, we have enough information to estimate all unknown cell probabilities. Because each of the marginal probability vectors must add up to one, the number of free parameters in the model is (I − 1) + (J − 1) + (K −I ). This is exactly like the two-way table

2.4.2.1 Mutual Independence

The simplest model is that ALL variables are independent of one another (A, B, C). The joint probabilities are equal to the product of the marginal probabilities, \(\pi_{ijk} = \pi_i \pi_j \pi_k\). In terms of odds ratios, the model (A, B, C) implies that the odds ratios in the marginal tables A × B, B × C, and A × C are equal to 1.

The chi-square independence test tests whether observed joint frequency counts \(O_{ijk}\) differ from expected frequency counts \(E_{ijk}\) under the independence model \(\pi_{ijk} = \pi_{i+} \pi_{+j} \pi_{k+}\). \(H_0\) is \(O_{ijk} = E_{ijk}\). This is essentially a one-way table chi-squre goodness-of-fit test.

For the Boy Scouts study, \(X^2\) and its associated p-value is

scout_e <- array(NA, dim(scout_ct))
for (i in 1:dim(scout_ct)[1]) {
  for (j in 1:dim(scout_ct)[2]) {
    for (k in 1:dim(scout_ct)[3]) {
      scout_e[i,j,k] <- (
        margin.table(scout_ct, 3)[k] *
        margin.table(scout_ct, 2)[j] * 
        margin.table(scout_ct, 1)[i]) /
        (sum(scout_ct))^2
    }
  }
}
scout_df <- (prod(dim(scout_ct)) - 1) - sum(dim(scout_ct) - 1)
scout_x2 <- sum((scout_ct - scout_e)^2 / scout_e)
print(scout_x2)
## [1] 215
pchisq(scout_x2, df = scout_df, lower.tail = FALSE)
## [1] 0.00000000000000000000000000000000000000000079

You can also get X^2 this way, although longer code.

x <- scout %>% group_by(delinquent) %>% 
  summarize(n = sum(count)) %>% ungroup() %>% mutate(p = n / sum(n))
## `summarise()` ungrouping output (override with `.groups` argument)
p_d <- x$p
names(p_d) <- x$delinquent

x <- scout %>% group_by(scout) %>% 
  summarize(n = sum(count)) %>% ungroup() %>% mutate(p = n / sum(n))
## `summarise()` ungrouping output (override with `.groups` argument)
p_b <- x$p
names(p_b) <- x$scout

x <- scout %>% group_by(SES) %>% 
  summarize(n = sum(count)) %>% ungroup() %>% mutate(p = n / sum(n))
## `summarise()` ungrouping output (override with `.groups` argument)
p_s <- x$p
names(p_s) <- x$SES

scout_e <- scout %>% 
  mutate(e = as.numeric(p_d[delinquent] * p_b[scout] * p_s[SES] * sum(count)))

x2 <- sum((scout_e$count - scout_e$e)^2 / scout_e$e)
print(x2)
## [1] 215

The deviance statistic is

scout_g2 <- 2 * sum(scout_ct * log(scout_ct / scout_e$e))
print(scout_g2)
## [1] 219
pchisq(scout_g2, df = scout_df, lower.tail = FALSE)
## [1] 0.00000000000000000000000000000000000000000013

Safe to say the mutual independence model does not fit.

chisq.test() does not test the complete independence model for k-way tables. Instead, conduct tests on all \({3 \choose 2} = 3\) embedded 2-way tables. For the Boy Scouts study, you can conduct a series of 2-way chi-sq tests.

SES*scout:

(scout_ses_scout <- margin.table(scout_ct, c(3, 2)))
##       scout
## SES     no yes
##   low  211  54
##   med  152 118
##   high  61 204
chisq.test(scout_ses_scout)
## 
## 	Pearson's Chi-squared test
## 
## data:  scout_ses_scout
## X-squared = 172, df = 2, p-value <0.0000000000000002

SES*delinquent

(scout_ses_delinquent <- margin.table(scout_ct, c(3, 1)))
##       delinquent
## SES     no yes
##   low  212  53
##   med  236  34
##   high 255  10
chisq.test(scout_ses_delinquent)
## 
## 	Pearson's Chi-squared test
## 
## data:  scout_ses_delinquent
## X-squared = 33, df = 2, p-value = 0.00000007

delinquent*scout

(scout_delinquent_scout <- margin.table(scout_ct, c(1, 2)))
##           scout
## delinquent  no yes
##        no  360 343
##        yes  64  33
chisq.test(scout_delinquent_scout)
## 
## 	Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  scout_delinquent_scout
## X-squared = 7, df = 1, p-value = 0.009

The odds ratio for scout delinquency vs non-scout delinquency is 0.54 with 95% CI (0.347, 0.845), so reject the null hypothesis that boy scout status and delinquent status are independent of one another (OR = 1), and thus that boy scout status, delinquent status, and socioeconomic status are not mutually independent.

(scout_or <- vcd::oddsratio(scout_delinquent_scout, log = FALSE))
##  odds ratios for delinquent and scout 
## 
## [1] 0.54
confint(scout_or)
##               2.5 % 97.5 %
## no:yes/no:yes  0.35   0.84

2.4.2.2 Joint Independence

The joint independence model is that two variables are jointly independent of a third (AB, C). The joint probabilities are equal to the product of the AB marginal probability and the C marginal probability, \(\pi_{ijk} = \pi_{ij} \pi_k\).

The degrees of freedom in the associated chi-sq test are the number of free parameters in the saturated model minus the number of free parameters in the joint independence model. \(df = (IJK-1) - ((IJ-1) + (K-1)) = (IJ-1)(K-1)\).

From the Boy Scouts study, suppose juvenile delinquency is the response variable. Test the null hypothesis that juvenile delinquency is independent of boy scout status and socioeconomic status.

scout_d_bs <- ftable(
  scout_ct, 
  row.vars = c("SES", "scout"), 
  col.vars = "delinquent"
)
print(scout_d_bs)
##            delinquent  no yes
## SES  scout                   
## low  no               169  42
##      yes               43  11
## med  no               132  20
##      yes              104  14
## high no                59   2
##      yes              196   8
chisq.test(scout_d_bs)
## 
## 	Pearson's Chi-squared test
## 
## data:  scout_d_bs
## X-squared = 33, df = 5, p-value = 0.000004

Notice how the contingency table is constructed to not indicate whether SES and scout are related.

2.4.2.3 Marginal Independence

The marginal independence model is that two variables are independent while ignoring the third (A, B). The joint probabilities are equal to the product of the two marginal probabilities, \(\pi_{ij} = \pi_{i} \pi_j\), just like a 2-way table.

The degrees of freedom in the associated chi-sq test are the number of free parameters in the saturated model minus the number of free parameters in the marginal independence model. \(df = (I - 1) + (J - 1) = (I-1)(J-1)\).

From the Boy Scouts study, suppose juvenile delinquency is the response variable. Test the null hypothesis that juvenile delinquency is independent of boy scout status and socioeconomic status.

ct <- xtabs(count ~ scout + SES, scout)
print(ct)
##      SES
## scout low med high
##   no  211 152   61
##   yes  54 118  204
chisq.test(ct)
## 
## 	Pearson's Chi-squared test
## 
## data:  ct
## X-squared = 172, df = 2, p-value <0.0000000000000002

2.4.2.4 Conditional Independence

The conditional independence model is that two variables are independent given a third (AB, AC). The joint probabilities are equal to the product of the conditional probabilities (B|A and C|A) and the marginal probability of A, \(\pi_{ijk} = \pi_{j|i}\pi_{k|i}\pi_{i++}\).

You can test for conditional independence (AB, AC) by individually testing the independence (B, C) for each level of A. \(X^2\) and \(G^2\) are the sum of the individual values and the degrees of freedom are \(I(J-1)(K-1)\). A better way to do it is with the Cochran-Mantel-Haenszel Test.

The degrees of freedom in the associated chi-sq test are the number of free parameters in the saturated model minus the number of free parameters in the joint independence model. \(df = (IJK-1) - ((IJ-1) + (K-1)) = (IJ-1)(K-1)\).

From the Boy Scouts study, the section on joint independence found that scouting and SES were jointly independent of delinquency, so it must be that either scouting is independent of delinquency, or SES is independent of delinquency, or both are independent of delinquency. The marginal test below establishes the independance of (scout, delinquent).

ct <- xtabs(count ~ scout + delinquent, scout)
print(ct)
##      delinquent
## scout  no yes
##   no  360  64
##   yes 343  33
chisq.test(ct, correct = FALSE)
## 
## 	Pearson's Chi-squared test
## 
## data:  ct
## X-squared = 7, df = 1, p-value = 0.006

The odds ratio 0.54 is a strong negative relationship between boy scout status and delinquency - boy scouts are 46% less likely (on the odds scale) to be delinquent than non-boy scouts.

ct <- xtabs(count ~ scout + delinquent, scout)
p <- prop.table(ct, margin = 1)
(or <- (p[2, 2] / p[2, 1]) / (p[1, 2] / p[1, 1]))
## [1] 0.54

However, one may ask is scouting and delinquency conditionally independent given SES? You can answer this question with three chi-square tests, one for each level of SES.

ct <- xtabs(count ~ scout + delinquent, scout[scout$SES == "low", ])
print(ct)
##      delinquent
## scout  no yes
##   no  169  42
##   yes  43  11
chisq.test(ct, correct = FALSE)
## 
## 	Pearson's Chi-squared test
## 
## data:  ct
## X-squared = 0.006, df = 1, p-value = 0.9
ct <- xtabs(count ~ scout + delinquent, scout[scout$SES == "med", ])
print(ct)
##      delinquent
## scout  no yes
##   no  132  20
##   yes 104  14
chisq.test(ct, correct = FALSE)
## 
## 	Pearson's Chi-squared test
## 
## data:  ct
## X-squared = 0.1, df = 1, p-value = 0.8
ct <- xtabs(count ~ scout + delinquent, scout[scout$SES == "high", ])
print(ct)
##      delinquent
## scout  no yes
##   no   59   2
##   yes 196   8
chisq.test(ct, correct = FALSE)
## Warning in chisq.test(ct, correct = FALSE): Chi-squared approximation may be
## incorrect
## 
## 	Pearson's Chi-squared test
## 
## data:  ct
## X-squared = 0.05, df = 1, p-value = 0.8

You can add up the three \(X^2\) values and run a chi-sq test with \(df = 3\). p = 0.984 indicating that the conditional independence model fits extremely well.

pchisq(0.0058145 + 0.10098 + 0.053447, df = 3, lower.tail = FALSE)
## [1] 0.98

The other way of testing for independence is the Cochran-Mantel-Haenszel test.

The Cochran-Mantel-Haenszel (CMH) test statistic

\[M^2 = \frac{\left[ \sum_k{(n_{11k} - \mu_{11k})}\right]^2}{\sum_k{Var(n_{11k})}}\]

For the Boy Scouts study, the test is

mantelhaen.test(scout_ct)
## 
## 	Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  scout_ct
## Mantel-Haenszel X-squared = 0.008, df = 1, p-value = 0.9
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.6 1.6
## sample estimates:
## common odds ratio 
##              0.98

X2= 0.008 (df = 1, p-value = 0.9287) indicating the conditional independence model is a good fit for this data, so do not reject the null hypothesis.

2.4.2.5 Homogenous Associations

Whereas the conditional independence model, (AB, AC) requires the BC odds ratios at each level of A to equal 1, the homogeneous association model, (AB, AC, BC), requires the BC odds ratios at each level of A to be identical, but not necessarily equal to 1.

The Breslow-Day statistic is

\[X^2 = \sum_{ijk}{\frac{\left(O_{ijk} - E_{ijk}\right)^2}{E_{ijk}}}\]

DescTools::BreslowDayTest(scout_ct)
## 
## 	Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  scout_ct
## X-squared = 0.2, df = 2, p-value = 0.9