Binary, Ordinal, Multinomial, and Poisson Approaches
This lab examines a single \(3 \times 3\) contingency table through the lens of four different regression models: binary logistic, proportional odds ordinal logistic, multinomial logistic, and Poisson log-linear. Each model makes different assumptions about the outcome, yet all draw inference from the same underlying table of counts. The goals are to:
We will also consider the role of covariate adjustment (age and CHD) and close with a Poisson log-linear formulation that unifies contingency table analysis with rate modeling and survival analysis.
The data form a \(3 \times 3\) table with a 3-level outcome y (0, 1, 2) and a 3-level predictor x (0, 1, 2).
# Recreate the dataset from the cell counts
d1 <- expand.grid(y = 0:2, x = 0:2, wt = NA)
d1$wt <- c(30, 15, 45, 18, 12, 60, 10, 8, 72)
# Expanded (one-row-per-subject) dataset
d2 <- data.frame(
x = rep(d1$x, d1$wt),
y = rep(d1$y, d1$wt)
)
# Verify the table (reversed so highest levels appear top-left)
table(factor(d2$x, levels = 2:0), factor(d2$y, levels = 2:0))
2 1 0
2 72 8 10
1 60 12 18
0 45 15 30We also create dichotomized and factor versions of the variables that will be needed across models.
d2$y.bin <- (d2$y >= 1) + 0
d2$x.bin <- (d2$x >= 1) + 0
d2$x.factor <- factor(d2$x)
d2$y.factor <- factor(d2$y)A utility function for computing logits by hand:
mylogit <- function(p) log(p / (1 - p))Here we collapse the 3-level outcome to binary (\(Y \ge 1\) vs. \(Y = 0\)) and the predictor to binary (\(X \ge 1\) vs. \(X = 0\)). This yields a \(2 \times 2\) table.
library(rms)Loading required package: Hmisc
Attaching package: 'Hmisc'The following object is masked _by_ '.GlobalEnv':
knitrSetThe following objects are masked from 'package:base':
format.pval, units# Fit via lrm
lrm(y.bin ~ x.bin, data = d2)Logistic Regression Model
lrm(formula = y.bin ~ x.bin, data = d2)
Model Likelihood Discrimination Rank Discrim.
Ratio Test Indexes Indexes
Obs 270 LR chi2 10.77 R2 0.060 C 0.617
0 58 d.f. 1 R2(1,270)0.036 Dxy 0.234
1 212 Pr(> chi2) 0.0010 R2(1,136.6)0.069 gamma 0.462
max |deriv| 1e-05 Brier 0.162 tau-a 0.079
Coef S.E. Wald Z Pr(>|Z|)
Intercept 0.6931 0.2236 3.10 0.0019
x.bin 0.9985 0.3038 3.29 0.0010 # Fit via glm (for comparison)
glm(y.bin ~ x.bin, family = "binomial", data = d2)
Call: glm(formula = y.bin ~ x.bin, family = "binomial", data = d2)
Coefficients:
(Intercept) x.bin
0.6931 0.9985
Degrees of Freedom: 269 Total (i.e. Null); 268 Residual
Null Deviance: 280.9
Residual Deviance: 270.2 AIC: 274.2# Marginal counts for the collapsed 2x2 table
table(factor(d2$y.bin, levels = 1:0), factor(d2$x.bin, levels = 1:0))
1 0
1 152 60
0 28 30# Intercept = logit(Pr(Y>=1 | X=0))
# Among X=0: 60 events out of 90
mylogit(60 / 90)[1] 0.6931472# Slope = log odds ratio
# OR = (a*d) / (b*c) from the 2x2 table
# X=0: Y>=1 has 60, Y=0 has 30
# X=1: Y>=1 has 152, Y=0 has 28 (pooling x=1 and x=2)
log(30 * 152 / (28 * 60))[1] 0.9985288The binary logistic model with a single binary predictor has 2 parameters (intercept + slope) for a \(2 \times 2\) table that has 2 free probabilities (one per row, since probabilities sum to 1 within rows). So it is saturated for the collapsed \(2 \times 2\) table — it reproduces those cell proportions exactly.
But the question of interest is whether the model is saturated for the original \(3 \times 3\) table, and the answer is no. The \(3 \times 3\) table has \(3 \times 2 = 6\) free conditional probabilities (2 per row). The binary logistic model estimates only 2 parameters. By collapsing both the outcome and predictor from 3 levels to 2, we have discarded information about the internal structure of the table. The model cannot recover the individual cell probabilities of the \(3 \times 3\) table — it only describes a coarsened summary. Whether that loss matters depends on whether the collapsed categories are scientifically meaningful.
The multinomial model treats all 3 outcome levels as nominal categories, using \(Y = 0\) as the reference. With x.factor (3 levels), the model has a separate set of log-odds for each non-reference outcome level at each predictor level.
Note: The multinomial model uses a log link on the ratio of category probabilities (i.e., log relative risk or log-odds relative to the reference category), not a logit link. Each equation models \(\log[\Pr(Y = j \mid X) / \Pr(Y = 0 \mid X)]\) as a linear function of \(X\).
\[ \log[\Pr(Y = 1 \mid X=x) / \Pr(Y = 0 \mid X=x)] = \beta_{01} + \beta_{11}I(X=1) + \beta_{21}I(X=2) \] * Isolate \(\beta_{21}\)
\[ \log[\Pr(Y = 1 \mid X=2) / \Pr(Y = 0 \mid X=2)] = \beta_{01} + \beta_{21} \] * \(x=0\)
\[ \log[\Pr(Y = 1 \mid X=0) / \Pr(Y = 0 \mid X=0)] = \beta_{01} \] $$ {21} = - = {01}
$$
library(nnet)
m_multi <- multinom(y ~ x.factor, data = d2)# weights: 12 (6 variable)
initial value 296.625318
iter 10 value 225.904554
final value 225.904535
convergedsummary(m_multi)Call:
multinom(formula = y ~ x.factor, data = d2)
Coefficients:
(Intercept) x.factor1 x.factor2
1 -0.6931419 0.2876678 0.4700011
2 0.4054601 0.7985043 1.5686261
Std. Errors:
(Intercept) x.factor1 x.factor2
1 0.3162270 0.4887619 0.5700878
2 0.2357023 0.3574597 0.4116369
Residual Deviance: 451.8091
AIC: 463.8091 The multinomial intercepts and slopes can be computed directly from the cell counts. For outcome level \(j\) relative to \(Y = 0\):
\[ \beta_{j0} = \log\left(\frac{n_{0j}}{n_{00}}\right), \qquad \beta_{jk} = \log\left(\frac{n_{kj}}{n_{k0}}\right) - \log\left(\frac{n_{0j}}{n_{00}}\right) \]
where \(n_{kj}\) is the count in predictor level \(k\), outcome level \(j\).
# Slopes for Y=2 (relative to Y=0)
# beta_21 (x=1 effect): log(n_12/n_10) - log(n_02/n_00)
log(72/90 / (30/90)) - log(45/90 / (30/90))[1] 0.4700036# Slopes for Y=1 (relative to Y=0)
log(8/90 / (10/90)) - log(15/90 / (30/90))[1] 0.4700036# Equivalently, beta_21 is a log odds ratio from the 2x2 sub-table
# comparing (x=2, y=1) vs (x=0, y=0):
with(d2[d2$x %in% c("0", "2") & d2$y %in% c("0", "1"), ],
table(factor(x, levels = c("2","0")), factor(y, levels = c("1","0"))))
1 0
2 8 10
0 15 30log(8 * 30 / (15 * 10)) # log OR[1] 0.4700036# This would be the easier way to obtain the other betas, e.g. beta_22
with(d2[d2$x %in% c("0", "2") & d2$y %in% c("0", "2"), ],
table(factor(x, levels = c("2","0")), factor(y, levels = c("2","0"))))
2 0
2 72 10
0 45 30log(72 * 30 / (45 * 10)) # log OR[1] 1.568616With 3 predictor levels and 3 outcome levels, the table has \(3 \times 2 = 6\) free probabilities (2 per row, since the 3 outcome probabilities sum to 1 within each row). The multinomial model with x.factor estimates \(2 \times 3 = 6\) parameters (an intercept and two slopes for each of the 2 non-reference outcome levels). The model is saturated — it will always reproduce the observed cell proportions exactly, regardless of the data.
The proportional odds model treats the outcome as ordinal and assumes a common odds ratio across all binary cutpoints of \(Y\).
m_ord <- lrm(y ~ x.factor, data = d2, x = TRUE)
m_ordLogistic Regression Model
lrm(formula = y ~ x.factor, data = d2, x = TRUE)
Model Likelihood Discrimination Rank Discrim.
Ratio Test Indexes Indexes
Obs 270 LR chi2 19.10 R2 0.083 C 0.637
0 58 d.f. 2 R2(2,270)0.061 Dxy 0.274
1 35 Pr(> chi2) <0.0001 R2(2,190.7)0.086 gamma 0.401
2 177 Brier 0.211 tau-a 0.139
max |deriv| 6e-05
Coef S.E. Wald Z Pr(>|Z|)
y>=1 0.6931 0.2108 3.29 0.0010
y>=2 0.0000 0.2049 0.00 1.0000
x.factor=1 0.6931 0.2974 2.33 0.0198
x.factor=2 1.3863 0.3310 4.19 <0.0001 The proportional odds model has two intercepts (one for each cutpoint of \(Y\)) and, with x.factor, two slopes (\(X = 1\) and \(X = 2\) effects). But the key constraint is that each slope is the same across cutpoints.
We can verify the proportional odds assumption by checking whether the log-odds ratios are the same at each cutpoint:
# At cutpoint Y >= 1 vs Y = 0:
mylogit(30/90) - mylogit(10/90)[1] 1.386294# At cutpoint Y >= 2 vs Y <= 1:
mylogit(45/90) - mylogit(18/90)[1] 1.386294# These should be equal if proportional odds holds
# (and in this dataset, they are — by construction)The proportional odds model with x.factor (3 levels) has \(2 + 2 = 4\) parameters (2 intercepts + 2 slopes). But the \(3 \times 3\) table has 6 free probabilities. So this model is not saturated in general — it has 2 residual degrees of freedom.
However, in this particular dataset, proportional odds holds exactly (the log-odds ratios are identical across cutpoints), so the model happens to reproduce the cell probabilities perfectly. This is a property of this dataset, not of the model itself.
# Predicted probabilities from both models — they agree here
newdata <- data.frame(x.factor = levels(d2$x.factor))
cat("Multinomial predicted probabilities:\n")Multinomial predicted probabilities:predict(m_multi, newdata = newdata, "probs") 0 1 2
1 0.3333339 0.16666782 0.4999983
2 0.2000014 0.13333304 0.6666656
3 0.1111106 0.08888874 0.8000006cat("\nOrdinal predicted probabilities:\n")
Ordinal predicted probabilities:predict(m_ord, newdata, type = "fitted") y>=1 y>=2
1 0.6666667 0.5000000
2 0.8000000 0.6666667
3 0.8888889 0.8000000The fact that these predictions agree is because \(2\hat\beta_1 = \hat\beta_2\) happens to hold in the ordinal model. This is a coincidence of this dataset, not a general result.
Should we adjust for age?
Key considerations: Is age a confounder? Is it associated with both \(X\) and \(Y\)? If so, adjustment may be warranted to obtain a less biased estimate of the \(X\)–\(Y\) association. However, if the study is designed so that \(X\) is randomized or balanced on age, adjustment may not be necessary (age could be a precision variable, which in the case of logistic regression tends to drive associations away from the null)
Should we adjust for CHD?
Key considerations: Is CHD a mediator on the causal pathway from \(X\) to \(Y\)? If so, adjusting for it would block part of the causal effect and bias the total effect estimate. Is CHD a collider? Adjusting for a collider can introduce spurious associations. The decision depends on the causal structure. A DAG is helpful here.
So far, we have modeled the conditional distribution of \(Y\) given \(X\) using logistic, ordinal, and multinomial models. Each conditions on the row margins (the number of subjects at each level of \(X\)) and models the probability of falling into each outcome category.
There is an alternative: instead of conditioning on the margins, we can model the cell counts directly using Poisson regression. This is the log-linear model approach to contingency tables.
Why would we do this? Several reasons:
We reshape the \(3 \times 3\) table into a “long” dataset with one row per cell (9 rows), a count variable, and factor variables for \(X\) and \(Y\).
The saturated Poisson log-linear model is:
\[ \log(\mu_{ij}) = \lambda + \alpha_i + \beta_j + \gamma_{ij} \]
where:
The interaction terms \(\gamma_{ij}\) are the parameters of scientific interest. They encode the same information as the log-odds ratios from the multinomial model.
# Cell-level data for the Poisson model
d_cells <- d1
d_cells$x.factor <- factor(d_cells$x)
d_cells$y.factor <- factor(d_cells$y)
d_cells$count <- d_cells$wt
d_cells y x wt x.factor y.factor count
1 0 0 30 0 0 30
2 1 0 15 0 1 15
3 2 0 45 0 2 45
4 0 1 18 1 0 18
5 1 1 12 1 1 12
6 2 1 60 1 2 60
7 0 2 10 2 0 10
8 1 2 8 2 1 8
9 2 2 72 2 2 72# Saturated Poisson log-linear model
m_pois <- glm(count ~ x.factor * y.factor, family = poisson, data = d_cells)
summary(m_pois)
Call:
glm(formula = count ~ x.factor * y.factor, family = poisson,
data = d_cells)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.4012 0.1826 18.629 < 2e-16 ***
x.factor1 -0.5108 0.2981 -1.713 0.086646 .
x.factor2 -1.0986 0.3651 -3.009 0.002624 **
y.factor1 -0.6931 0.3162 -2.192 0.028385 *
y.factor2 0.4055 0.2357 1.720 0.085388 .
x.factor1:y.factor1 0.2877 0.4888 0.589 0.556135
x.factor2:y.factor1 0.4700 0.5701 0.824 0.409689
x.factor1:y.factor2 0.7985 0.3575 2.234 0.025494 *
x.factor2:y.factor2 1.5686 0.4116 3.811 0.000139 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 1.4144e+02 on 8 degrees of freedom
Residual deviance: -2.8866e-15 on 0 degrees of freedom
AIC: 62.671
Number of Fisher Scoring iterations: 3The main effects for x.factor capture how the row totals differ — they absorb the marginal distribution of the predictor. The main effects for y.factor capture how the column totals differ — they absorb the marginal distribution of the outcome. Neither set of main effects tells us about the association between \(X\) and \(Y\).
The interaction terms are where the action is. Each interaction parameter is a log odds ratio comparing a specific (\(X\), \(Y\)) cell to the reference cell, relative to what would be expected under independence.
# Compare interaction parameters to multinomial slopes
cat("Poisson interaction coefficients:\n")Poisson interaction coefficients:coef(m_pois)[grep(":", names(coef(m_pois)))]x.factor1:y.factor1 x.factor2:y.factor1 x.factor1:y.factor2 x.factor2:y.factor2
0.2876821 0.4700036 0.7985077 1.5686159 cat("\nMultinomial coefficients:\n")
Multinomial coefficients:coef(m_multi) (Intercept) x.factor1 x.factor2
1 -0.6931419 0.2876678 0.4700011
2 0.4054601 0.7985043 1.5686261The Poisson interaction terms and the multinomial regression coefficients encode the same log-odds ratios — they are just organized differently. The multinomial model parameterizes them as slopes within each outcome equation; the Poisson model parameterizes them as interaction terms in a single equation.
This equivalence is sometimes called the Birch result or the Poisson trick. The key insight:
If you model cell counts as independent Poisson random variables with a log-linear model, and then condition on the observed row totals (the marginal counts of \(X\)), the conditional distribution of the cell counts within each row is multinomial. Since the row main effects \(\alpha_i\) perfectly absorb the row margins, the Poisson likelihood factors into a part that depends only on the main effects and a part that depends only on the interaction — and the interaction part is exactly the multinomial likelihood.
In practice, this means:
The saturated Poisson log-linear model for a \(3 \times 3\) table has \(3 \times 3 = 9\) parameters (intercept + 2 row effects + 2 column effects + 4 interactions) for 9 cell counts. It is saturated — it reproduces every cell count exactly.
# Fitted values should exactly match observed counts
cbind(observed = d_cells$count, fitted = fitted(m_pois)) observed fitted
1 30 30
2 15 15
3 45 45
4 18 18
5 12 12
6 60 60
7 10 10
8 8 8
9 72 72If we drop the interaction terms (the independence model), we have \(1 + 2 + 2 = 5\) parameters for 9 cells, leaving 4 residual degrees of freedom. This reduced model tests the hypothesis of independence between \(X\) and \(Y\) — it is the log-linear analogue of the chi-squared test.
# Independence model (no interaction)
m_pois_indep <- glm(count ~ x.factor + y.factor, family = poisson, data = d_cells)
summary(m_pois_indep)
Call:
glm(formula = count ~ x.factor + y.factor, family = poisson,
data = d_cells)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.962e+00 1.570e-01 18.865 < 2e-16 ***
x.factor1 -1.581e-12 1.491e-01 0.000 1.0000
x.factor2 4.399e-12 1.491e-01 0.000 1.0000
y.factor1 -5.051e-01 2.140e-01 -2.360 0.0183 *
y.factor2 1.116e+00 1.513e-01 7.374 1.65e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 141.442 on 8 degrees of freedom
Residual deviance: 19.096 on 4 degrees of freedom
AIC: 73.767
Number of Fisher Scoring iterations: 4# Likelihood ratio test for independence
anova(m_pois_indep, m_pois, test = "LRT")Analysis of Deviance Table
Model 1: count ~ x.factor + y.factor
Model 2: count ~ x.factor * y.factor
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 4 19.096
2 0 0.000 4 19.096 0.0007525 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In many applications of Poisson regression, an offset term appears:
\[ \log(\mu_i) = \log(t_i) + \mathbf{x}_i'\boldsymbol{\beta} \]
where \(t_i\) is the exposure (person-time, population size, etc.). The offset standardizes the counts so that we model a rate (\(\mu_i / t_i\)) rather than a raw count. Without the offset, we would confound “more events” with “more opportunity to observe events.”
In our contingency table, no offset is needed. Why?
Every subject contributes exactly one observation — they enter the study, get classified into one cell, and that is it. There is no variation in the opportunity to be counted. Every person had the same “exposure”: one unit of observation, one chance to land in a cell.
The implicit offset is \(\log(1) = 0\) for every cell, and so it vanishes from the model. The row and column main effects handle the marginal totals, which is the contingency table analogue of what an offset does in the rate context.
To make this concrete, consider what happens if we do use an offset equal to the log of the row totals:
# With offset = log(row total), we model conditional probabilities
d_cells$row_total <- ave(d_cells$count, d_cells$x.factor, FUN = sum)
m_pois_offset <- glm(count ~ x.factor * y.factor + offset(log(row_total)),
family = poisson, data = d_cells)
# The interaction terms are unchanged
cat("Without offset (interaction terms):\n")Without offset (interaction terms):coef(m_pois)[grep(":", names(coef(m_pois)))]x.factor1:y.factor1 x.factor2:y.factor1 x.factor1:y.factor2 x.factor2:y.factor2
0.2876821 0.4700036 0.7985077 1.5686159 cat("\nWith offset (interaction terms):\n")
With offset (interaction terms):coef(m_pois_offset)[grep(":", names(coef(m_pois_offset)))]x.factor1:y.factor1 x.factor2:y.factor1 x.factor1:y.factor2 x.factor2:y.factor2
0.2876821 0.4700036 0.7985077 1.5686159 The interaction terms — the association parameters — are identical with or without the offset. The offset only changes the main effect estimates, which are nuisance parameters for the association question. When the offset equals the log row total, the model explicitly conditions on the row margins and the main effects now represent the conditional log-probabilities rather than the marginal log-counts. But the association is the same either way.
This is the Birch result in action: inference about association is invariant to whether you treat the margins as fixed (offset / multinomial) or random (Poisson with main effects).
The Poisson distribution assumes that the variance equals the mean: \(\text{Var}(Y_i) = \mu_i\). In many applications of Poisson regression — particularly with individual-level count data — this assumption is suspect. Counts often exhibit overdispersion (variance exceeding the mean) due to unobserved heterogeneity, clustering, or other sources of extra-Poisson variation. In those settings, the model-based standard errors from glm() are too small, and it is common practice to use robust (sandwich) standard errors to obtain valid inference without relying on the variance assumption.
In the contingency table setting, this concern does not arise. Here is why: the saturated Poisson log-linear model has one parameter per cell count — 9 parameters for 9 observations. The fitted values reproduce the observed counts exactly. There are zero residual degrees of freedom and no residual deviance. The mean–variance relationship is not testable because there is no residual variation left to assess. The model fits the data perfectly by construction.
# Residual deviance and df for the saturated model
cat("Residual deviance:", deviance(m_pois), "\n")Residual deviance: -2.88658e-15 cat("Residual df:", df.residual(m_pois), "\n")Residual df: 0 # For comparison, the independence model does have residual df
# and we can check for overdispersion there
cat("\nIndependence model:\n")
Independence model:cat("Residual deviance:", deviance(m_pois_indep), "\n")Residual deviance: 19.09603 cat("Residual df:", df.residual(m_pois_indep), "\n")Residual df: 4 cat("Dispersion estimate:", deviance(m_pois_indep) / df.residual(m_pois_indep), "\n")Dispersion estimate: 4.774007 For the independence model (which is not saturated), overdispersion could in principle be relevant — if the residual deviance substantially exceeds the residual degrees of freedom, the Poisson variance assumption may be violated. In that case, one could either use family = quasipoisson (which estimates a dispersion parameter) or compute robust standard errors using the sandwich estimator.
# Robust SEs via sandwich estimator (for the independence model)
library(sandwich)
library(lmtest)Loading required package: zoo
Attaching package: 'zoo'The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Attaching package: 'lmtest'The following object is masked from 'package:rms':
lrtestcat("Model-based SEs (independence model):\n")Model-based SEs (independence model):coeftest(m_pois_indep)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.9618e+00 1.5700e-01 18.8653 < 2.2e-16 ***
x.factor1 -1.5814e-12 1.4907e-01 0.0000 1.00000
x.factor2 4.3992e-12 1.4907e-01 0.0000 1.00000
y.factor1 -5.0509e-01 2.1404e-01 -2.3598 0.01828 *
y.factor2 1.1157e+00 1.5130e-01 7.3743 1.653e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\nRobust (sandwich) SEs:\n")
Robust (sandwich) SEs:coeftest(m_pois_indep, vcov = vcovHC(m_pois_indep, type = "HC0"))
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.9618e+00 3.1512e-01 9.3991 < 2.2e-16 ***
x.factor1 -1.5814e-12 1.9993e-01 0.0000 1.00000
x.factor2 4.3992e-12 2.6999e-01 0.0000 1.00000
y.factor1 -5.0509e-01 2.8352e-01 -1.7815 0.07483 .
y.factor2 1.1157e+00 2.6819e-01 4.1601 3.182e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The robust standard errors protect against misspecification of the variance function. When the mean model is correct but the variance assumption is wrong, the point estimates from Poisson regression remain consistent (they solve the correct estimating equation), but the standard errors need correction. This is the same logic as using robust SEs in ordinary least squares when heteroscedasticity is present.
For the saturated model, robust standard errors are not needed — the model is exact. But the concept becomes important when you move to non-saturated Poisson models (e.g., with continuous predictors, or the independence model) or to the rate and survival settings discussed above, where the Poisson variance assumption is a real constraint. The sandwich estimator is the standard tool for relaxing that constraint while keeping the convenient Poisson regression framework.
| Model | Outcome Treatment | Predictor Treatment | Parameters | Free Cells | Saturated |
|---|---|---|---|---|---|
| Binary Logistic | Dichotomized (Y >= 1 vs Y = 0) | Binary (X >= 1 vs X = 0) | 2 | 2 | Yes (for 2x2 only; No for 3x3) |
| Proportional Odds | Ordinal (cumulative logits) | Factor (3 levels) | 4 | 6 | No (in general) |
| Multinomial | Nominal (3 categories, Y=0 ref) | Factor (3 levels) | 6 | 6 | Yes |
| Poisson Log-Linear | Cell counts (9 cells) | Factor (3 levels) | 9 | 9 | Yes |
Binary Logistic: Simple, familiar, and well-understood. But collapsing the outcome loses information — we can no longer distinguish between outcome levels 1 and 2. The choice of where to dichotomize is often arbitrary and can affect conclusions.
Proportional Odds: Preserves the ordinal structure and is parsimonious (fewer parameters). The proportional odds assumption constrains the model, which is a strength when the assumption holds (more power, simpler interpretation) and a weakness when it does not (biased estimates, poor fit). Not saturated for \(> 2\) predictor levels, so model fit should be checked.
Multinomial: Makes no assumptions about the outcome structure — fully flexible. Always saturated when the predictor is categorical. But this flexibility comes at a cost: more parameters to estimate, and no borrowing of strength across outcome levels. With continuous predictors or sparse data, the multinomial model can be unstable.
Poisson Log-Linear: Provides a unified framework that connects to rate modeling and survival analysis. The saturated model is equivalent to the multinomial. The independence model provides a global test of association. The offset mechanism generalizes naturally to settings with unequal exposure. The main disadvantage for contingency table analysis per se is that the parameterization (main effects + interactions) is less directly interpretable than the conditional probabilities of the multinomial model — though the interaction terms are the same log-odds ratios.