BSAN 750 Lab 3: Linear Regression with Boston Housing

Foundations of Machine Learning and AI

How are housing values associated with neighborhood characteristics?

This lab adapts the earlier R regression lab using Boston Housing, the dataset available in R as MASS::Boston. We progress from exploratory analysis to model fitting, inference, categorical predictors, interactions, diagnostics, prediction, and model assessment.

Keep Boston.csv in the same folder as this notebook. In Colab, upload both files. The data are read locally; no R installation or dataset download is required when running the lab.

Learning objectives

1. Fit and interpret simple and multiple linear regressions.
2. Obtain coefficient tests, confidence intervals, fitted values, and residuals.
3. Represent categorical predictors and interpret an interaction.
4. Distinguish confidence intervals for a conditional mean from prediction intervals for a new observation.
5. Diagnose model assumptions and explore an improved functional form.
6. Evaluate predictions with training/test samples and cross-validation.

How to work: RUN cells are worked examples. PRACTICE cells leave room for your code. Answer the questions before running each section, then interpret the output. The final bootstrap extension is optional.

1. Set up Python

Before coding: Which package stores tables? Which provides coefficient tests and confidence intervals for a regression?

We use pandas for tables, numpy for numerical operations and sampling, matplotlib for figures, and statsmodels for regression. If a package is missing, run this once in a separate notebook cell:

%pip install numpy pandas matplotlib statsmodels

We implement the split and cross-validation explicitly so you can see which rows are used for fitting and evaluation.

# RUN: Import the tools used in this lab.
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
from IPython.display import display
from sklearn.model_selection import train_test_split

pd.set_option("display.precision", 3)
plt.rcParams.update({"figure.figsize": (8, 5), "axes.grid": True, "grid.alpha": 0.2})
print("Setup complete.")

2. Load and understand Boston Housing

Before coding: What does one row represent? Is the response an individual house price?

The dataset contains 506 area-level observations and 14 variables. The response medv is the median value of owner-occupied homes, measured in thousands of dollars. Consequently, a prediction interval for a new row concerns another area’s median value, not the sale price of an individual house.

Variable	Meaning
medv	Median owner-occupied home value (thousands of dollars), the response
rm	Average rooms per dwelling
crim	Per-capita crime rate
lstat	Percentage described by the original source as lower-status population
chas	1 if the tract borders the Charles River; 0 otherwise
rad	Highway-accessibility index
zn	Share of residential land zoned for large lots
indus	Share of land in non-retail business use
nox	Nitrogen oxides concentration
age	Share of older owner-occupied homes
dis	Weighted distance to employment centers
tax	Property-tax rate
ptratio	Pupil–teacher ratio
black	Historical race-derived measure, \(1000(B_k-0.63)^2\), where \(B_k\) is the area’s Black population share

This historical dataset includes socially loaded variables and area-level measurements.

Sources: MASS data documentation and Rdatasets CSV copy. The CSV also includes rownames, an export identifier that we use only as a row index.

# RUN: Upload the CSV only if it is missing in Colab.
data_path = Path("Boston.csv")
if not data_path.exists():
    try:
        from google.colab import files
    except ImportError:
        raise FileNotFoundError("Put Boston.csv in the notebook's working folder.")
    files.upload()

# rownames is an identifier, not a predictor.
boston = pd.read_csv(data_path, index_col="rownames")
boston.index.name = "Observation"
print("Rows and variables:", boston.shape)
display(boston.head())
display(boston.dtypes.rename("Data type"))
display(boston.isna().sum().rename("Missing values"))

Explore distributions and relationships

Before coding: Which variables might help predict housing values? Which numeric columns represent categories or codes?

Restrict a scatterplot matrix to a few quantitative variables so it remains readable. Category codes such as chas should be summarized separately.

# RUN: Summarize the data and visualize selected numeric relationships.
display(boston.describe().T)
display(boston["chas"].value_counts().sort_index().rename("River indicator counts"))
display(boston["rad"].value_counts().sort_index().rename("Accessibility-index counts"))

