Lab 4 — Building a price predictor: penalties, pipelines, cross-validation

Author

Hicham Janati — FSML

We use similar data as in Lab 3 with a different goal: a real estate agency wants to predict the price of the next house on the market.

We use a larger version of the Ames dataset. The agency gives you 500 sales, and we keep 100 of them aside to test our models.

Boxes as always: 🔮 predict before running (no AI) · 🧠 explain in your own words · ✍️ short task, boilerplate given · 🤖 a prompt that sends you to find something the lab did not hand you.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.preprocessing import StandardScaler, PolynomialFeatures, OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV, train_test_split

BLUE, RED, GREEN, ORANGE, GRAY = "#1f6fb2", "#c92a1a", "#1c7a43", "#e8871a", "#8a8a8a"
plt.rcParams["figure.figsize"] = (6.5, 3.6)
plt.rcParams["axes.spines.top"] = plt.rcParams["axes.spines.right"] = False

scores = {}          # every test RMSE of the lab lands here, compared at the end

The error metric of the whole lab is the root mean squared error of predictions \(\hat y_i\) against the true values \(y_i\):

\[\text{RMSE} \;=\; \sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(y_i - \hat y_i\right)^2}.\]

✍️ Your turn. Write the function rmse(y_true, y_pred) that returns the RMSE of two arrays.

# TODO: define rmse(y_true, y_pred)
    # TODO: return the RMSE

print(rmse(np.array([1, 2, 3]), np.array([1, 2, 5])))      # should print 1.1547...

Part 1 — What the numbers alone can do

df = pd.read_csv("https://hichamjanati.github.io/media/teaching/SMV/ames_prices.csv")
print(df.shape)
df.columns
df.dtypes

🧠 Explain. What does the dtype of a column tell you? What kind of values do the object columns hold here?

✍️ Your answer:

…

select_dtypes keeps the columns of the given types. SalePrice is the response, not a feature, so it leaves the numeric list:

num_cols = df.select_dtypes("number").columns.drop("SalePrice").tolist()
cat_cols = df.select_dtypes("object").columns.tolist()
print(f"{len(num_cols)} numeric features · {len(cat_cols)} categorical features")

The response is the log price, as in Lab 3. train_test_split shuffles the houses and cuts them into a training part and a test part of 100 houses, which we only use to evaluate.

y = np.log(df["SalePrice"])
X = df.drop(columns="SalePrice")

df_train, df_test, y_train, y_test = train_test_split(X, y, test_size=100, random_state=0)
print(f"train {len(df_train)} houses · test {len(df_test)} houses")

✍️ Your turn. In this part we only use the numeric features. Build X_train and X_test: the training and the test frames restricted to num_cols.

# TODO: X_train, X_test = ...
print(X_train.shape, X_test.shape)

1.1 A first model: every number, nothing else

Every scikit-learn model works the same way: create the object, then .fit(X, y) learns its parameters and .predict(X) returns predictions. LinearRegression is ordinary least squares (OLS), as in Lab 3.

ols = LinearRegression()
ols.fit(X_train, y_train)

print(f"intercept {ols.intercept_:.3f}")
pd.Series(ols.coef_, index=num_cols).head()
ols.predict(X_train)[:5]

✍️ Your turn. Compute the RMSE of ols on the training rows and on the test rows.

# TODO: rmse_train = ...
# TODO: rmse_test = ...

scores["numbers only, OLS"] = rmse_test
print(f"train RMSE {rmse_train:.3f}   test RMSE {rmse_test:.3f}")

🧠 Explain. Your test RMSE is about 0.130, measured in log price. How do you interpret this number?

✍️ Your answer:

…

🧠 Explain. The train and test RMSE are close together. What does that tell you about what limits this model?

✍️ Your answer:

…

1.2 More complexity: products of columns

A linear model in these columns says a square meter is worth the same in a mansion and in a shed, and that a garage is worth the same with or without a basement. Interactions relax exactly that: for every pair of columns we add their product, so the effect of one may depend on the level of the other. No powers, just products.

inter = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
I_train = inter.fit_transform(X_train)
I_test = inter.transform(X_test)
print(f"{I_train.shape[1]} columns for {I_train.shape[0]} training houses")

