# 1. Omitted Variable Bias
### What is Omitted Variable Bias?
Imagine trying to evaluate which factors have an impact on someone's hourly wage,
such as age, education, and work experience. Your formula works reasonably well,
but you forgot to include one important factor: whether the person identifies as
male or female. Because wages differ systematically between men and women in the
real world, the model tasked with evaluating the impact of each factor will be
"confused." It tries to fit *everyone* with a single average pattern, which ends
up being wrong in different directions for each group.
This is called **omitted variable bias**: when an important factor is left out of
a model, the model's predictions become systematically skewed, often in ways that
hurt specific groups.
In this tutorial, we use a synthetic (artificially generated) wage dataset to see
exactly what happens when we **leave sex out** of a wage prediction model.
```{code-cell}
:tags: ["remove_input", "remove_output"]
# Install packages when running in JupyterLite (Pyodide) via Thebe.
# In a regular Python environment micropip does not exist, so the except branch runs instead.
try:
import micropip
await micropip.install(['scikit-learn', 'matplotlib', 'pandas', 'numpy'])
except ImportError:
pass
```
```{code-cell}
# ── Import the tools we need ─────────────────────────────────────────────
import pandas as pd # loading and handling data
import numpy as np # numerical operations
import matplotlib.pyplot as plt # creating charts
from sklearn.model_selection import train_test_split # splitting data into train/test sets
from sklearn.linear_model import LinearRegression # our wage prediction model
from sklearn.metrics import mean_squared_error # measuring prediction errors
from IPython.display import display, HTML
import warnings
warnings.filterwarnings('ignore')
import json as _json
# ── Load the dataset ───────────────────────────────────────────────────────
# wages.csv was generated by generate_wages.py
df = pd.read_csv("wages.csv")
print(f"Loaded {len(df)} workers with {len(df.columns)} variables each.")
df.head()
```
```{code-cell}
# ── Prepare the data for modelling ────────────────────────────────────────
# The columns we will use as inputs to predict wages
feature_cols = ['age', 'education', 'experience', 'hours_per_week', 'tenure', 'sex']
# Human-readable labels used in the interactive widget below
feature_labels = {
'age': 'Age',
'education': 'Education (years)',
'experience': 'Experience (years)',
'hours_per_week': 'Hours/Week',
'tenure': 'Tenure (years)',
'sex': 'Sex',
}
# Split the dataset: 70% for training the model, 30% for testing it
df_train, df_test = train_test_split(df, test_size=0.3, random_state=42)
df_test = df_test.copy() # avoids a pandas warning when we add columns later
# Separate the input columns from the outcome we want to predict
X_train = df_train[feature_cols]
X_test = df_test[feature_cols]
y_train = df_train['wage']
y_test = df_test['wage']
# RMSE (Root Mean Squared Error) = average prediction error in $/hr — lower is better
def rmse(y_true, y_pred):
return np.sqrt(mean_squared_error(y_true, y_pred))
print(f"Training set: {len(df_train)} workers | Test set: {len(df_test)} workers")
```
## Why does sex matter for wage prediction?
Before building any models, let's look at how wages are spread across the two groups.
The chart below is a **histogram**: each bar shows how many people earn wages in a given
range. When the two groups have very different distributions, a model that ignores group
membership will struggle to predict wages accurately for either group.
```{code-cell}
:tags: ["remove_input"]
fig, ax = plt.subplots(figsize=(9, 5))
female_wages = df.loc[df['sex'] == 0, 'wage']
male_wages = df.loc[df['sex'] == 1, 'wage']
ax.hist(female_wages, bins=40, alpha=0.55, color='coral', label='Female', density=True)
ax.hist(male_wages, bins=40, alpha=0.55, color='steelblue', label='Male', density=True)
ax.axvline(female_wages.mean(), color='coral', linestyle='--', linewidth=1.8,
label=f'Female mean: ${female_wages.mean():.2f}/hr')
ax.axvline(male_wages.mean(), color='steelblue', linestyle='--', linewidth=1.8,
label=f'Male mean: ${male_wages.mean():.2f}/hr')
ax.set_xlabel('Hourly Wage ($/hr)')
ax.set_ylabel('Density')
ax.set_title('Wage Distribution by Sex')
ax.legend()
plt.tight_layout()
plt.show()
print(f"Mean wage — Female: ${female_wages.mean():.2f}/hr | Male: ${male_wages.mean():.2f}/hr")
print(f"Gap: ${male_wages.mean() - female_wages.mean():.2f}/hr")
```
The two wage distributions are clearly shifted relative to each other. A model that
cannot tell men from women will apply a single "average" wage formula to everyone
which will systematically over-predict wages for women and under-predict for men.
## Training Two Models
We now train two **linear regression** models on the same data, attempting to find
a linear relationship between the input variables (like age or education) and wages.
- **Model WITH sex**: uses age, education, experience, hours per week, tenure, *and* sex
- **Model WITHOUT sex**: uses everything *except* sex using `.drop(columns=['sex'])`
```{code-cell}
# ── Model WITH sex ─────────────────────────────────────────────────────────
lr_with = LinearRegression()
lr_with.fit(X_train, y_train)
# ── Model WITHOUT sex ──────────────────────────────────────────────────────
# Dropping a column from a DataFrame is straightforward:
lr_without = LinearRegression()
lr_without.fit(X_train.drop(columns=['sex']), y_train)
# ── Store predictions and residuals as new columns in the test table ───────
# This lets us filter by sex easily in the next steps — no index tricks needed.
df_test['pred_with'] = lr_with.predict(X_test)
df_test['pred_without'] = lr_without.predict(X_test.drop(columns=['sex']))
df_test['resid_with'] = df_test['pred_with'] - df_test['wage']
df_test['resid_without'] = df_test['pred_without'] - df_test['wage']
print(f"Model WITH sex: RMSE = {rmse(y_test, df_test['pred_with']):.3f} $/hr")
print(f"Model WITHOUT sex: RMSE = {rmse(y_test, df_test['pred_without']):.3f} $/hr")
```
The overall error is already higher when sex is omitted. But the overall number
hides *who* the model is getting wrong. Let's break the error down by sex group.
```{code-cell}
# Filter the test set by sex
female = df_test[df_test['sex'] == 0]
male = df_test[df_test['sex'] == 1]
female_rmse_with = rmse(female['wage'], female['pred_with'])
male_rmse_with = rmse(male['wage'], male['pred_with'])
female_rmse_without = rmse(female['wage'], female['pred_without'])
male_rmse_without = rmse(male['wage'], male['pred_without'])
print("Prediction error by group (RMSE in $/hr):")
print(f" Model WITH sex Female: {female_rmse_with:.2f} Male: {male_rmse_with:.2f}")
print(f" Model WITHOUT sex Female: {female_rmse_without:.2f} Male: {male_rmse_without:.2f}")
```
The bar chart below makes the group-level differences easy to see at a glance.
Each pair of bars shows the prediction error for women (left) and men (right),
comparing the model that includes sex (blue) with the one that does not (coral).
```{code-cell}
fig, ax = plt.subplots(figsize=(9, 5))
x = np.arange(2)
width = 0.35
bars1 = ax.bar(x - width/2, [female_rmse_with, male_rmse_with], width,
label='With sex', color='steelblue', alpha=0.8)
bars2 = ax.bar(x + width/2, [female_rmse_without, male_rmse_without], width,
label='Without sex', color='coral', alpha=0.8)
ax.set_ylabel('RMSE ($/hr)')
ax.set_title('Prediction Error by Sex:\nControlling vs. Not Controlling for Sex')
ax.set_xticks(x)
ax.set_xticklabels(['Female', 'Male'])
ax.legend()
ax.set_ylim(0, max(female_rmse_without, male_rmse_without) * 1.25)
for bar in list(bars1) + list(bars2):
h = bar.get_height()
ax.annotate(f'{h:.2f}', xy=(bar.get_x() + bar.get_width()/2, h),
xytext=(0, 3), textcoords="offset points", ha='center', va='bottom', fontsize=10)
plt.tight_layout()
plt.show()
```
## Residual Analysis
While the RMSE gives us a general idea of how large the errors are, we still
don't know if the model is over-predicting or under-predicting wages for each
group. To investigate that, we look at **residuals**, the difference between
what the model predicted and what the actual wage was:
> **Residual = Predicted wage − Actual wage**
- A **positive** residual means the model predicted *too high* (over-prediction).
- A **negative** residual means the model predicted *too low* (under-prediction).
- Residuals clustered near **zero** mean the model is unbiased on average.
```{code-cell}
# Plot residual histograms side by side
# 'female' and 'male' already contain the residual columns from the training step
fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharey=True)
for ax, col, title in [
(axes[0], 'resid_with', 'Model WITH sex'),
(axes[1], 'resid_without', 'Model WITHOUT sex'),
]:
ax.hist(female[col], bins=40, alpha=0.55, color='coral',
label=f'Female (mean {female[col].mean():+.2f})', density=True)
ax.hist(male[col], bins=40, alpha=0.55, color='steelblue',
label=f'Male (mean {male[col].mean():+.2f})', density=True)
ax.axvline(female[col].mean(), color='coral', linestyle='--', linewidth=1.8)
ax.axvline(male[col].mean(), color='steelblue', linestyle='--', linewidth=1.8)
ax.axvline(0, color='black', linewidth=1, alpha=0.4) # perfect prediction line
ax.set_xlabel('Residual ($/hr) [predicted − actual]')
ax.set_ylabel('Density')
ax.set_title(title)
ax.legend(fontsize=9)
plt.suptitle('Residual Distributions by Sex', fontsize=13, y=1.02)
plt.tight_layout()
plt.show()
print("Mean residuals:")
print(f" WITH sex Female: {female['resid_with'].mean():+.3f} Male: {male['resid_with'].mean():+.3f}")
print(f" WITHOUT sex Female: {female['resid_without'].mean():+.3f} Male: {male['resid_without'].mean():+.3f}")
```
When sex is included, the average residual for both groups is close to zero,
the model is not systematically wrong for either group. When sex is omitted,
the model cannot account for the wage gap, so it over-predicts for women
(positive residuals) and under-predicts for men (negative residuals) by roughly
$1.7/hr in each direction.
## Interactive Feature Explorer
So far we compared just two configurations: all variables included, versus all
variables minus sex. But what happens when you remove other variables or several
at once?
Use the checkboxes below to choose which variables to include, then click
**Train & Evaluate** to see how the prediction errors shift for women and men.
Try removing education, or tenure, and see if the bias changes direction.
```{code-cell}
:tags: ["remove_input", "remove_output"]
_FCOLS = ['age', 'education', 'experience', 'hours_per_week', 'tenure', 'sex']
_FLBLS = {
'age': 'Age', 'education': 'Education (years)',
'experience': 'Experience (years)', 'hours_per_week': 'Hours/Week',
'tenure': 'Tenure (years)', 'sex': 'Sex',
}
_BINS = np.arange(-15, 16, 1.0) # 30 bins of width 1
_bb_data = {}
for _mask in range(1, 64): # all non-empty subsets of 6 features
_sel = [_FCOLS[_i] for _i in range(6) if _mask & (1 << _i)]
_key = ','.join(_sel)
_lr = LinearRegression()
_lr.fit(df_train[_sel], y_train)
_preds = _lr.predict(df_test[_sel])
_resid = _preds - df_test['wage'].values
_fm = df_test['sex'].values == 0
_mm = df_test['sex'].values == 1
_fr = _resid[_fm]
_mr = _resid[_mm]
_fh, _ = np.histogram(_fr, bins=_BINS, density=True)
_mh, _ = np.histogram(_mr, bins=_BINS, density=True)
_bb_data[_key] = {
'rmse': round(float(rmse(y_test, _preds)), 3),
'f_rmse': round(float(rmse(y_test[_fm], _preds[_fm])), 3),
'm_rmse': round(float(rmse(y_test[_mm], _preds[_mm])), 3),
'f_mean': round(float(_fr.mean()), 3),
'm_mean': round(float(_mr.mean()), 3),
'fh': [round(float(v), 5) for v in _fh],
'mh': [round(float(v), 5) for v in _mh],
}
_BB_BINS_LEFT = [float(v) for v in _BINS[:-1]] # left edge of each bin
```
```{code-cell}
:tags: ["remove_input"]
_js = _json.dumps(_bb_data)
_lbls = _json.dumps(_FLBLS)
_bins = _json.dumps(_BB_BINS_LEFT)
display(HTML(f"""
Select features to include:
"""))
```
## Key Observations
By comparing the two models, we can draw several important conclusions about
how omitting a variable introduces bias:
- **Without sex, the model treats everyone the same**: It applies a single
"one size fits all" wage formula, which over-predicts wages for women and
under-predicts wages for men — by roughly $5/hr in each direction.
- **Residuals show which direction the model is wrong**: RMSE only tells you
*how wrong* the model is. Residual histograms reveal *who* the model
systematically favours or disfavours — a clear sign of omitted variable bias.
- **Including sex fixes the problem**: When the model knows about sex, it can
learn group-specific patterns and its errors shrink roughly by half for each group.
- **Other variables matter too**: Use the interactive tool above to explore what
happens when you remove education, tenure, or other variables — some omissions
affect women more, others affect men more.