1. Simpson’s Paradox

Simpson’s Paradox is a statistical phenomenon when the association between two variables appears, disappears or reverses when the population is divided into sub-populations. For example, admissions data from UC-Berkeley showed that males are more likely to be admitted than females. However, when breaking the admissions data down by departments, males were less likely to be admitted than females.

Here’s the data that captured admissions of males and females by departments.

import pandas as pd

df = pd.read_csv('./data/berkeley-admission.csv')
df.shape

(24, 4)

df

	department	sex	admission	n
0	A	Male	Admitted	512
1	A	Male	Rejected	313
2	A	Female	Admitted	89
3	A	Female	Rejected	19
4	B	Male	Admitted	353
5	B	Male	Rejected	207
6	B	Female	Admitted	17
7	B	Female	Rejected	8
8	C	Male	Admitted	120
9	C	Male	Rejected	205
10	C	Female	Admitted	202
11	C	Female	Rejected	391
12	D	Male	Admitted	139
13	D	Male	Rejected	279
14	D	Female	Admitted	131
15	D	Female	Rejected	244
16	E	Male	Admitted	53
17	E	Male	Rejected	138
18	E	Female	Admitted	94
19	E	Female	Rejected	299
20	F	Male	Admitted	22
21	F	Male	Rejected	351
22	F	Female	Admitted	24
23	F	Female	Rejected	317

We can just look at the percentage of admitted and rejected students (regardless of department or sex).

df \
    .groupby(['admission']) \
    .agg('sum') \
    .assign(p=lambda d: d['n'] / d['n'].sum())

	n	p
admission
Admitted	1756	0.387895
Rejected	2771	0.612105

When admissions was broken down by sex, females were admitted at 30% and males at 45%. The data, when viewed from this angle, suggested that admission was biased in favor of males.

df \
    .groupby(['sex', 'admission']) \
    .agg('sum') \
    .reset_index() \
    .assign(p=lambda d: d.apply(lambda r: r['n'] / d[d['sex']==r['sex']]['n'].sum(), axis=1))

	sex	admission	n	p
0	Female	Admitted	557	0.303542
1	Female	Rejected	1278	0.696458
2	Male	Admitted	1199	0.445394
3	Male	Rejected	1493	0.554606

If we apply a \(\chi^2\) test of independence, we will see that the sex and admissions are significantly correlated.

from scipy.stats import chi2_contingency

o = df.groupby(['sex', 'admission']).agg('sum').reset_index(drop=True).values.reshape(2, 2)

chi2, p, dof, e = chi2_contingency(o)

print(f'X^2 = {chi2}')
print(f'p = {p}')
print(f'dof = {dof}')

X^2 = 91.87943233270713
p = 9.212151878377876e-22
dof = 1

Here are the observed cell counts.

pd.DataFrame(o,
             columns=pd.MultiIndex.from_tuples([('Admission', 'Admitted'), ('Admission', 'Rejected')]),
             index=pd.MultiIndex.from_tuples([('Sex', 'Female'), ('Sex', 'Male')])
)

		Admission
		Admitted	Rejected
Sex	Female	557	1278
	Male	1199	1493

These are the expected cell counts.

pd.DataFrame(e,
             columns=pd.MultiIndex.from_tuples([('Admission', 'Admitted'), ('Admission', 'Rejected')]),
             index=pd.MultiIndex.from_tuples([('Sex', 'Female'), ('Sex', 'Male')])
)

		Admission
		Admitted	Rejected
Sex	Female	711.787055	1123.212945
	Male	1044.212945	1647.787055

Here comes Simpson’s Paradox. What happens when we partition admission by sex and department? Notice that in only 2 departments do males have higher admission rates than females. Whereas before, it seemed admission was biased in favor of males, now, broken down additionally by department, admission bias is reversed to be in favor of females.

def highlight(s):
    color = 'background-color: rgb(0, 255, 0, 0.18)'
    if s.iloc[0] > s.iloc[2]:
        return [color, '', '', '']
    else:
        return ['', '', color, '']

(df.pivot(index='department', columns=['sex', 'admission'], values='n') / \
    df.pivot(index='department', columns=['sex', 'admission'], values='n').sum()) \
    .style \
    .apply(highlight, axis=1)

sex	Male		Female
admission	Admitted	Rejected	Admitted	Rejected
department
A	0.427023	0.209645	0.159785	0.014867
B	0.294412	0.138647	0.030521	0.006260
C	0.100083	0.137307	0.362657	0.305947
D	0.115930	0.186872	0.235189	0.190923
E	0.044204	0.092431	0.168761	0.233959
F	0.018349	0.235097	0.043088	0.248044

We can also use log-linear models to see what the relationships between admission, sex and department really is. Below, we have

• 1 independence model,

• 3 joint independence models,

• 3 conditional independence models and

• 1 homogenous model.

The homogenous model (AB, AC, BC) has the lowest deviance \(G^2\) at 20.2.

from patsy import dmatrices
import statsmodels.api as sm

dept_levels = ['F', 'A', 'B', 'C', 'D', 'E']

def get_deviance(m, f):
    y, X = dmatrices(f, df, return_type='dataframe')
    r = sm.GLM(y, X, family=sm.families.Poisson()).fit()

    return {'model': m, 'df': r.df_resid, 'deviance': r.deviance}

