# Section 4.5 — Model selection#

This notebook contains the code examples from Section 4.5 Model selection from the No Bullshit Guide to Statistics.

## Notebook setup#

# load Python modules
import os
import numpy as np
import pandas as pd
import seaborn as sns

# Figures setup
import matplotlib.pyplot as plt
plt.clf()  # needed otherwise sns.set_theme doesn't work
from plot_helpers import RCPARAMS
# RCPARAMS.update({'figure.figsize': (10, 3)})   # good for screen
RCPARAMS.update({'figure.figsize': (3, 2.5)})  # good for print
sns.set_theme(
context="paper",
style="whitegrid",
palette="colorblind",
rc=RCPARAMS,
)

%config InlineBackend.figure_format = 'retina'

# Where to store figures
DESTDIR = "figures/lm/modelselection"

from ministats.utils import savefigure

<Figure size 640x480 with 0 Axes>

from scipy.stats import bernoulli
from scipy.stats import norm
import statsmodels.formula.api as smf


## Pattern 1: the fork#

### Example 1: simple fork#

Consider a dataset where there is zero causal effect between $$X$$ and $$Y$$.

np.random.seed(41)

n = 200
ws = norm(0,1).rvs(n)
xs = 2*ws + norm(0,1).rvs(n)
ys = 0*xs + 3*ws + norm(0,1).rvs(n)

df1 = pd.DataFrame({"x": xs, "w": ws, "y": ys})


First let’s fit a native model y ~ x to see what happens.

lm1a = smf.ols("y ~ 1 + x", data=df1).fit()
lm1a.params

Intercept   -0.118596
x            1.268185
dtype: float64

from ministats import plot_reg
plot_reg(lm1a)

filename = os.path.join(DESTDIR, "simple_fork_y_vs_x.pdf")
savefigure(plt.gcf(), filename)

Saved figure to figures/lm/modelselection/simple_fork_y_vs_x.pdf
Saved figure to figures/lm/modelselection/simple_fork_y_vs_x.png


Recall the true strength of the causal association $$X \to Y$$ is zero, so the $$\widehat{\beta}_{\tt{x}}=1.268$$ is a wrong estimate.

Now let’s fit another model that controls for the common cause $$W$$, and thus removes the confounding.

lm1b = smf.ols("y ~ 1 + x + w", data=df1).fit()
lm1b.params

Intercept   -0.089337
x            0.077516
w            2.859363
dtype: float64


The model lm1b correctly recovers causal association close to zero.

### Example 2: student’s memory capacity#

np.random.seed(43)

n = 100
mems = norm(5,1).rvs(n)
efforts = 20 - 3*mems + norm(0,1).rvs(n)
scores = 10*mems + 2*efforts + norm(0,1).rvs(n)

df2 = pd.DataFrame({"mem": mems,
"effort": efforts,
"score": scores})


First let’s fit a model that doesn’t account for the common cause.

lm2a = smf.ols("score ~ 1 + effort", data=df2).fit()
lm2a.params

Intercept    64.275333
effort       -0.840428
dtype: float64


We see a negative effect $$\widehat{\beta}_{\tt{effort}}=-0.84$$, the true effect is $$+2$$.

Now we fit a model that includes mem, i.e., we control for the common cause confounder.

lm2b = smf.ols("score ~ 1 + effort + mem", data=df2).fit()
lm2b.params

Intercept    2.650981
effort       1.889825
mem          9.568862
dtype: float64


The model lm2b correctly recovers causal association $$1.89$$, which is close to the true value $$2$$.

TODO explain

### Example 2R: random assignment of the effort variable#

Suppose students are randomly assigned into a low effort (5h/week) and high effort (15h/week) groups by flipping a coin.

np.random.seed(47)

n = 300
mems = norm(5,1).rvs(n)
coins = bernoulli(p=0.5).rvs(n)
efforts = 5*coins + 15*(1-coins)
scores = 10*mems + 2*efforts + norm(0,1).rvs(n)

df2r = pd.DataFrame({"mem": mems,
"effort": efforts,
"score": scores})


The effect of the random assignment is to decouple effort from mem, thus removing the association, as we can see if we compare the correlations in the original dataset df2 and the randomized dataset df2r.

# non-randomized              # with random assignment
df2.corr()["mem"]["effort"],  df2r.corr()["mem"]["effort"]

(-0.9452404905554658, 0.037005646395469494)


Randomization allows us to recover the correct estimate, even without including the common cause.

lm2r = smf.ols("score ~ 1 + effort", data=df2r).fit()
lm2r.params

Intercept    49.176242
effort        2.077639
dtype: float64


Note we can also including mem in the model, and that doesn’t hurt.

