Calculating odds ratio when there is a confounder: contingency table vs logistic regression

For a 2x2 contingency table,

$$\begin{array}{ccc} \hline & Y=0 & Y=1 \\ \hline T=0 & n_{00} & n_{01} \\ \hline T=1 & n_{10} & n_{11} \\ \hline \end{array}$$

calculating the odds ratio is straightforward: $$\displaystyle \text{OR}=\frac{n_{11} / n_{10}}{n_{01}/n_{00}}$$.

This odds ratio is also the same as that one calculated using logistic regression: one can prove that $$\displaystyle e^{\beta_0} = \frac{n_{01}}{n_{00}}$$ and $$\displaystyle e^{\beta_T}=\frac{n_{11} / n_{10}}{n_{01}/n_{00}}$$. For details, please see Appendix 2 below.

However, I am not sure if such an agreement holds when there is a confounder:

$$\begin{array}{ccc} \hline & & Y=0 & Y=1 \\ \hline Z=0 & T=0 & n_{000} & n_{001} \\ \hline & T=1 & n_{010} & n_{011} \\ \hline Z=1 & T=0 & n_{100} & n_{101} \\ \hline & T=1 & n_{110} & n_{111} \\ \hline \end{array}$$

In this case, my thought of calculating the odds ratio using the table is in the following way:

Conditional on $$Z=0$$, $$\displaystyle \text{OR}_0 = \frac{n_{011}/n_{010}}{n_{001}/n_{000}}$$. This can be thought as the "individual treatment effect" for people of $$Z=0$$.

Conditional on $$Z=1$$, $$\displaystyle \text{OR}_1 = \frac{n_{111}/n_{110}}{n_{101}/n_{100}}$$. This can be thought as the "individual treatment effect" for people of $$Z=1$$.

Then the "average treatment effect" is $$\Pr(Z=0) \, \text{OR}_0 + \Pr(Z=1) \, \text{OR}_1$$, where $$\Pr(Z=0)$$ can be estimated as the sum of the first two rows over the sum of all rows, and similarly for $$\Pr(Z=1)$$.

However, such an estimate of odds ratio does not agree with the logistic regression estimate (which by the way has no analytical solution to $$e^{\beta}$$ in this case). I am demonstrating my calculations on the "kidney stone data" at Appendix 1.

My questions:

1. Is my way of calculating the odds ratio correct for the 4x2 contingency table with a confounder?

2. If indeed there is no agreement between the two ways of calculating the odds ratio when there is a confounder, what justifications do we have for using logistic regression to "adjust" or "correct" for confounders and estimate the "true" causal effect of the treatment variable?

3. Is there an intuitive explanation why logistic regression agrees with contingency table at all for the 2x2 table? To me, logistic regression is a parametric model with rigid assumptions while contingency table seems nonparametric with minimal assumptions.

Appendix 1: Kidney stone data

$$\begin{array}{ccc} \hline & & Y=0 & Y=1 \\ \hline Z=0 & T=0 & 6 & 81 \\ \hline & T=1 & 36 & 234 \\ \hline Z=1 & T=0 & 71 & 192 \\ \hline & T=1 & 25 & 55 \\ \hline \end{array}$$

Ignoring the confounder $$Z$$, we have the typical 2x2 table:

$$\begin{array}{cc} \hline & Y = 0 & Y = 1 \\ \hline T=0 & 77 & 273 \\ \hline T=1 & 61 & 289 \\ \hline \end{array}$$

And both ways of calculating the odds ratio give 1.336.

But for the original 4x2 table, my way of calculating the odds ratio gives 0.328, but logistic regression gives 0.700 (which agrees with numerically solving the system of nonlinear equations at the end of the post). (In both cases, the odds ratio decreases from above 1 to below 1, hence the "Simpson's paradox".)

import numpy as np
from sklearn.linear_model import LogisticRegression

X = np.array([(0,0)]*87 + [(0,1)]*270 + [(1,0)]*263 + [(1,1)]*80)
y = np.array([0]*6+[1]*81 + [0]*36+[1]*234 + [0]*71+[1]*192 + [0]*25+[1]*55)

f = LogisticRegression(C=1e9).fit(X, y)

np.exp(f.intercept_), np.exp(f.coef_)
# (array([9.91206217]), array([[0.28349989, 0.69962134]]))


Appendix 2: Derivation

For the 2x2 contingency table: $$\displaystyle \hat{y} = \frac{1}{1+e^{-(\beta_0 + \beta_T T)}}$$

$$\displaystyle \hat{\beta} = \arg \min_{\beta} \hat{R} = \arg \min_{\beta} \sum_i L(y_i, \hat{y}_i) = \arg \min_{\beta} \sum_i \Big(y_i \log \hat{y}_i + (1-y_i) \log (1-\hat{y}_i) \Big) \\= \arg \min_{\beta} \Big( n_{00} \big(1 - \frac{1}{1+e^{-\beta_0}}\big) + n_{01} \frac{1}{1+e^{-\beta_0}} + n_{10} \big(1 - \frac{1}{1+e^{-(\beta_0+\beta_T)}}\big) + n_{11} \frac{1}{1+e^{-(\beta_0+\beta_T)}} \Big)$$