🔮 Predict. There are now more columns than houses. Before running: is β̂ still the unique minimizer of the squared error? What will the training and the test RMSE look like?

✍️ Your answer:

…

✍️ Your turn. Fit a LinearRegression on the interaction columns and compute its RMSE on the training and on the test rows.

# TODO: ols_inter = ...
# TODO: rmse_train_inter = ...
# TODO: rmse_test_inter = ...

scores["interactions, OLS"] = rmse_test_inter
print(f"train RMSE {rmse_train_inter:.3f}   test RMSE {rmse_test_inter:.3f}")

🧠 Explain. The training RMSE fell and the test RMSE rose sharply. Explain how a model that fits the training houses better than the previous one can do so much worse on new ones, and what having more columns than rows did to the estimator.

✍️ Your answer:

…

1.3 A penalty, and why it needs standardized columns

The lecture’s answer is not to delete columns but to shrink them:

\[\hat\beta_\lambda \;=\; \operatorname*{argmin}_{\beta}\; \frac{1}{n}\|\mathbf{y} - \mathbb{X}\beta\|^2 \;+\; \lambda\|\beta\|^2 .\]

With \(\lambda > 0\) the problem has a unique solution again, even with more columns than rows. But the penalty adds up \(\beta_1^2 + \dots + \beta_p^2\), and the size of a coefficient depends on the unit of its column: record the lot area in hectares instead of m² and its coefficient becomes 10 000 times larger, so the penalty hits it much harder although nothing about the houses changed. A penalty on the coefficients only treats the columns alike once they are on a common scale: they must be standardized.

StandardScaler does that: it centers each column and divides it by its standard deviation. It is a transformer: .fit(X) learns the means and standard deviations, .transform(X) applies them.

scaler = StandardScaler()
scaler.fit(I_train)
S_train = scaler.transform(I_train)

print(f"means {S_train.mean(axis=0)[:3].round(2)}   "
      f"standard deviations {S_train.std(axis=0)[:3].round(2)}")

🧠 Explain. Which means and standard deviations should be used to standardize the test rows, and why?

✍️ Your answer:

…

S_test = scaler.transform(I_test)

Ridge(alpha=10) is the penalized least squares above, with λ called alpha. It has the same .fit / .predict API as LinearRegression.

✍️ Your turn. Fit Ridge(alpha=10) on the standardized interaction columns and compute its RMSE on the training and on the test rows.

# TODO: ridge = ...
# TODO: rmse_train_ridge = ...
# TODO: rmse_test_ridge = ...

print(f"ridge λ=10: train RMSE {rmse_train_ridge:.3f}   test RMSE {rmse_test_ridge:.3f}")

A Pipeline chains the scaler and the model into one estimator: .fit fits the scaler on the rows it is given, transforms them and fits the model; .predict transforms with the stored scaler and predicts. In the next section the whole thing is refitted on dozens of subsets of the training rows, and the scaler has to be recomputed inside each one.

pipe = Pipeline([("scale", StandardScaler()), ("model", Ridge(alpha=10))])
pipe.fit(I_train, y_train)
print(f"pipeline, ridge λ=10: test RMSE {rmse(y_test, pipe.predict(I_test)):.3f}")

1.4 Choosing λ without touching the test set

🔮 Predict. λ = 10 was picked out of the air. The quickest way to do better: try 25 values and keep the one with the smallest test RMSE. The lecture forbids it. Write down what that number would stop measuring, and why the procedure gets worse the longer the list of candidates.

✍️ Your answer:

…

The correct protocol splits the training set again: fit on four fifths, score on the fifth, rotate, average. This is cross-validation, and GridSearchCV runs it for every value of a grid. Start from the estimator to tune:

ridge_pipe = Pipeline([("scale", StandardScaler()), ("model", Ridge())])

GridSearchCV takes four things:

  • the estimator to tune;
  • a dictionary of the hyperparameter values to try; a parameter inside a pipeline is named by its step name, two underscores, then the parameter: "model__alpha";
  • a scoring to maximize, so the RMSE is passed as its negative, "neg_root_mean_squared_error";
  • cv, the number of folds.

Once fitted, it refits the best value on the whole training set.