lm2rm = smf.ols("score ~ 1 + effort + mem", data=df2r).fit()
lm2rm.params

Intercept    -0.343982
effort        2.002295
mem          10.063906
dtype: float64


## Pattern 2: the pipe#

### Example 3: simple pipe#

np.random.seed(42)

n = 300
xs = norm(0,1).rvs(n)
ms = xs + norm(0,1).rvs(n)
ys = ms + norm(0,1).rvs(n)

df3 = pd.DataFrame({"x": xs, "m": ms, "y": ys})


First we fit the model y ~ 1 + x that doesn’t include the variable $$M$$, which is the correct thing to do.

lm3a = smf.ols("y ~ 1 + x", data=df3).fit()
lm3a.params

Intercept    0.060272
x            0.922122
dtype: float64

from ministats import plot_reg
plot_reg(lm3a)

ax = plt.gca()
# ax.set_title("Regression plot")
filename = os.path.join(DESTDIR, "simple_pipe_y_vs_x.pdf")
savefigure(plt.gcf(), filename)

Saved figure to figures/lm/modelselection/simple_pipe_y_vs_x.pdf
Saved figure to figures/lm/modelselection/simple_pipe_y_vs_x.png


However, if we were choose to include $$M$$ in the model, we get a completely different result.

lm3b = smf.ols("y ~ 1 + x + m", data=df3).fit()
lm3b.params

Intercept    0.081242
x           -0.005561
m            0.965910
dtype: float64

from ministats import plot_partreg

plot_partreg(lm3b, "x");

# ALT
# from statsmodels.graphics.api import plot_partregress
# with plt.rc_context({"figure.figsize":(3, 2.5)}):
#     plot_partregress("y", "x",  exog_others=["m"], data=df3, obs_labels=False)
#     ax = plt.gca()
#     ax.set_title("Partial regression plot")
#     ax.set_xlabel("e(x | m)")
#     ax.set_ylabel("e(y | m)")
#     filename = os.path.join(DESTDIR, "simple_pipe_partregress_y_vs_x.pdf")
#     savefigure(plt.gcf(), filename)

# # OLD ALT. using plot_lm_partial not a good plot,
# #      since it still shows the trend
# from ministats import plot_lm_partial
# sns.scatterplot(data=df3, x="x", y="y")
# plot_lm_partial(lm3b, "x")
# filename = os.path.join(DESTDIR, "simple_pipe_y_vs_x_given_m.pdf")
# savefigure(plt.gcf(), filename)


### Example 4: student competency as a mediator#

Students score as a function of effort. We assume the final scores for the course was mediated through improvements in competency with the material (comp).

np.random.seed(42)

n = 200
efforts = norm(9,2).rvs(n)
comps = 2*efforts + norm(0,1).rvs(n)
scores = 3*comps + norm(0,1).rvs(n)

df4 = pd.DataFrame({"effort": efforts,
"comp": comps,
"score": scores})


The causal graph in this situation is an instance of the mediator pattern, so the correct modelling decision is not to include the comp variable, which allows us to recover the correct parameter.

lm4a = smf.ols("score ~ 1 + effort", data=df4).fit()
lm4a.params

Intercept   -0.538521
effort       6.079663
dtype: float64


However, if we make the mistake of including the variable comp in the model, we will obtain a zero parameter for the effort, which is a wrong conclusion.

lm4b = smf.ols("score ~ 1 + effort + comp", data=df4).fit()
lm4b.params

Intercept    0.545934
effort      -0.030261
comp         2.979818
dtype: float64


It even appears that effort has a small negative effect.

## Pattern 3: the collider#

### Example 5: simple collider#

np.random.seed(42)

n = 300
xs = norm(0,1).rvs(n)
ys = norm(0,1).rvs(n)
zs = (xs + ys >= 1.7).astype(int)

df5 = pd.DataFrame({"x": xs, "y": ys, "z": zs})


Note y is completely independent of x. And if we were to fit the model y ~ x we get the correct result.

lm5a = smf.ols("y ~ 1 + x", data=df5).fit()
lm5a.params

Intercept   -0.021710
x           -0.039576
dtype: float64

plot_reg(lm5a)

ax = plt.gca()
filename = os.path.join(DESTDIR, "simple_collider_y_vs_x.pdf")
savefigure(plt.gcf(), filename)

Saved figure to figures/lm/modelselection/simple_collider_y_vs_x.pdf
Saved figure to figures/lm/modelselection/simple_collider_y_vs_x.png


Including $$Z$$ in the model incorrectly shows $$X \to Y$$ association.

lm5b = smf.ols("y ~ 1 + x + z", data=df5).fit()
lm5b.params