formulas = {
    '(A, B, C)': 'n ~ admission + sex + C(department, levels=dept_levels)',
    '(A, BC)': 'n ~ admission + sex + C(department, levels=dept_levels) + sex * C(department, levels=dept_levels)',
    '(B, AC)': 'n ~ admission + sex + C(department, levels=dept_levels) + admission * C(department, levels=dept_levels)',
    '(C, AB)': 'n ~ admission + sex + C(department, levels=dept_levels) + admission * sex',
    '(AC, BC)': 'n ~ admission + sex + C(department, levels=dept_levels) + admission * C(department, levels=dept_levels) + sex * C(department, levels=dept_levels)',
    '(AB, BC)': 'n ~ admission + sex + C(department, levels=dept_levels) + admission * sex + sex * C(department, levels=dept_levels)',
    '(AB, AC)': 'n ~ admission + sex + C(department, levels=dept_levels) + admission * sex + admission * C(department, levels=dept_levels)',
    '(AB, AC, BC)': 'n ~ (admission + sex + C(department, levels=dept_levels))**2'
}

result_df = pd.DataFrame([get_deviance(m, f) for m, f in formulas.items()]).set_index('model')
result_df

	df	deviance
model
(A, B, C)	16	2097.116655
(A, BC)	11	876.743972
(B, AC)	11	1242.058607
(C, AB)	15	2003.390913
(AC, BC)	6	21.685924
(AB, BC)	10	783.018230
(AB, AC)	10	1148.332865
(AB, AC, BC)	5	20.217417

If we perform a likehood ratio test of the homogenous model with all other ones, we see that the conditional independence model (AC, BC) does not differ from the homogenous one (p = 0.23). This conditional independence model may be viewed as “simpler” than the homogenous one, and favored. So what? Well, (AC, BC) means that admission and sex are independent given the department, and whatever association we saw was spurious between sex and admission.

from scipy.stats import chi2

def do_likelihood_test(model1, model2='(AB, AC, BC)'):
    chi_sq = result_df.loc[model1].deviance - result_df.loc[model2].deviance
    dof = result_df.loc[model1].df - result_df.loc[model2].df

    return 1 - chi2.cdf(chi_sq, dof)

pd.DataFrame(
    [(m, do_likelihood_test(m)) for m in result_df.index.drop('(AB, AC, BC)')],
    columns=['model', 'p-value']
)

	model	p-value
0	(A, B, C)	0.000000
1	(A, BC)	0.000000
2	(B, AC)	0.000000
3	(C, AB)	0.000000
4	(AC, BC)	0.225581
5	(AB, BC)	0.000000
6	(AB, AC)	0.000000

How about propensity score matching? We can treat sex as the treatment T and admission as the outcome O with department as the feature X.

import itertools

def to_records(r):
    return [(r['department'], r['sex'], r['admission']) for _ in range(r['n'])]

raw_df = pd.DataFrame(
            itertools.chain(*df.apply(to_records, axis=1)),
            columns=['department', 'sex', 'admission']) \
        .assign(
            sex=lambda d: d['sex'].apply(lambda s: 1 if s == 'Male' else 0),
            admission=lambda d: d['admission'].apply(lambda s: 1 if s == 'Admitted' else 0)
        )

raw_df.shape

(4527, 3)

raw_df.head()

	department	sex	admission
0	A	1	1
1	A	1	1
2	A	1	1
3	A	1	1
4	A	1	1

Let’s learn the Inverse Probability Treatment Weight IPTW.

import numpy as np
from sklearn.linear_model import LogisticRegression

dept_levels = ['F', 'A', 'B', 'C', 'D', 'E']

y, X = dmatrices('sex ~ C(department, levels=dept_levels)', raw_df, return_type='dataframe')
X = X.iloc[:,1:]
y = np.ravel(y.values)

model = LogisticRegression(random_state=37)
model.fit(X, y)

LogisticRegression(random_state=37)

def get_propensity(r):
    return 1 / r['propensity_score'] if r['T'] == 1 else 1 / (1 - r['propensity_score'])

psm_df = pd.DataFrame({
            'T': y,
            'propensity_score': model.predict_proba(X)[:,1]
        }) \
        .assign(
            inverse_weight=lambda d: d.apply(get_propensity, axis=1),
            sampling_weight=lambda d: d['inverse_weight'] / d['inverse_weight'].sum()
        )

And then we can use the IPTW to do weighted logistic regression.

y, X = dmatrices('admission ~ sex + C(department, levels=dept_levels) - 1', raw_df, return_type='dataframe')
X = X.iloc[:,1:]
y = np.ravel(y.values)
w = psm_df['inverse_weight']

model = LogisticRegression()
model.fit(X, y, w)

coef = pd.concat([
    pd.Series(model.intercept_, ['intercept']),
    pd.Series(model.coef_[0], X.columns)
])

As you can see, sex (being male) has a negative average treatment effect ATE on admission. It decreases the log odds by 0.22.

coef

intercept                              -2.436449
C(department, levels=dept_levels)[A]    3.491369
C(department, levels=dept_levels)[B]    3.177035
C(department, levels=dept_levels)[C]    1.941869
C(department, levels=dept_levels)[D]    1.879617
C(department, levels=dept_levels)[E]    1.481090
sex                                    -0.218998
dtype: float64

Or, sex (being male), decreases the odds of admission by about 20%.

np.exp(coef)

intercept                                0.087471
C(department, levels=dept_levels)[A]    32.830876
C(department, levels=dept_levels)[B]    23.975554
C(department, levels=dept_levels)[C]     6.971766
C(department, levels=dept_levels)[D]     6.550998
C(department, levels=dept_levels)[E]     4.397738
sex                                      0.803323
dtype: float64