Worked Example: Student Grade Prediction

Step 1 — Frame the problem, with the affected party named

Target: the student's final module mark, on a 0–100 scale.

Horizon: predicted at the end of week four, using only information available by then.

Decision supported: which students are invited to an optional study-skills programme. The invitation is not compulsory and the programme is not recorded on a transcript.

Cost of error: asymmetric in a specific way. A false positive (a student flagged who would have been fine) costs a wasted invitation and a small amount of the student's time. A false negative (a student not flagged who then struggles) costs the student access to support at the moment it would have helped. Both matter, but the second is worse.

Who is affected: every first-year student in the cohort. The students are not customers; they are the subjects of a decision made about them by an institution that holds power over their education.

Step 2 — Audit the features, with proxy variables in view

available_columns = [
    # prior attainment
    "prior_qualification_points",   # UCAS-style points
    "entry_maths_grade",

    # engagement (weeks 1–4)
    "vle_logins_weeks_1_4",
    "assignment_1_mark",
    "attendance_pct_weeks_1_4",

    # context
    "term_time_work_hours",
    "caring_responsibilities",
    "commute_minutes",
    "home_postcode_area",           # reduced to outward code
    "secondary_school_type",        # state / independent / FE college

    # demographics
    "age_at_entry",
    "gender",
    "ethnicity",
    "disability_status",
    "eal_flag",                     # English as an additional language

    # NOT available at week 4
    "final_mark",                   # the target
    "exam_board_flag",              # set at the end of the year
]

Applying the leakage test from Tutorial 7 removes final_mark (the target) and exam_board_flag (set at the end of the year). That part is straightforward.

The harder decisions concern the context features. None of them is a protected attribute. All of them are available at week four. But three of them — term_time_work_hours, home_postcode_area and secondary_school_type — are strongly associated with socioeconomic background, and one of them, caring_responsibilities, is associated with gender. They are proxy variables.

Two features are excluded on fairness grounds, and four protected attributes are retained for analysis but excluded from the model. Those four retained columns are essential — a model whose fairness cannot be measured is a model whose fairness cannot be claimed.

Step 3 — Split, and check representativeness

X_train, X_temp, y_train, y_temp = train_test_split(
    X, y, test_size=0.40, random_state=42, stratify=df["ethnicity"]
)
X_val, X_test, y_val, y_test = train_test_split(
    X_temp, y_temp, test_size=0.50, random_state=42,
    stratify=X_temp["ethnicity"]
)

# Keep the sensitive attributes aligned with the splits for later analysis
sensitive_train = X_train[PROTECTED].copy()
sensitive_val   = X_val[PROTECTED].copy()

# Then drop them from the model features
X_train_model = X_train.drop(columns=PROTECTED)
X_val_model   = X_val.drop(columns=PROTECTED)

Two details here are new.

• Stratified splitting. With a smaller dataset, an unstratified random split can produce a validation set with almost no students from a minority group, which makes subgroup metrics meaningless. Stratifying on the sensitive attribute preserves the group proportions across splits.
• Sensitive attributes are separated before dropping. They are removed from the model's feature set, but kept alongside so that the analysis in Step 7 is possible. If you delete them from the dataframe entirely at the start, you cannot check the model afterwards.

It is also worth checking whether the training data resembles the population the model will serve. If the cohort is 62% female and the training data is 71%, the model will be less well calibrated for male students — not because of any modelling error, but because there is less data to learn from.

Step 4 — Explore, and look for the structure that will drive the fairness result

print(y_train.describe())
print(f"skewness: {y_train.skew():.2f}")

# Distribution of the sensitive attributes in the training set
for col in PROTECTED:
    print(f"\n{col}:")
    print(X_train[col].value_counts(normalize=True).round(3))

# Do the proxy features differ across groups?
print("\nmean term_time_work_hours by ethnicity:")
print(X_train.groupby("ethnicity")["term_time_work_hours"].mean().round(1))

print("\nmean secondary_school_type == 'independent' by ethnicity:")
print((X_train.groupby("ethnicity")["secondary_school_type"]
       .apply(lambda s: (s == "independent").mean())).round(3))