Setting the partial derivative of empirical risk $$\hat{R}$$ with respect to parameters to zero, one has:

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_0} = 0 \Rightarrow \frac{-n_{00}}{1+e^{-\beta_0}} + \frac{n_{01}}{1+e^{\beta_0}} + \frac{ -n_{10}}{1+e^{-(\beta_0+\beta_T)}} + \frac{n_{11}}{1+e^{\beta_0+\beta_T}} = 0$$

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_T} = 0 \Rightarrow \frac{- n_{10}}{1+e^{-(\beta_0+\beta_T)}} + \frac{n_{11}}{1+e^{\beta_0+\beta_T}} = 0$$

Solving them, one has: $$\displaystyle e^{\beta_0} = \frac{n_{01}}{n_{00}}$$ and $$\displaystyle e^{\beta_T}=\frac{n_{11} / n_{10}}{n_{01}/n_{00}}$$.

For the 4x2 contingency table with a confounder: $$\displaystyle \hat{y} = \frac{1}{1+e^{-(\beta_0+\beta_Z Z + \beta_T T)}}$$. Repeating the same procedure one has the following system of nonlinear equations:

Let $$\displaystyle A_1 = \frac{-n_{000}}{1+e^{-\beta_0}} + \frac{n_{001}}{1+e^{\beta_0}}$$,

$$\displaystyle A_2 = \frac{-n_{010}}{1+e^{-(\beta_0+\beta_T)}} + \frac{n_{011}}{1+e^{\beta_0+\beta_T}}$$,

$$\displaystyle A_3 = \frac{-n_{100}}{1+e^{-(\beta_0+\beta_Z)}} +\frac{n_{101}}{1+e^{\beta_0+\beta_Z}}$$,

$$\displaystyle A_4 = \frac{-n_{110}}{1+e^{-(\beta_0+\beta_Z+\beta_T)}} + \frac{n_{111}}{1+e^{\beta_0+\beta_Z+\beta_T}}$$,

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_0} = 0 \Rightarrow A_1 + A_2 + A_3 + A_4 = 0$$

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_Z} = 0 \Rightarrow A_3 + A_4 = 0$$

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_T} = 0 \Rightarrow A_2 + A_4 = 0$$

And the system of equations does not seem to have an analytical solution.

Edits:

@Noah suggested adding an interaction term to the model, and it works.

Logistic regression model: $$\displaystyle \hat{y} = \frac{1}{1+e^{-(\beta_0 + \beta_Z Z + \beta_T T + \beta_{ZT} Z*T)}}$$.

The terms $$A_1$$ to $$A_3$$ are the same for this model, while $$A_4$$ becomes $$\displaystyle \frac{-n_{110}}{1+e^{-(\beta_0+\beta_Z+\beta_T + \beta_{ZT})}} + \frac{n_{111}}{1+e^{\beta_0+\beta_Z+\beta_T+\beta_{ZT}}}$$.

And again, we have:

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_0} = 0 \Rightarrow A_1 + A_2 + A_3 + A_4 = 0$$

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_Z} = 0 \Rightarrow A_3 + A_4 = 0$$

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_T} = 0 \Rightarrow A_2 + A_4 = 0$$

$$\displaystyle \frac{\partial \hat{R}}{\partial \beta_{ZT}} = 0 \Rightarrow A_4 = 0$$

In this case, the system of equations have analytical solutions:

$$\displaystyle e^{\beta_0} = \frac{n_{001}}{n_{000}}, \, e^{\beta_Z} = \frac{n_{101}/n_{100}}{n_{001}/n_{000}}, \, e^{\beta_T} = \frac{n_{011}/n_{010}}{n_{001}/n_{000}}, \, e^{\beta_{ZT}} = \frac{(n_{111}/n_{110}) (n_{001}/n_{000})}{(n_{011}/n_{010}) (n_{101}/n_{100})}$$.

• See my answer here. If you include the T by Z interaction in the logistic regression model, you should be able to reconstruct the conditional odds ratios in the table.
– Noah
Commented Feb 19, 2020 at 22:04
• @Noah Adding a new column of $Z$-by-$T$ interaction to the matrix $X$ gives me an $e^{\beta_T}$ that is the same as $\text{OR}_0$, the odds ratio computed from the table conditional on $Z=0$. Is this what you meant? Commented Feb 19, 2020 at 22:26
• Yes. $e^{\beta_T + \beta_Z}$ should give you the odds ratio computed from the table for $Z=1$.
– Noah
Commented Feb 19, 2020 at 22:41
• @Noah It doesn't seem so though: Using the kidney stone data, $e^{\beta_Z}=0.200$ and $e^{\beta_T}$= 0.481, while $\text{OR}_0=0.481$ and $\text{OR}_1=0.168$. Commented Feb 19, 2020 at 22:47
• That appears to be a quirk of Python which uses 1 as the reference category for $Z$. That's an issue of labeling, not an estimation problem. Also, $e^{\beta_Z}$ is the odds ratio for the confounder at $T=0$, which is verified by the table.
– Noah
Commented Feb 19, 2020 at 23:11

2. Logistic regression makes the assumption that the odds ratio in your two $$Z$$ groups are the same. If so, then it is the most efficient method, in the sense it reaches the Cramer-Rao bound asymptotically and is consistent. If not, then it is "wrong", but still approximates the population causal odds ratio, as you calculated (see some work by Greenland in the book).