pd.plotting.scatter_matrix(
    boston[["medv", "rm", "crim", "lstat"]],
    figsize=(10, 10), diagonal="hist", alpha=0.35, s=12
)
plt.suptitle("Selected Boston Housing variables", y=1.01)
plt.show()

Your turn

1. How many observations and variables are in the data? Are values missing?
2. Which variables are quantitative, and which require care because they are codes?
3. Describe the relationships between medv and rm, and between medv and lstat.
4. How does the distribution of crim look like?

3. Simple linear regression: housing value and rooms

Before coding: Predict the sign of the rooms coefficient and write the model equation.

\[medv=\beta_0+\beta_1 rm+\epsilon.\]

We begin with a plain scatterplot. Then we fit OLS and add a line in a separate figure.

# RUN: Draw the plain scatterplot.
fig, ax = plt.subplots()
ax.scatter(boston["rm"], boston["medv"], color="#4B5563", s=20, alpha=0.5)
ax.set(xlabel="Average rooms per dwelling (rm)", ylabel="Median home value ($1,000s)",
       title="Boston housing value and average rooms")
fig.tight_layout()
plt.show()

Fit the model and add its line

Before coding: What do ols() and fit() each do? Why use a sorted grid to draw the fitted line?

The formula specifies the response and predictor. ols() sets up the model, and fit() estimates coefficients. An intercept is included automatically.

# RUN: Estimate the model, inspect its summary, and plot fitted values.
simple_model = smf.ols("medv ~ rm", data=boston).fit()
print(simple_model.summary())

rooms_grid = np.linspace(boston["rm"].min(), boston["rm"].max(), 120)
line_data = pd.DataFrame({"rm": rooms_grid})
fig, ax = plt.subplots()
ax.scatter(boston["rm"], boston["medv"], color="#4B5563", s=20, alpha=0.5)
ax.plot(rooms_grid, simple_model.predict(line_data), color="#2563EB",
        linewidth=2, label="Fitted regression")
ax.set(xlabel="Average rooms per dwelling (rm)", ylabel="Median home value ($1,000s)",
       title="Simple regression: medv ~ rm")
ax.legend()
fig.tight_layout()
plt.show()

Extract coefficients, tests, and confidence intervals

Before coding: What null hypothesis does the rooms p-value test? Does the intercept have a meaningful interpretation here?

The slope is the expected change in area median value, in thousands of dollars, associated with a one-unit increase in average rooms under this model. It is not an estimate of what adding a room to one particular house would do to its selling price.

# RUN: Build a compact coefficient table.
coefficient_ci = simple_model.conf_int(alpha=0.05)
display(pd.DataFrame({
    "Estimate": simple_model.params,
    "Std. error": simple_model.bse,
    "t statistic": simple_model.tvalues,
    "p-value": simple_model.pvalues,
    "95% lower": coefficient_ci[0],
    "95% upper": coefficient_ci[1],
}))

Fitted values, residuals, and MSE

Before coding: Verify \(e_i=y_i-\hat y_i\). Why is residual MSE not the average squared training error?

Let \(k\) be the number of estimated coefficients, including the intercept.

• Residual MSE estimates \(\sigma^2\) using \(SSE/(n-k)\); Python returns it as mse_resid.
• Average squared training error uses \(SSE/n\).
• Residual standard error is \(\sqrt{SSE/(n-k)}\), analogous to sigma in R.

Neither is a direct measure of prediction error on new observations.

# RUN: Extract fitted values and residuals and calculate fit measures.
fitted_table = boston[["rm", "medv"]].copy()
fitted_table["Fitted medv"] = simple_model.fittedvalues
fitted_table["Residual"] = simple_model.resid
display(fitted_table.head())