final mark:  mean 63.4   std 13.8   min 21   max 96
skewness: -0.21

ethnicity:
  group_a   0.612
  group_b   0.198
  group_c   0.112
  group_d   0.078

mean term_time_work_hours by ethnicity:
  group_a   8.2
  group_b   9.1
  group_c  13.4
  group_d  15.7

mean secondary_school_type == 'independent' by ethnicity:
  group_a   0.241
  group_b   0.152
  group_c   0.071
  group_d   0.042

This is the structure that will produce a fairness problem. Students in group D work nearly twice as many hours as students in group A, and are six times less likely to have attended an independent school. Both features are in the model. Both are predictive of a lower mark. Neither is a protected attribute.

The model has no idea that group D exists. It is about to learn that working long hours and attending a state school predict lower marks — and it will then apply that learning uniformly, which means it will flag group D students more often than group A students, for reasons that have nothing to do with their individual ability.

Step 5 — Build the pipeline and fit

numeric_cols = [
    "prior_qualification_points", "entry_maths_grade",
    "vle_logins_weeks_1_4", "assignment_1_mark",
    "attendance_pct_weeks_1_4", "term_time_work_hours",
    "commute_minutes", "age_at_entry",
]
categorical_cols = [
    "secondary_school_type", "caring_responsibilities",
]

numeric_pipe = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("scale",  StandardScaler()),
])
categorical_pipe = Pipeline([
    ("impute", SimpleImputer(strategy="most_frequent")),
    ("onehot", OneHotEncoder(handle_unknown="ignore", drop="first")),
])

preprocess = ColumnTransformer([
    ("num", numeric_pipe,     numeric_cols),
    ("cat", categorical_pipe, categorical_cols),
])

model = Pipeline([
    ("preprocess", preprocess),
    ("regressor",  LinearRegression()),
])

model.fit(X_train_model, y_train)

Nothing here is new. The pipeline is structurally identical to Tutorial 7's, with a smaller feature set.

Step 6 — Evaluate overall, then per subgroup

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

def group_metrics(y_true, y_pred, label):
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    mae  = mean_absolute_error(y_true, y_pred)
    r2   = r2_score(y_true, y_pred)
    bias = (y_true - y_pred).mean()          # mean residual — systematic error
    n    = len(y_true)
    print(f"{label:<18} n={n:>4}  RMSE {rmse:>5.2f}  MAE {mae:>5.2f}  "
          f"R² {r2:>6.3f}  bias {bias:>+6.2f}")

y_pred_val = model.predict(X_val_model)

print("=== OVERALL ===")
group_metrics(y_val, y_pred_val, "all students")

print("\n=== BY ETHNICITY ===")
for group in sorted(sensitive_val["ethnicity"].unique()):
    mask = sensitive_val["ethnicity"] == group
    group_metrics(y_val[mask], y_pred_val[mask], f"  {group}")

print("\n=== BY DISABILITY STATUS ===")
for group in sorted(sensitive_val["disability_status"].unique()):
    mask = sensitive_val["disability_status"] == group
    group_metrics(y_val[mask], y_pred_val[mask], f"  {group}")

=== OVERALL ===
all students       n= 180  RMSE  9.42  MAE  7.51  R²  0.714  bias  +0.00

=== BY ETHNICITY ===
  group_a          n= 110  RMSE  8.88  MAE  7.10  R²  0.742  bias  +0.62
  group_b          n=  36  RMSE  9.71  MAE  7.68  R²  0.702  bias  −0.51
  group_c          n=  20  RMSE 10.94  MAE  8.72  R²  0.618  bias  −1.84
  group_d          n=  14  RMSE 12.31  MAE  9.96  R²  0.412  bias  −6.21

=== BY DISABILITY STATUS ===
  no               n= 152  RMSE  9.28  MAE  7.42  R²  0.722  bias  +0.18
  yes              n=  28  RMSE 10.02  MAE  7.98  R²  0.668  bias  −0.94