alphas = np.logspace(-1, 5, 25)
ridge_cv = GridSearchCV(ridge_pipe, {"model__alpha": alphas},
                        scoring="neg_root_mean_squared_error", cv=5)
ridge_cv.fit(I_train, y_train)
print(ridge_cv.best_params_)
ridge_cv.best_estimator_

Every fold score is stored in cv_results_, a dictionary that reads best as a data frame:

results = pd.DataFrame(ridge_cv.cv_results_)
results[["param_model__alpha", "mean_test_score", "std_test_score"]].head()
cv_rmse = -results["mean_test_score"]
plt.semilogx(alphas, cv_rmse, "o-", color=BLUE, ms=4)
plt.axvline(ridge_cv.best_params_["model__alpha"], color=RED, lw=1)
plt.xlabel("λ"); plt.ylabel("5-fold CV RMSE"); plt.title("ridge on the interactions")
plt.show()
print(f"λ* = {ridge_cv.best_params_['model__alpha']:.4g}   CV RMSE {cv_rmse.min():.3f}")

The lasso replaces \(\|\beta\|^2\) by \(\|\beta\|_1 = \sum_j|\beta_j|\), whose corner sets coefficients exactly to zero: it shrinks and selects at once.

✍️ Your turn. Do the same for the lasso: tune Lasso() on the interaction columns over alphas_l, then compare its CV RMSE with the ridge one.

alphas_l = np.logspace(-2, 2, 20)
# TODO: lasso_pipe = ...
# TODO: lasso_cv = GridSearchCV(...)
# TODO: lasso_cv.fit(...)

kept = int((lasso_cv.best_estimator_["model"].coef_ != 0).sum())
print(f"ridge CV RMSE {-ridge_cv.best_score_:.3f}")
print(f"lasso CV RMSE {-lasso_cv.best_score_:.3f}   "
      f"α* = {lasso_cv.best_params_['model__alpha']:.4g}   columns kept {kept} of {I_train.shape[1]}")

Now, and only now, the test set:

scores["interactions, ridge"] = rmse(y_test, ridge_cv.predict(I_test))
scores["interactions, lasso"] = rmse(y_test, lasso_cv.predict(I_test))
for k in ["numbers only, OLS", "interactions, OLS", "interactions, ridge", "interactions, lasso"]:
    print(f"{k:22s} test RMSE {scores[k]:.3f}")

🧠 Explain. Did regularization help on the interactions? Compare with OLS on the 32 numeric columns.

✍️ Your answer:

…

Part 2 — The columns that are not numbers

Now we turn to the categorical features: the neighborhood, the zoning, the kitchen grade, the roof, the foundation. Lab 3 turned one such column into dummies by hand. Here there are 28 of them.

2.1 One-hot encoding

OneHotEncoder turns a categorical column into one 0/1 column per level. Like StandardScaler, it is a transformer: .fit learns the levels of each column, .transform encodes them.

✍️ Your turn. Fit a OneHotEncoder() on the categorical columns of the training set and transform them. Compare the number of encoded columns with the number of levels: was a level dropped?

# TODO: ohe = ...
# TODO: C_train = ...

levels = sum(df_train[c].nunique() for c in cat_cols)
print(f"{levels} levels in the training set · {C_train.shape[1]} encoded columns")

✍️ Your answer:

…

🧠 Explain. Lab 3 dropped one level per categorical column. Why is it better not to drop one with ridge? Hint: look at the coefficients and how they are penalized.

✍️ Your answer:

…

✍️ Your turn. Transform the categorical columns of the test set with the same encoder.

# TODO: C_test = ...

🧠 Explain. What does the error say, and why does it happen on the test set and not on the training set?

✍️ Your answer:

…

With handle_unknown="ignore", an unknown level is encoded as all zeros: none of that column’s dummies is switched on. Here is one of the three test houses with a wood-shake roof, a material no training house has:

ohe = OneHotEncoder(handle_unknown="ignore").fit(df_train[cat_cols])
C_test = ohe.transform(df_test[cat_cols])

row = df_test.loc[df_test["RoofMatl"] == "WdShake", cat_cols].head(1)
encoded = pd.Series(ohe.transform(row).toarray()[0], index=ohe.get_feature_names_out())
encoded[encoded.index.str.startswith("RoofMatl_")]

2.2 Numbers and categories together