sse = np.sum(simple_model.resid ** 2)
n = int(simple_model.nobs)
k = len(simple_model.params)  # Includes the intercept.
display(pd.Series({
    "n": n,
    "Number of coefficients (k)": k,
    "Residual degrees of freedom": simple_model.df_resid,
    "Residual MSE: SSE / (n - k)": sse / (n - k),
    "Training MSE: SSE / n": sse / n,
    "Residual standard error": np.sqrt(simple_model.mse_resid),
    "R-squared": simple_model.rsquared,
    "Adjusted R-squared": simple_model.rsquared_adj,
}))
print("Residual identity holds:",
      np.allclose(simple_model.resid, boston["medv"] - simple_model.fittedvalues))

Your turn

1. Write the fitted equation and interpret the rooms slope in dollars.
2. Is that slope significantly different from zero at the 5% level?
3. Compute the predicted median value for an area with rm = 6 by hand, then verify with predict().
4. Explain why predictions at rm = 15 would be extrapolation.

# PRACTICE: Calculate and verify the fitted value at rm = 6.
# Hints: simple_model.params["Intercept"], simple_model.params["rm"]
# Check with simple_model.predict(pd.DataFrame({"rm": [6]})).

# YOUR CODE HERE

4. Multiple regression

Before coding: Could the rooms coefficient change after adjusting for other neighborhood characteristics?

\[medv=\beta_0+\beta_1 crim+\beta_2 rm+\beta_3 lstat+\epsilon.\]

Each slope now describes an association while holding the other included predictors constant.

# RUN: Fit the three-predictor model.
multiple_formula = "medv ~ crim + rm + lstat"
multiple_model = smf.ols(multiple_formula, data=boston).fit()
print(multiple_model.summary())

5. Categorical predictors

A binary predictor: proximity to the Charles River

Before coding: If chas is coded 0/1, how many additional coefficients are needed? What does its coefficient compare?

# RUN: Add the river indicator to the multiple regression.
river_formula = "medv ~ crim + rm + lstat + chas"
river_model = smf.ols(river_formula, data=boston).fit()
print(river_model.summary())

Question: How to interpret the coefficient for chas?

More than two categories: highway-accessibility index

Before coding: How many distinct rad values occur? If we treat them as categories, how many indicator coefficients are required?

To illustrate categorical encoding, treat each observed rad level separately using C(rad), with level 1 as the reference. This avoids assuming that a one-unit step has the same association everywhere. It also uses more coefficients than a single numeric slope, so this is a modeling choice to evaluate.

# RUN: Set a fixed category order from the observed accessibility codes.
boston["rad"] = pd.Categorical(boston["rad"])
print("Accessibility levels:", rad_levels)

# code below is to illustrate how the dummy variables are converted from the original categorical variable
indicators = pd.get_dummies(boston["rad"], prefix="rad", drop_first=True, dtype=int)
display(pd.concat([boston[["rad"]], indicators], axis=1).head())

# fit model
categorical_formula = "medv ~ crim + rm + lstat + chas + C(rad)"
categorical_model = smf.ols(categorical_formula, data=boston).fit()
print(categorical_model.summary())

Question: How to interpret the coefficients for rad? Is this variable statistically significant?

Your turn

1. Interpret the chas coefficient in the categorical model.
2. Interpret one C(rad) coefficient relative to level 1.
3. Count all estimated coefficients, including the intercept, and verify \(df_{resid}=n-k\).
4. How would changing the reference category affect coefficient labels and fitted values?

# PRACTICE: Verify the residual degrees of freedom.
# Hints: categorical_model.nobs, len(categorical_model.params)

# YOUR CODE HERE

# PRACTICE: Change the reference category.
categorical_formula1 = "medv ~ crim + rm + lstat + chas + C(rad, Treatment(reference=24))"
categorical_model1 = smf.ols(categorical_formula2, data=boston).fit()
print(categorical_model.summary())

6. An interaction: rooms × river proximity

Before coding: Could the rooms association differ between areas bordering the river and other areas?

rm * chas includes rm, chas, and their product. We keep crime, lstat, and the accessibility indicators as controls.

\[medv=\beta_0+\beta_1 rm+\beta_2 chas+ \beta_3(rm\times chas)+\text{other predictor terms}+\epsilon.\]