The overall model looks reasonable: RMSE 9.42 marks, R² 0.71, no overall bias. A report that stopped here would be misleading in three specific ways.

• R² ranges from 0.74 to 0.41 across ethnic groups. The model explains less than half the variance for group D. It understands that group substantially less well than the others.
• Group D has a mean residual of −6.21 marks. The model systematically under-predicts group D — it predicts marks that are on average six marks lower than the students actually achieve. This is not random error; it is a consistent directional error, and it is the most serious finding in the table.
• Group D has the smallest sample. Fourteen students in the validation set. The estimates for this group are themselves uncertain, and a responsible analysis would report confidence intervals and note the small n rather than treating the numbers as precise.

Step 7 — Trace the bias to its source

A subgroup difference is a symptom. To act on it, the cause has to be identified. Two checks:

# 1. How does the model distribute support across groups?
threshold = np.percentile(y_pred_val, 10)     # bottom 10% of predictions
flagged = y_pred_val <= threshold

print(f"threshold: {threshold:.1f} marks\n")
print("group   n    flagged   rate    actual_bottom10   over_flag")
for group in sorted(sensitive_val["ethnicity"].unique()):
    mask   = (sensitive_val["ethnicity"] == group).values
    n_grp  = mask.sum()
    n_flag = (flagged & mask).sum()
    # what fraction of this group is genuinely in the bottom 10% of actual marks?
    actual = (y_val <= np.percentile(y_val, 10)).values
    actual_rate = (actual & mask).sum() / n_grp
    print(f"{group:<7} {n_grp:>3}   {n_flag:>4}    "
          f"{n_flag/n_grp:>5.1%}   {actual_rate:>5.1%}          "
          f"{n_flag/n_grp - actual_rate:>+6.1%}")

# 2. Refit without the flagged proxy features and compare subgroup bias
proxy_cols = ["term_time_work_hours", "secondary_school_type"]
model_noproxy = Pipeline([
    ("preprocess", ColumnTransformer([
        ("num", numeric_pipe,
         [c for c in numeric_cols if c not in proxy_cols]),
        ("cat", categorical_pipe,
         [c for c in categorical_cols if c not in proxy_cols]),
    ])),
    ("regressor", LinearRegression()),
]).fit(X_train_model, y_train)

y_pred_np = model_noproxy.predict(X_val_model)

print("\n=== group D bias, with and without proxy features ===")
for label, preds in [("with proxies", y_pred_val), ("without proxies", y_pred_np)]:
    mask = (sensitive_val["ethnicity"] == "group_d").values
    bias = (y_val[mask] - preds[mask]).mean()
    r2   = r2_score(y_val[mask], preds[mask])
    print(f"{label:<18} bias {bias:>+6.2f}   R² {r2:>6.3f}")

threshold: 50.5 marks

group   n    flagged   rate    actual_bottom10   over_flag
group_a 110      6       5.5%      4.5%            +1.0%
group_b  36      4      11.1%      8.3%            +2.8%
group_c  20      4      20.0%     15.0%            +5.0%
group_d  14      4      28.6%     21.4%            +7.2%

=== group D bias, with and without proxy features ===
with proxies       bias  -6.21   R²  0.412
without proxies    bias  -2.08   R²  0.298

Two findings, and they point in opposite directions.

• The support is distributed unequally. Group D students are flagged at 28.6% versus 5.5% for group A. Part of that is genuine need — 21.4% of group D are in the actual bottom 10%, compared with 4.5% of group A. But the over-flagging is also unequal: +7.2 percentage points for group D, +1.0 for group A.
• Removing the proxy features reduces the bias but also the accuracy. Group D's mean residual shrinks from −6.21 to −2.08 — a substantial improvement. But R² for group D falls from 0.41 to 0.30, and the overall model's R² falls as well. The model is now fairer and less useful.

This is the trade-off from Section 3.2, made concrete with numbers. There is no version of this model that is both as accurate and as fair as the alternatives. Which one to deploy depends on which cost is judged more serious — and that is not a statistical question.

Step 8 — Read the coefficients, with the correlation structure in mind