Intercept   -0.170415
x           -0.244277
z            1.639660
dtype: float64

from ministats import plot_partreg

with plt.rc_context({"figure.figsize":(3, 2.5)}):
plot_partreg(lm5b, "x")

ax = plt.gca()
ax.set_title("Partial regression plot")
# ax.set_xlabel("e(x | z)")
# ax.set_ylabel("e(y | z)")
filename = os.path.join(DESTDIR, "simple_collider_partregress_y_vs_x.pdf")
savefigure(plt.gcf(), filename)

Saved figure to figures/lm/modelselection/simple_collider_partregress_y_vs_x.pdf
Saved figure to figures/lm/modelselection/simple_collider_partregress_y_vs_x.png


If we choose to control for z using stratification,we’ll find there is a negative association, in both sub-samples.

df5z0 = df5[df5["z"]==0]
lm5z0 = smf.ols("y ~ 1 + x", data=df5z0).fit()
lm5z0.params

Intercept   -0.164016
x           -0.197910
dtype: float64

df5z1 = df5[df5["z"]==1]
lm5z1 = smf.ols("y ~ 1 + x", data=df5z1).fit()
lm5z1.params

Intercept    2.032647
x           -0.666737
dtype: float64

sns.regplot(x="x", y="y", data=df5z0, ci=None, label="z=0")
sns.regplot(x="x", y="y", data=df5z1, ci=None, line_kws={"linestyle":"dashed"}, label="z=1")
plt.legend();

ax = plt.gca()
ax.set_title("Regression plots")
filename = os.path.join(DESTDIR, "simple_collider_y_vs_x_groupby_z.pdf")
savefigure(plt.gcf(), filename)

Saved figure to figures/lm/modelselection/simple_collider_y_vs_x_groupby_z.pdf
Saved figure to figures/lm/modelselection/simple_collider_y_vs_x_groupby_z.png


### Example 6: student GPA as a collider#

np.random.seed(46)

n = 300
efforts = norm(9,2).rvs(n)
scores = 5*efforts + norm(10,10).rvs(n)
gpas = 0.1*efforts + 0.02*scores + norm(0.6,0.3).rvs(n)

df6 = pd.DataFrame({"effort": efforts,
"score": scores,
"gpa": gpas})


When we don’t adjust for the collider gpa, we obtain the correct effect of effort on scores.

lm6a = smf.ols("score ~ 1 + effort", data=df6).fit()
lm6a.params

Intercept    12.747478
effort        4.717006
dtype: float64


But with adjustment for gpa reduces the effect significantly.

lm6b = smf.ols("score ~ 1 + effort + gpa", data=df6).fit()
lm6b.params

Intercept     1.464950
effort        1.177600
gpa          16.624215
dtype: float64


### Example 6SB: selection bias#

Suppose the study was conducted only on students who obtained an award, which is given to students with GPA 3.4 or greater.

df6sb = df6[df6["gpa"]>=3.4]
lm6sb = smf.ols("score ~ 1 + effort", data=df6sb).fit()
lm6sb.params

Intercept    78.588495
effort       -0.126503
dtype: float64


A negative causal association appears, which is very misleading.

## Complete example#

scary conclusion: there is approx. 5% reduction in lung capacity for a typical 15 year old who smokes

## Discussion#

### Metrics-based variable selection#

# load the doctors data set


Fit the short model

formula2 = "score ~ 1 + alc + weed + exrc"
lm2 = smf.ols(formula2, data=doctors).fit()
# lm2.params


Fit long model with useful loc variable.

formula3 = "score ~ 1 + alc + weed + exrc + C(loc)"
lm3 = smf.ols(formula3, data=doctors).fit()
# lm3.params


Fit long model with useless variable permit.

formula2p = "score ~ 1 + alc + weed + exrc + permit"
lm2p = smf.ols(formula2p, data=doctors).fit()
# lm2p.params


#### Comparing metrics#

lm2.rsquared_adj, lm3.rsquared_adj, lm2p.rsquared_adj

(0.8390497506713147, 0.8506062566189324, 0.8407115273797523)

lm2.aic, lm3.aic, lm2p.aic

(1103.2518084235273, 1092.5985552344712, 1102.6030626936558)


The value lm2p.aic = 1102.60 is smaller than lm2.aic, but not not that different.

lm2.bic, lm3.bic, lm2p.bic

(1115.4512324525256, 1107.8478352707189, 1117.8523427299035)


#### F-test for the submodel#

F, p, _ = lm3.compare_f_test(lm2)
F, p

(12.758115596295502, 0.00047598123084921035)


The $$p$$-value is smaller than $$0.05$$, so we conclude that adding the variable loc improves the model.