We want the scaler on the numeric columns and the encoder on the categorical ones. A ColumnTransformer applies each transformer to its own columns and glues the results side by side. It is a transformer too, with the same .fit / .transform API.

prep = ColumnTransformer([("num", StandardScaler(), num_cols),
                          ("cat", OneHotEncoder(handle_unknown="ignore"), cat_cols)],
                         verbose_feature_names_out=False)
Z_train = prep.fit_transform(df_train)
print(f"encoded design: {Z_train.shape[1]} columns for {Z_train.shape[0]} training houses")

🔮 Predict. Will a ridge on this design beat the 32 plain numeric columns? Why?

✍️ Your answer:

…

Like the scaler in 1.3, the whole transformer goes inside the pipeline: every cross-validation refit then redoes the scaling and the encoding on its own rows.

✍️ Your turn. Put prep and a Ridge() in one pipeline, then choose λ over np.logspace(-2, 4, 25) by 5-fold cross-validation on the training set.

# TODO: full_pipe = ...
# TODO: grid = GridSearchCV(...)
# TODO: grid.fit(...)

scores["numbers + categories, ridge"] = rmse(y_test, grid.predict(df_test))
print(f"λ* = {grid.best_params_['model__alpha']:.4g}   "
      f"CV RMSE {-grid.best_score_:.3f}   "
      f"test RMSE {scores['numbers + categories, ridge']:.3f}")
median_price = np.exp(y_train).median()          # the median training price, for a dollar reading
table = pd.DataFrame({"test RMSE": pd.Series(scores)})
table["typical error"] = (100 * (np.exp(table["test RMSE"]) - 1)).round(1).astype(str) + " %"
table[f"on a {median_price:,.0f} $ house"] = (median_price * (np.exp(table["test RMSE"]) - 1)).round(-2).map("±{:,.0f} $".format)
table.sort_values("test RMSE")

🧠 Explain. Write the agency’s three-sentence report: which model you ship, what it will cost them in dollars on a typical house, and one sentence on why the biggest gain of the lab came from columns you had ignored rather than from any penalty.

✍️ Your answer:

…

2.3 Which columns does the lasso keep?

✍️ Your turn. Build the same pipeline with Lasso(max_iter=20_000) in place of Ridge() and tune α over alphas_full. Then read the coefficients of the best model: the fitted pipeline is in best_estimator_, and a linear model stores its coefficients in its coef_ attribute.

alphas_full = np.logspace(-4, 0, 20)
# TODO: lasso_full = GridSearchCV(...)
# TODO: lasso_full.fit(...)

best = lasso_full.best_estimator_
# TODO: coefs to be defined as the Series of the best model's coefficients, indexed by column name

kept = coefs[coefs != 0]
print(f"α* = {lasso_full.best_params_['model__alpha']:.4g}   "
      f"test RMSE {rmse(y_test, lasso_full.predict(df_test)):.3f}   "
      f"{len(kept)} of {len(coefs)} columns kept")
kept.abs().sort_values(ascending=False).head(10)

🧠 Explain. The model is log(price) = β0 + β1 x1 + β2 x2 + … Suppose β1 = 0.1 and x1 increases by 1, all else equal. How is the price affected?

✍️ Your answer:

…

🧠 Explain. Which columns does the lasso keep, and which matter most?

✍️ Your answer:

…

The plot refits the lasso at every α of the grid and draws the coefficients of the eight columns with the largest ones at α*:

Z_train = prep.fit_transform(df_train)
path = pd.DataFrame([Lasso(alpha=a, max_iter=20_000).fit(Z_train, y_train).coef_ for a in alphas_full],
                    index=alphas_full, columns=prep.get_feature_names_out())
top = kept.abs().sort_values(ascending=False).index[:8]

ax = path[top].plot(logx=True, figsize=(7, 4))
ax.axvline(lasso_full.best_params_["model__alpha"], color=GRAY, ls="--", lw=1)
ax.axhline(0, color="black", lw=0.5)
ax.set_xlabel("α"); ax.set_ylabel("coefficient"); ax.legend(fontsize=7)
plt.show()

🧠 Explain. Read the plot: what happens to the coefficients as α grows, and which columns are the last to leave the model?

✍️ Your answer:

…