The non-river rooms slope is \(\beta_1\); the river slope is \(\beta_1+\beta_3\). The river difference at \(r\) rooms is \(\beta_2+\beta_3(r-6)\). Thus the chas main effect is the group difference specifically at six rooms. Do not assume that an interaction must be statistically significant.

# RUN: Fit the interaction and compare it with the additive model.
interaction_formula = "medv ~ crim + lstat + C(rad) + rm * chas"
interaction_model = smf.ols(interaction_formula, data=boston).fit()
print(interaction_model.summary())

b = interaction_model.params
display(pd.Series({
    "Non-river rooms slope": b["rm"],
    "River rooms slope": b["rm"] + b["rm:chas"],
    "Difference in slopes": b["rm:chas"],
    "River difference": b["chas"],
}))

Questions: 1. Interpret the interaction coefficient and its p-value. 2. Why keep both main effects when including their interaction?

7. Prediction and uncertainty

Before coding: What input columns are needed to predict from the three-predictor model? Which interval should be wider?

get_prediction().summary_frame() returns both a confidence interval for the conditional mean (mean_ci_lower/upper) and a prediction interval for one new area-level response (obs_ci_lower/upper). These pointwise 95% intervals rely on the fitted model and its assumptions.

# RUN: Predict for two illustrative areas with different average room counts.
new_areas = pd.DataFrame({
    "crim": [0.1, 0.1],
    "rm": [6.0, 7.0],
    "lstat": [10.0, 10.0],
})
prediction_summary = multiple_model.get_prediction(new_areas).summary_frame(alpha=0.05)
display(pd.concat([new_areas, prediction_summary], axis=1))

Confidence Interval vs. Prediction Interval

For a given set of predictor values, regression can provide two different types of intervals:

• Confidence interval for the mean response: estimates the range for the average value of \(Y\) among all observations with these predictor values.
• Prediction interval: estimates the range for one new individual observation of \(Y\) with these predictor values.

Conceptually,

\[ \text{Confidence interval: } \hat{Y} \pm \text{uncertainty in estimating the mean} \]

\[ \text{Prediction interval: } \hat{Y} \pm \text{uncertainty in estimating the mean + individual variation} \]

Therefore, the prediction interval is always wider than the confidence interval because predicting an individual observation involves additional uncertainty.

Visualize confidence and prediction intervals

Before coding: Why does the prediction interval include more uncertainty than the confidence interval?

Return to the simple rooms model for a two-dimensional figure. The prediction interval still concerns an area’s medv, not one house price.

# RUN: Draw both interval bands for the simple regression.
bands = simple_model.get_prediction(line_data).summary_frame(alpha=0.05)
fig, ax = plt.subplots()
ax.fill_between(rooms_grid, bands["obs_ci_lower"].to_numpy(),
                bands["obs_ci_upper"].to_numpy(), color="#DBEAFE",
                label="95% prediction interval")
ax.fill_between(rooms_grid, bands["mean_ci_lower"].to_numpy(),
                bands["mean_ci_upper"].to_numpy(), color="#93C5FD",
                label="95% confidence interval for mean")
ax.scatter(boston["rm"], boston["medv"], color="#4B5563", s=15, alpha=0.35)
ax.plot(rooms_grid, bands["mean"], color="#2563EB", linewidth=2)
ax.set(xlabel="Average rooms per dwelling", ylabel="Median home value ($1,000s)",
       title="Uncertainty in the simple rooms model")
ax.legend()
fig.tight_layout()
plt.show()

8. Model diagnostics

Before coding: What patterns suggest nonlinearity, unequal error variance, or influential cases?

The following helper creates Python analogues of the four usual R plot(lm(...)) diagnostics. You can reuse it without memorizing every plotting line.

Plot	What to inspect
Residuals vs. fitted	Curvature or other systematic patterns
Normal Q-Q	Deviations from the normal reference line
Scale-location	Changing residual spread
Residuals vs. leverage	Unusual predictor combinations together with large residuals

Normality concerns model errors, not the distributions of predictors or raw home values. Independence depends on data collection; a random split does not remove possible spatial dependence among nearby areas.