F, p, _ = lm2p.compare_f_test(lm2)
F, p

(2.5857397307381698, 0.10991892492566509)


The $$p$$-value is greater than $$0.05$$, so we conclude that adding the variable permit doesn’t improve the model.

# CUT MATERIAL#

## Bonus examples#

### Bonus example 1#

np.random.seed(1897)
n = 1000
iqs = norm(100,15).rvs(n)
ltimes = 200 - iqs + norm(0,1).rvs(n)
tscores = 0.5*iqs + 0.1*ltimes + norm(0,1).rvs(n)

bdf2 = pd.DataFrame({"iq":iqs,
"ltime": ltimes,
"tscore": tscores})

lm2a = smf.ols("tscore ~ 1 + ltime", data=bdf2).fit()
lm2a.params

Intercept    99.602100
ltime        -0.395873
dtype: float64

lm2b = smf.ols("tscore ~ 1 + ltime + iq", data=bdf2).fit()
lm2b.params

Intercept    3.580677
ltime        0.081681
iq           0.482373
dtype: float64


### Bonus example 2#

np.random.seed(1896)
n = 1000
learns = norm(0,1).rvs(n)
knows = 5*learns + norm(0,1).rvs(n)
undstds = 3*knows + norm(0,1).rvs(n)

bdf1 = pd.DataFrame({"learn":learns,
"know": knows,
"undstd": undstds})

blm1a = smf.ols("undstd ~ 1 + learn", data=bdf1).fit()
blm1a.params

Intercept    -0.045587
learn        14.890022
dtype: float64

blm1b = smf.ols("undstd ~ 1 + learn + know", data=bdf1).fit()
blm1b.params

Intercept   -0.036520
learn        0.130609
know         2.975806
dtype: float64


### Bonus example 3#

np.random.seed(42)
n = 1000

ntwrks = norm(0,1).rvs(n)
comps = norm(0,1).rvs(n)
boths = ((ntwrks > 1) | (comps > 1))
lucks = bernoulli(0.05).rvs(n)
proms = (1 - lucks)*boths + lucks*(1 - boths)

bdf3 = pd.DataFrame({"ntwrk": ntwrks,
"comp": comps,
"prom": proms})


Without adjusting for the collider prom, there is almost no effect of the network ability on competence.

blm3a = smf.ols("comp ~ 1 + ntwrk", data=bdf3).fit()
blm3a.params

Intercept    0.071632
ntwrk       -0.041152
dtype: float64


But with adjustment for prom, there seems to be a negative effect.

blm3b = smf.ols("comp ~ 1 + ntwrk + prom", data=bdf3).fit()
blm3b.params

Intercept   -0.290081
ntwrk       -0.239975
prom         1.087964
dtype: float64


The false negative effect can also appear from sampling bias, e.g. if we restrict our analysis only to people who were promoted.

df3proms = bdf3[bdf3["prom"]==1]
blm3c = smf.ols("comp ~ ntwrk", data=df3proms).fit()
blm3c.params

Intercept    0.898530
ntwrk       -0.426244
dtype: float64


### Bonus example 4: benefits of random assignment#

np.random.seed(1896)
n = 1000

iqs = norm(100,15).rvs(n)
groups = bernoulli(p=0.5).rvs(n)
ltimes = 80*groups + 120*(1-groups)
tscores = 0.5*iqs + 0.1*ltimes + norm(0,1).rvs(n)

bdf4 = pd.DataFrame({"iq":iqs,
"ltime": ltimes,
"tscore": tscores})

# non-randomized            # random assignment
bdf2.corr()["iq"]["ltime"],  bdf4.corr()["iq"]["ltime"]

(-0.9979264589333364, -0.020129851374243963)

bdf4.groupby("ltime")["iq"].mean()

ltime
80     99.980423
120    99.358152
Name: iq, dtype: float64

blm4a = smf.ols("tscore ~ 1 + ltime", data=bdf4).fit()
blm4a.params

Intercept    50.688293
ltime         0.091233
dtype: float64

blm4b = smf.ols("tscore ~ 1 + ltime + iq", data=bdf4).fit()
blm4b.params

Intercept   -0.303676
ltime        0.099069
iq           0.503749
dtype: float64

blm4c = smf.ols("tscore ~ 1 + iq", data=bdf4).fit()
blm4c.params

Intercept    9.896117
iq           0.501168
dtype: float64


#### Side-quest (causal influence of mem)#

Suppose we’re interested in … mem wasn’t randomized, but because we de

blm8c = smf.ols("score ~ 1 + mem", data=df2r).fit()
blm8c.params

Intercept    17.382804
mem          10.430159
dtype: float64