regressor = model.named_steps["regressor"]
features  = model.named_steps["preprocess"].get_feature_names_out()

coef_df = pd.DataFrame({
    "feature": features,
    "coef":    regressor.coef_,
})

# Standardised features: coefficients are per standard deviation
# Convert to per-original-unit by dividing by the feature std dev
stds = X_train_model[numeric_cols].std()
for name, sd in stds.items():
    key = f"num__{name}"
    coef_df.loc[coef_df["feature"] == key, "coef"] /= sd

print(coef_df.sort_values("coef", key=abs, ascending=False).to_string(index=False))

# How stable are these? Refit on 20 bootstrap samples and look at the spread.
boots = []
for seed in range(20):
    idx = np.random.RandomState(seed).choice(len(X_train_model), len(X_train_model), replace=True)
    m = clone(model).fit(X_train_model.iloc[idx], y_train.iloc[idx])
    boots.append(m.named_steps["regressor"].coef_)
boots = np.array(boots)

coef_df["std_over_boot"] = boots.std(axis=0)
coef_df["stability"] = coef_df["coef"].abs() / coef_df["std_over_boot"].replace(0, np.nan)

print("\ncoefficient stability (|coef| / bootstrap sd):")
print(coef_df[["feature", "coef", "std_over_boot", "stability"]]
      .sort_values("stability").to_string(index=False))

              feature   coef
num__assignment_1_mark  0.412
num__attendance_pct_weeks_1_4  0.287
num__prior_qualification_points  0.196
num__vle_logins_weeks_1_4  0.121
num__term_time_work_hours  -0.108
cat__secondary_school_type_independent  0.094
num__entry_maths_grade  0.063
num__commute_minutes  -0.041
num__age_at_entry  0.019

coefficient stability (|coef| / bootstrap sd):
              feature   coef   sd_over_boot   stability
num__age_at_entry  0.019        0.031          0.61
num__commute_minutes  -0.041        0.038          1.08
num__entry_maths_grade  0.063        0.041          1.54
cat__secondary_school_type_independent  0.094        0.052          1.81
num__term_time_work_hours  -0.108        0.044          2.45
num__vle_logins_weeks_1_4  0.121        0.039          3.10
num__prior_qualification_points  0.196        0.036          5.44
num__attendance_pct_weeks_1_4  0.287        0.029          9.90
num__assignment_1_mark  0.412        0.022         18.73

Three coefficients are stable enough to report with confidence: assignment 1 mark, attendance, and prior qualification points. They are the core of the model and their signs and magnitudes are robust.

Two coefficients are not stable: age_at_entry and commute_minutes. Their bootstrap standard deviations are comparable to their point estimates, meaning that a different sample of the same population could easily produce a coefficient of the opposite sign. They should not appear in any report, and they should not be used to justify a decision. If the module leader asks "does commuting hurt marks?", the honest answer is that this dataset cannot say.

Step 9 — Evaluate on the test set, once

y_pred_test = model.predict(X_test_model)

print("=== TEST SET ===")
group_metrics(y_test, y_pred_test, "all students")

print("\nby ethnicity:")
for group in sorted(X_test["ethnicity"].unique()):
    mask = (X_test["ethnicity"] == group).values
    if mask.sum() >= 5:
        group_metrics(y_test[mask], y_pred_test[mask], f"  {group}")

=== TEST SET ===
all students       n= 180  RMSE  9.58  MAE  7.63  R²  0.701  bias  +0.11

by ethnicity:
  group_a          n= 110  RMSE  8.95  MAE  7.18  R²  0.729  bias  +0.44
  group_b          n=  36  RMSE  9.88  MAE  7.79  R²  0.688  bias  −0.62
  group_c          n=  20  RMSE 11.20  MAE  8.91  R²  0.601  bias  −2.01
  group_d          n=  14  RMSE 12.64  MAE 10.12  R²  0.398  bias  −6.48

The test results confirm the validation findings, including the group D bias. The pattern is not an artefact of the validation split — it is a genuine property of the model.

Step 10 — Write the model card