# RUN: Define a reusable four-panel diagnostic figure.
def regression_diagnostics(model, title, fitted_label="Fitted medv ($1,000s)"):
    influence = model.get_influence()
    fitted = np.asarray(model.fittedvalues)
    residuals = np.asarray(model.resid)
    standardized = influence.resid_studentized_internal
    leverage = influence.hat_matrix_diag
    cooks_d = influence.cooks_distance[0]
    scale = np.sqrt(np.abs(standardized))
    fig, axes = plt.subplots(2, 2, figsize=(12, 8))

    # LOWESS draws a smooth trend to help reveal patterns.
    for ax, x, y, xlabel, ylabel in [
        (axes[0, 0], fitted, residuals, fitted_label, "Residuals"),
        (axes[1, 0], fitted, scale, fitted_label, "Sqrt(|standardized residual|)"),
        (axes[1, 1], leverage, standardized, "Leverage", "Standardized residuals"),
    ]:
        ax.scatter(x, y, color="#6B7280", s=30)
        smooth = sm.nonparametric.lowess(y, x, frac=0.7)
        ax.plot(smooth[:, 0], smooth[:, 1], color="#C8102E", linewidth=1.5)
        ax.set(xlabel=xlabel, ylabel=ylabel)

    axes[0, 0].axhline(0, color="black", linestyle="--", linewidth=0.8)
    axes[1, 1].axhline(0, color="black", linestyle="--", linewidth=0.8)
    sm.qqplot(standardized, line="s", ax=axes[0, 1])
    axes[0, 0].set_title("Residuals vs. fitted")
    axes[0, 1].set(title="Normal Q-Q", ylabel="Standardized residuals")
    axes[1, 0].set_title("Scale-location")
    axes[1, 1].set_title("Residuals vs. leverage")

    # Cook's-distance contours combine residual size with leverage.
    k = len(model.params)  # Coefficients, including the intercept.
    h_max = min(0.95, max(0.3, leverage.max() * 1.15))
    h_grid = np.linspace(0.001, h_max, 250)
    y_limit = max(3, np.abs(standardized).max() * 1.25)
    for distance in [0.5, 1.0]:
        contour = np.sqrt(distance * k * (1 - h_grid) / h_grid)
        axes[1, 1].plot(h_grid, contour, "--", color="#94A3B8",
                        label=f"Cook's D = {distance}")
        axes[1, 1].plot(h_grid, -contour, "--", color="#94A3B8")
    # Label the original observation identifiers.
    for i in np.argsort(cooks_d)[-3:]:
        axes[1, 1].annotate(str(model.model.data.row_labels[i]), (leverage[i], standardized[i]),
                            xytext=(4, 4), textcoords="offset points", fontsize=8)
    axes[1, 1].set(xlim=(0, h_max), ylim=(-y_limit, y_limit))
    axes[1, 1].legend(fontsize=8)
    fig.suptitle(title, fontsize=14)
    fig.tight_layout()
    plt.show()

Diagnose the three-predictor model

Before coding: Given the scatterplot matrix, where do you expect the linear specification to struggle? Does a high \(R^2\) guarantee satisfactory diagnostics?

# RUN: Inspect diagnostics for the model from the old R lab.
regression_diagnostics(multiple_model, "Boston: medv ~ crim + rm + lstat")

Explore a nonlinear relationship

Before coding: Would a straight line capture the relationship between lstat and housing value? Can a model with a squared predictor still be a linear regression?

A quadratic predictor keeps the model linear in its coefficients. In a formula, I(lstat ** 2) means include the arithmetic square as a predictor. Keep the linear term as well.

We fit a quadratic extension of the categorical model. Compare its diagnostics and fit with the simpler specification; an extra term does not automatically fix every assumption.

# RUN: Add a quadratic lstat term and inspect the result.
quadratic_formula = categorical_formula + " + I(lstat ** 2)"
quadratic_model = smf.ols(quadratic_formula, data=boston).fit()
print(quadratic_model.summary())
display(sm.stats.anova_lm(categorical_model, quadratic_model))
regression_diagnostics(quadratic_model, "Boston: categorical model + quadratic lstat")

