vote_data = matrix(c(175, 16, 54, 188), 2, 2, byrow=TRUE) # Table 11.1
colnames(vote_data) = c('D_2008', 'R_2008')
rownames(vote_data) = c('D_2004', 'R_2004')
vote_data## D_2008 R_2008
## D_2004 175 16
## R_2004 54 188Null is marginal homogeneity \(H_0: \pi_{1+} = \pi_{+1}\)
Difference is \(\delta=\pi_{1+} - \pi_{+1}\)
Estimate
\[ \widehat{\delta} = \widehat{\pi}_{1+} - \widehat{\pi}_{+1} \]
\[ SE(\widehat{\delta}) = \sqrt{\frac{N_{12} + N_{21}}{N_{++}^2}} \]
Each subject produces a 2x2 table with rows being years, columns being political parties.
E.g. a Bush supporter in 2004 who changed votes to Obama would be
voter = matrix(c(0,1,1,0), 2, 2)
rownames(voter) = c(2004, 2008)
colnames(voter) = c('D', 'R')
voter## D R
## 2004 0 1
## 2008 1 0We can thus form a \(2x2x433\) table representing each voter’s data.
Marginal homogeneity \(\equiv\) conditional independence in \(2x2x433\) table.
Let \(N\) be a \(2x2xK\) table.
\(N_{11k} \vert \text{marginals of } N[:,:,k] \sim \text{Hypergeometric}(\dots)\)
Moments:
\[ \begin{aligned} \mu_{11k} &= E[N_{11k} \vert \text{marginals of } N[:,:,k]] \\ &= \frac{N_{1+k}N_{+1k}}{N_{++k}} \\ \sigma^2_{11k} &= \text{Var}(N_{11k} \vert \text{marginals of } N[:,:,k]) \\ &= \frac{N_{1+k}N_{2+k}N_{+1k}{N_{+2k}}}{N^2_{++k}(N_{++k}-1)} \end{aligned} \]
\[ CMH = \frac{[\sum_k (N_{11k} - \mu_{11k})]^2}{\sum_k \sigma^2_{11k}} \]
\[ CMH \overset{k \to \infty}{\to} \chi^2_1 \]
Cochran-Mantel-Haenszel with one stratum per subject (433 subjects) on \(2x2x433\) table.
\(X\): whether vote is Republican or Democrat
\(Y\): year of vote
\(Z\): subject label
Applying CMH: no variation when both votes are R or D: only discordant pairs contribute.
pi_2008 = (175 + 54) / (175 + 16 + 54 + 188) # total is 433
pi_2004 = (175 + 16) / (175 + 16 + 54 + 188)
delta = pi_2004 - pi_2008
SE = sqrt(54 + 16) / (175 + 16 + 54 + 188)
Z = delta / SE
Z## [1] -4.541869We’ve observed 70 discordant pairs (people who changed votes).
Under \(H_0\), \(N_{12} | N_{12} + N_{21} \sim \text{Binomial}(N_{12} + N_{21}, 1/2)\)
Exact p-value can be computed using e.g. pbinom.
## [1] 2.926978e-06## [1] -4.541869## [1] 5.575776e-06 5.853956e-06Data can be represented as \(Y_i=(Y_{i,1}, Y_{i,2}), 1 \leq i \leq 433\) where \(Y_{i,1}\) are votes in 2004 and \(Y_{i,2}\) are votes in 2008 (with, say, \(Y_{i,1}=1\) if vote is D).
Might model these as
\[ \text{logit}(P(Y_{i,j}=1)) = \begin{cases} \alpha & j=1 \\ \alpha + \beta & j=2 \end{cases} \]
Might model these as \[ \text{logit}(P(Y_{i,j}=1)) = \begin{cases} \alpha_i & j=1 \\ \alpha_i + \beta & j=2 \end{cases} \]
This has 434 parameters!
Can’t expect to be easy to fit, or standard MLE theory to work.
Eliminate \(\alpha_i\)’s by conditioning – this is yields a conditional likelihood similar to those we saw when discussing pseudolikelihood.
\[ \begin{aligned} L(\alpha_1, \dots, \alpha_{433}, \beta) &\propto \exp \biggl(\biggl(\sum_i \alpha_i (Y_{i,1}+Y_{i,2})\biggr) \\ & \qquad + \beta \biggl(\sum_i Y_{i,2}\biggr)\biggr) \end{aligned} \]
Sufficient statistic for \(\alpha_i\) is \(Y_{i,1}+Y_{i,2}\)
Conditioning on \((Y_{i,1}+Y_{i,2})_{1 \leq i \leq n}\) eliminates \(\alpha_i\) from the likelihood.
No information when \(Y_{i,1}+Y_{i,2}\) is 0 or 2. Focuses again on discordant pairs…
For discordant pairs, let
\[ Y_k^* = \begin{cases} 1 & Y_{k,2} = 1 \\ 0 & Y_{k,1} = 1 \end{cases} \]
Conditional likelihood \[ \prod_k \left[\frac{e^{\beta}}{1 + e^{\beta}}\right]^{Y^*_k} \left[\frac{1}{1 + e^{\beta}}\right]^{1 - Y^*_k} \]
This is logistic likelihood intercept only (one sample Binomial likelihood) \[ \widehat{\beta} = \text{logit}(\hat{\pi^*}), \qquad \widehat{\pi}^* = \frac{N_{21}}{N_{21} + N_{12}} \]
Suppose the pairs have \(p\)-dimensional features \((X_{i,1}, X_{i,2})\).
Similar conditioning and calculation yields product over discordant pairs \[ \prod_k \left[\frac{e^{X_{k,2}^T\beta}}{e^{X_{k,2}^T\beta} + e^{X_{k,2}^T\beta}}\right]^{Y^*_k} \left[\frac{e^{X_{k,1}^T\beta}}{e^{X_{k,2}^T\beta} + e^{X_{k,2}^T\beta}}\right]^{1-Y^*_k} \]
A standard logistic likelihood with features \(X_k^* = X_{k,2}-X_{k,1}\).
An \(I \times I\) table often has (essentially) same row names as column names.
There are several models of
\[ \pi_{ij} \qquad 1 \leq i, j \leq I \]
or, thinking of log-linear model
\[ \mu_{ij} = \exp(\lambda_{ij}) \qquad 1 \leq i, j \leq I \]
\[ \pi_{i+} = \pi_{+i} \qquad 1 \leq i \leq I \]
or
\[ \mu_{i+} = \mu_{+i} \qquad 1 \leq i \leq I \]
\[ \log \mu_{ij} = \lambda + \lambda_i + \lambda_j + \lambda_{ij} \]
with \(\lambda_{ij} = \lambda_{ji}\)
Note: symmetry implies marginal homogeneity (converse not generally true).
Model fit
\[ \widehat{\mu}_{ij} = \frac{N_{ij} + N_{ji}}{2} \]
\[ X^2 = \sum_{(i,j):i < j} \frac{(N_{ij} - N_{ji})^2}{N_{ij} + N_{ji}} \]
\[ \log \mu_{ij} = \lambda + \lambda^R_i + \lambda^C_j + \lambda_{ij} \]
with \(\lambda_{ij} = \lambda_{ji}\) but \(\lambda^R\) not restricted to agree with \(\lambda^C\).
Model: \[ \log \mu_{ij} = \lambda + \lambda^R_i + \lambda^C_j + \lambda^I_{i} \delta_{ij} \]
Exact independence is \(\lambda^I \equiv 0\).
A lowish dimensional departure from independence.
## [,1] [,2] [,3] [,4]
## [1,] 22 2 2 0
## [2,] 5 7 14 0
## [3,] 0 2 36 0
## [4,] 0 1 17 10carcinoma = data.frame(N=as.numeric(agreement), # Table 11.8
A=factor(c(rep(1,4), rep(2,4), rep(3,4), rep(4,4))),
B=factor(rep(1:4,4)))
indep_fit = glm(N ~ A + B, family=poisson, data=carcinoma)
indep_resid = matrix(resid(indep_fit), 4, 4)
indep_resid # positive on diagonal## [,1] [,2] [,3] [,4]
## [1,] 5.0439297 -0.4140653 -4.2770661 -2.099233
## [2,] -0.4002454 2.2177506 -0.3128405 -2.099233
## [3,] -4.1701116 -1.0460256 2.6798251 -2.537849
## [4,] -3.5796032 -1.2657260 0.1540102 3.676325departure = factor((carcinoma$A == carcinoma$B) * as.numeric(carcinoma$A))
quasi_indep_fit = glm(N ~ A + B + departure, family=poisson, data=carcinoma)
quasi_indep_resid = matrix(resid(quasi_indep_fit), 4, 4)
quasi_indep_resid # residuals are 0 on diagonal## [,1] [,2] [,3] [,4]
## [1,] 5.161914e-08 1.1958607 -0.7480157 -6.169760e-05
## [2,] 1.486609e+00 0.0000000 -0.6636398 -1.395289e-04
## [3,] -1.236669e+00 0.6307957 0.0000000 -7.928944e-05
## [4,] -1.932544e+00 -1.3509531 1.0251702 0.000000e+00## Analysis of Deviance Table
##
## Model 1: N ~ A + B
## Model 2: N ~ A + B + departure
## Resid. Df Resid. Dev Df Deviance
## 1 9 117.957
## 2 5 13.178 4 104.78\[ \kappa = \frac{\sum_i \pi_{ii} - \sum_i \pi_{i+} \pi_{+i}}{1 - \sum_i \pi_{i+} \pi_{+i}} \]
row_margin = apply(agreement, 2, sum) / sum(agreement)
col_margin = apply(agreement, 1, sum) / sum(agreement)
kappa_hat = (sum(diag(agreement)) / sum(agreement) - sum(row_margin * col_margin)) / (1 - sum(row_margin * col_margin))
kappa_hat## [1] 0.4930056baseball = matrix(c(NA,12, 6,10,10,
6,NA, 9,11,13,
12, 9,NA,12, 9,
8, 7, 6,NA,12,
8, 5, 9, 6,NA),
5, 5, byrow=TRUE)
rownames(baseball) = c('Boston',
'New York',
'Tampa Bay',
'Toronto',
'Baltimore')
colnames(baseball) = rownames(baseball)
baseball## Boston New York Tampa Bay Toronto Baltimore
## Boston NA 12 6 10 10
## New York 6 NA 9 11 13
## Tampa Bay 12 9 NA 12 9
## Toronto 8 7 6 NA 12
## Baltimore 8 5 9 6 NALet \(\Pi_{ij}\) denote probability team \(i\) beats team \(j\).
Model
\[ \text{logit}(\Pi_{ij}) = \beta_i - \beta_j \]
##
## Attaching package: 'BradleyTerry2'## The following object is masked _by_ '.GlobalEnv':
##
## baseball## player1 player2 win1 win2
## 1 Boston New York 12 6
## 2 Boston Tampa Bay 6 12
## 3 Boston Toronto 10 8
## 4 Boston Baltimore 10 8
## 5 New York Tampa Bay 9 9
## 6 New York Toronto 11 7
## 7 New York Baltimore 13 5
## 8 Tampa Bay Toronto 12 6
## 9 Tampa Bay Baltimore 9 9
## 10 Toronto Baltimore 12 6BradleyTerry2::BTm(cbind(win1, win2),
player1=player1,
player2=player2,
data=processed_baseball, refcat='Baltimore')## Bradley Terry model fit by glm.fit
##
## Call: BradleyTerry2::BTm(outcome = cbind(win1, win2), player1 = player1,
## player2 = player2, refcat = "Baltimore", data = processed_baseball)
##
## Coefficients:
## ..Boston ..New York ..Tampa Bay ..Toronto
## 0.4539 0.4991 0.6354 0.2287
##
## Degrees of Freedom: 10 Total (i.e. Null); 6 Residual
## Null Deviance: 13.18
## Residual Deviance: 7.703 AIC: 48.66Can break down games into home and away
Let \(\Pi^*_{ij}\) denote probability home team \(i\) beats away team \(j\)
Model \[ \text{logit}(\Pi^*_{ij}) = \alpha + \beta_i - \beta_j. \]