Investigate influential observations

Before coding: Does high leverage necessarily mean the response is unusual? Should a large Cook’s distance automatically trigger deletion?

Leverage measures unusual predictor combinations. Cook’s distance summarizes the effect an observation can have on the fitted coefficients. The rule \(D_i>4/n\) is a screening heuristic, not a formal test or deletion rule. Examine records and context and assess sensitivity before changing the data.

# RUN: List observations with the largest Cook's distances.
influence = quadratic_model.get_influence()
influence_table = boston[["medv", "rm", "crim", "lstat", "chas", "rad"]].copy()
influence_table["Standardized residual"] = influence.resid_studentized_internal
influence_table["Leverage"] = influence.hat_matrix_diag
influence_table["Cook's D"] = influence.cooks_distance[0]
influence_table["Flag: D > 4/n"] = influence_table["Cook's D"] > 4 / len(boston)
display(influence_table.sort_values("Cook's D", ascending=False).head(10))

Compare model fit

Before coding: Must training \(R^2\) improve when predictors are added? Does improved in-sample fit guarantee better prediction?

All models below use the same response and observations. Higher adjusted \(R^2\) and lower residual MSE, AIC, or BIC can favor a specification, but these are not test errors. Predictors should not be removed mechanically by individual p-values, especially when interactions or squared terms are present.

# RUN: Compare all full-data models in one table.
full_data_models = {
    "Rooms only": simple_model,
    "Three predictors": multiple_model,
    "Additive + river": river_model,
    "Categorical accessibility": categorical_model,
    "Rooms × river interaction": interaction_model,
    "Quadratic lstat": quadratic_model,
}
fit_comparison = pd.DataFrame([
    {
        "Model": name,
        "Coefficients": len(model.params),
        "R-squared": model.rsquared,
        "Adjusted R-squared": model.rsquared_adj,
        "Residual MSE": model.mse_resid,
        "AIC": model.aic,
        "BIC": model.bic,
    }
    for name, model in full_data_models.items()
]).set_index("Model")
display(fit_comparison)

Your turn: another possible improvement

1. Identify a diagnostic problem that remains in the quadratic model.
2. Consider modeling np.log(medv) if spread increases with fitted values.
3. Fit a log-response version of the three-predictor model and inspect diagnostics.
4. Explain why exponentiating a fitted log value does not generally estimate the conditional arithmetic mean of medv.

When evaluating a log-response model, compute prediction errors on the original medv scale. A retransformation correction, if used, must be estimated using only the relevant training observations. Do not compare raw residual MSE across models with different response scales.

# PRACTICE: Explore a log-response model.
# Hint: smf.ols("np.log(medv) ~ crim + rm + lstat", data=boston).fit()
# Then call regression_diagnostics(your_model, "Log-response model",
#                                  fitted_label="Fitted log(medv)").
# The fitted_label argument sets the correct response scale on the axes.

# YOUR CODE HERE

9. Evaluate prediction: train/test split

Before coding: Why must we refit on the training sample instead of reusing a model fitted to all observations?

Use an 80/20 split. We specify candidate formulas, compare them through cross-validation inside the training sample, and evaluate the chosen model once on the test sample.

The full-data analyses above demonstrate inference and diagnostics. Because we have already explored the complete dataset, the split below illustrates the evaluation procedure rather than providing a fully untouched model-development audit. In a new project, reserve test observations before model-guided exploration. Random splits also assess a different task from predicting a new city, later period, or geographically separate region.

# Fix candidate specifications before examining test performance.
candidate_formulas = {
    "Rooms only": "medv ~ rm",
    "Three predictors": multiple_formula,
    "Categorical accessibility": categorical_formula,
    "Rooms × river interaction": interaction_formula,
    "Quadratic lstat": quadratic_formula,
}

# random split 80/20 training testing
train, test = train_test_split(boston, test_size=0.2, random_state=1234)

# train models on training sample and compute MSPE, RMSPE and MAPE on testing sample
test_results = []
for name, formula in candidate_formulas.items():
    model = smf.ols(formula, data=train).fit()
    predictions = model.predict(test)
    actual_test = test["medv"]
    errors = actual_test - predictions
    test_results.append({
        "Model": name,
        "Training MSE": model.mse_resid,
        "Training R-squared": model.rsquared,
        "Test MSPE": np.mean(errors ** 2),
        "Test RMSPE": np.sqrt(np.mean(errors ** 2)),
        "Test MAPE": np.mean(np.abs(errors)),
        "Test R-squared": 1 - np.sum(errors ** 2) / np.sum(
            (actual_test - actual_test.mean()) ** 2
        ),
    })

display(pd.DataFrame(test_results).set_index("Model").sort_values("Test RMSPE"))

10. Ten-fold cross-validation

Before coding: How often does each training observation serve as validation data? Why use identical folds for all candidates?

At the end of the cross-validation process, we need to get CV score, which can be any metrics we want to evaluate.

Note, CV-\(R^2\) uses all validation predictions together; it is not the average of ten fold-specific \(R^2\) values.

# RUN: Create the same ten folds for every candidate model.
fold_rng = np.random.default_rng(751)
folds = np.array_split(fold_rng.permutation(len(train)), 10)
cv_rows = []

# Define a helper function for cross-validation for linear regression
def run_cross_validation(data, candidate_formulas, folds):
    cv_results_rows = []
    for name, formula in candidate_formulas.items():
        oof_predictions = np.empty(len(data))
        for validation_rows in folds:
            fitting_rows = np.setdiff1d(np.arange(len(data)), validation_rows)
            fold_train = data.iloc[fitting_rows]
            fold_valid = data.iloc[validation_rows]
            # All fitting occurs on this fold's training observations.
            fold_model = smf.ols(formula, data=fold_train).fit()
            oof_predictions[validation_rows] = fold_model.predict(fold_valid)

        actual = data["medv"].to_numpy()
        errors = actual - oof_predictions
        cv_results_rows.append({
            "Model": name,
            "CV_MSPE": np.mean(errors ** 2),
            "CV_RMSPE": np.sqrt(np.mean(errors ** 2)),
            "CV_MAPE": np.mean(np.abs(errors)),
            "CV_R-squared": 1 - np.sum(errors ** 2) / np.sum((actual - actual.mean()) ** 2),
        })
    return pd.DataFrame(cv_results_rows).set_index("Model")

cv_results = run_cross_validation(train, candidate_formulas, folds).sort_values("CV_RMSPE")
display(cv_results)
selected_name = cv_results.index[0]
print("Selected using training CV:", selected_name)

Evaluate the selected model on the test sample

Before coding: How do MSE, RMSE, and MAE differ in units? Can test \(R^2\) be negative?

• Test MSE (MSPE) is average squared prediction error, measured in squared thousands of dollars.
• Test RMSE and MAE are measured in thousands of dollars.
• MAE is not MAPE, which is a percentage metric.
• Test \(R^2=1-SSE_{test}/SST_{test}\) can be negative.

We also report a baseline that predicts the training mean for every test row. This differs from the test-mean reference used in the definition of test \(R^2\). Do not switch models after inspecting the test results.

# RUN: Refit the selected model on all training observations.

test_results = []

for name, formula in candidate_formulas.items():
    # Refit each candidate model on all training observations
    model_on_full_train = smf.ols(formula, data=train).fit()
    # Make predictions on the test set
    predictions = model_on_full_train.predict(test).to_numpy()

    # Calculate errors and metrics
    actual_test = test["medv"].to_numpy()
    errors = actual_test - predictions

    test_results.append({
        "Model": name,
        "Test MSPE": np.mean(errors ** 2),
        "Test RMSPE": np.sqrt(np.mean(errors ** 2)),
        "Test MAPE": np.mean(np.abs(errors)),
        "Test R-squared": 1 - np.sum(errors ** 2) / np.sum(
            (actual_test - actual_test.mean()) ** 2
        ),
    })

display(pd.DataFrame(test_results).set_index("Model").sort_values("Test RMSPE"))

Your turn

1. Does CV favor the same specification as the full-data fit measures?
2. Why can residual MSE differ from test MSPE?
3. Why is repeatedly changing the random seed to find a favorable test result inappropriate?

11. Optional extension: bootstrap the rooms slope

Before coding: How does resampling with replacement differ from splitting rows into training and test samples?

Resample the 506 rows 500 times and refit the simple rooms model. The 2.5th and 97.5th percentiles form an approximate 95% percentile-bootstrap interval.

This illustrates sampling variability under an independent-row resampling assumption. It does not correct functional-form problems or spatial dependence. This full-data inference exercise is separate from test-set evaluation.

# RUN (OPTIONAL): Bootstrap the simple-regression rooms coefficient.
bootstrap_rng = np.random.default_rng(752)
bootstrap_slopes = []
for repetition in range(500):
    sample_rows = bootstrap_rng.integers(0, len(boston), size=len(boston))
    bootstrap_sample = boston.iloc[sample_rows]
    bootstrap_model = smf.ols("medv ~ rm", data=bootstrap_sample).fit()
    bootstrap_slopes.append(bootstrap_model.params["rm"])

bootstrap_slopes = np.asarray(bootstrap_slopes)
lower, upper = np.quantile(bootstrap_slopes, [0.025, 0.975])

## Show bootsrap std and CI, and compare with the theoretical estimtes.
theoretical_ci = simple_model.conf_int(alpha=0.05).loc["rm"]

display(pd.Series({
    "Original rooms slope": simple_model.params["rm"],
    "Bootstrap standard error": bootstrap_slopes.std(ddof=1),
    "95% bootstrap lower": lower,
    "95% bootstrap upper": upper,
    "Theoretical standard error": simple_model.bse["rm"],
    "95% theoretical lower": theoretical_ci[0],
    "95% theoretical upper": theoretical_ci[1]
}))

## show the Bootstrap distribution
fig, ax = plt.subplots()
ax.hist(bootstrap_slopes, bins=25, color="#93C5FD", edgecolor="white")
ax.axvline(simple_model.params["rm"], color="#2563EB",
           linewidth=2, label="Original estimate")
ax.axvline(lower, color="#C8102E", linestyle="--", label="95% percentile endpoints")
ax.axvline(upper, color="#C8102E", linestyle="--")
ax.axvline(theoretical_ci[0], color="#16A34A", linestyle=":", label="95% theoretical CI")
ax.axvline(theoretical_ci[1], color="#16A34A", linestyle=":")
ax.set(xlabel="Estimated rooms slope ($1,000s per room)", ylabel="Frequency",
       title="Bootstrap distribution of the rooms coefficient")
ax.legend()
fig.tight_layout()
plt.show()

Lab 3 checklist

Restart the kernel and run the notebook from top to bottom. Then write a short interpretation covering:

1. What one observation and the response variable represent.
2. A simple-regression slope versus a slope adjusted for other predictors.
3. Binary coding, the reference category for rad, and residual degrees of freedom.
4. The rooms × river interaction and its conditional interpretation.
5. Confidence intervals versus prediction intervals.
6. One diagnostic finding and an appropriate model improvement.
7. Cross-validation, test error, and the limits of the evaluation design.

Python reference

Task	Python
Fit OLS	smf.ols(formula, data=df).fit()
Regression summary	model.summary()
Coefficients and tests	model.params, model.bse, model.pvalues
Coefficient intervals	model.conf_int()
Fitted values and residuals	model.fittedvalues, model.resid
Residual variance estimate	model.mse_resid
Categorical predictor	C(rad) in the formula
Main effects plus interaction	rm_c * chas
Squared predictor	I(lstat ** 2)
Predictions	model.predict(new_data)
Mean and observation intervals	model.get_prediction(new_data).summary_frame()
Nested-model comparison	sm.stats.anova_lm(reduced, full)
Influence measures	model.get_influence()

Use this workflow as preparation for the regression analysis in Homework 3.