Lab 3 — Reading regression output like a statistician

Author

Hicham Janati — FSML

You’re a statistician analyzing house sales in Ames, Iowa (2006–2010). ames_housing.csv holds one row per sale: the price and eight covariates. It is a cleaned subset of the data published by De Cock (2011) and used in the Kaggle competition House Prices. Prices are in dollars, areas in m².

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. AI is fine for syntax; the written interpretations must be yours.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats

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

def resid_vs_fitted(fit, ax, title=""):
    """Residuals against fitted values, with the spread in three fitted-value terciles."""
    r, f = fit.resid, fit.fittedvalues
    ax.scatter(f, r, s=8, color=GRAY, alpha=0.6)
    ax.axhline(0, color=RED, lw=1)
    terc = pd.qcut(f, 3, labels=False)
    sds = [r[terc == k].std() for k in range(3)]
    ax.set_title(f"{title}\nresidual sd by fitted tercile: " +
                 ", ".join(f"{s:.3g}" for s in sds), fontsize=9)
    ax.set_xlabel("fitted"); ax.set_ylabel("residual")

def qq(fit, ax, title=""):
    """QQ plot of the standardized residuals against Gaussian quantiles."""
    z = fit.get_influence().resid_studentized_internal
    z = z[np.isfinite(z)]                   # a row with leverage 1 has no studentized residual
    sm.qqplot(z, line="45", ax=ax, markersize=2)
    ax.set_title(f"{title}\nexcess kurtosis {stats.kurtosis(z):.1f}, skew {stats.skew(z):.2f}",
                 fontsize=9)

def resid_vs_x(fit, x, ax, title="", bins=10):
    """Residuals against one covariate, with the mean residual per decile of x."""
    ax.scatter(x, fit.resid, s=6, color=GRAY, alpha=0.4)
    ax.axhline(0, color=RED, lw=1)
    b = pd.qcut(x, bins, labels=False, duplicates="drop")
    ax.plot(x.groupby(b).median(), fit.resid.groupby(b).mean(), "o-", color=BLUE, ms=4)
    ax.set_title(title, fontsize=10); ax.set_xlabel(x.name); ax.set_ylabel("residual")

Part 1 — The data, and one column that is not a number

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

column	meaning
`SalePrice`	sale price, dollars
`GrLivArea`	above-ground living area, m²
`OverallQual`	overall material and finish quality, rated 1 (very poor) to 10 (very excellent)
`YearBuilt`	construction year
`LotArea`	lot size, m²
`GarageCars`	garage capacity, cars
`BedroomAbvGr`	bedrooms above ground
`PoolArea`	pool area, m² (0 when there is no pool)
`BldgType`	building type, five categories (see below)

The five building types: 1Fam is a detached single-family house. 2fmCon is a two-family conversion: a house built for one family and later split into two units (two kitchens, two entrances, one building). Duplex was built from the start as two side-by-side units under one roof. A townhouse is a row house sharing walls with its neighbors: Twnhs is an inside unit (a neighbor on both sides, windows only front and back), TwnhsE an end unit (one shared wall, one free side, usually a bigger lot).

df.describe().round(1)

1.1 Dummy variables

BldgType is text. A regression needs numbers, and the wrong fix is to code the five types 1, 2, 3, 4, 5: that would say a duplex is “one more” than a two-family conversion, and a townhouse “twice” a duplex. The right fix is one column per category, holding 1 if the house is of that type and 0 otherwise: dummy variables, also called one-hot encoding.

print(df["BldgType"].value_counts())
X_dum = pd.get_dummies(df["BldgType"], dtype=float)
X_dum.iloc[[0, 5, 9]]                               # three houses, five columns

🔮 Predict. The five dummy columns of one house always sum to 1. The model also has an intercept column, all ones. What does the lecture say about a design matrix (matrix X) whose columns are linearly dependent? So how many of the five columns can we keep?

✍️ Your answer:

…

D = pd.get_dummies(df, columns=["BldgType"], drop_first=True, dtype=float)
features = [c for c in D.columns if c != "SalePrice"]
print(features)

drop_first=True removes BldgType_1Fam: a house is now a single-family house when all the other kept dummies are 0. This is the reference level.

🧠 Explain. In a regression of price on these features, what exactly will the coefficient of BldgType_Duplex measure? Say what “other covariates” are held fixed.

✍️ Your answer:

…

1.2 Look before you fit: the pair plot

🔮 Predict. Of the four covariates in the plot below, which correlates most with price? And name one pair of covariates you expect to be correlated with each other, with a one-line reason.

✍️ Your answer:

…

numeric = ["SalePrice", "OverallQual", "GrLivArea", "LotArea", "YearBuilt"]
pd.plotting.scatter_matrix(df[numeric], figsize=(10, 10), s=4, alpha=0.4,
                           color=GRAY, hist_kwds={"bins": 30, "color": GRAY})
plt.suptitle("pair plot: scatter of each pair, histogram on the diagonal", y=0.92)
plt.show()
df[numeric].corr().round(2)

🧠 Explain. Look at the pair plot. Which two covariates are most correlated? What does it mean for the coefficients of these two covariates?

✍️ Your answer:

…

🧠 Explain. The diagonal: which two histograms are far from symmetric? Keep them in mind for Part 3.

✍️ Your answer:

…

Part 2 — First fit: the output table, line by line

X = sm.add_constant(D[features])
y = D["SalePrice"]
fit_raw = sm.OLS(y, X).fit()
print(fit_raw.summary().tables[1])
print(f"R² = {fit_raw.rsquared:.3f}")

🧠 Explain. PoolArea has p = 0.83. Does that mean a swimming pool adds nothing to the price of a house?

✍️ Your answer:

…

🧠 Explain. Your manager wants to remove every non-significant variable and refit. Is this justified in general?

✍️ Your answer:

…

🧠 Explain. Exactly 5 of the 1 456 houses have a pool. Explain the size of that standard error in one sentence. Then say what you would need in order to answer “does a pool add value”, and whether any amount of cleverness on this file could get you there.

✍️ Your answer:

…

🧠 Explain. LotArea is significant, p < 0.001, but its coefficient is 8.9 $ per m² against 789 $ per m² of living area, 89 times smaller. Does that make lot area practically unimportant to a house price? Answer yes or no, in one sentence, before you run the next cell.

✍️ Your answer:

…

A coefficient is a price per one more m², and no house differs from another by a single m². A room is about 20 m² of living area; lots differ by hundreds of m². The statistical version of “how much does this column actually move around” is to shift each column by one standard deviation, its typical variation: coef × column std. We do this for all numerical variables:

numerical_features = ["GrLivArea", "LotArea", "PoolArea"]
sd = D[numerical_features].std()
impact_numerical = (fit_raw.params[numerical_features] * sd)
print(impact_numerical.round(2))

Part 3 — Diagnostics, one plot at a time

3.1 The spread of the residuals

fig, ax = plt.subplots()
resid_vs_fitted(fit_raw, ax, "raw price model")
plt.show()

🧠 Explain. Describe the plot. Which columns of the table you just read stop being trustworthy?

✍️ Your answer:

…

3.2 A usual fix

The funnel says: the error the model makes scales with the price. An agent who misjudges a house by 10 % is off by 15 000 $ on a 150 000 $ house and by 50 000 $ on a 500 000 $ house. The noise is multiplicative:

\[\text{price} = e^{\beta_0 + \beta^\top x}\cdot e^{\varepsilon} \quad\Longleftrightarrow\quad \log \text{price} = \beta_0 + \beta^\top x + \varepsilon .\]

On the log scale the same 10 % error is the same 0.1 everywhere: a constant spread is exactly (A3).

✍️ Your turn. Refit the model of Part 2 on log(SalePrice): same X, logy in place of y. Print the coefficient table, then draw the residuals against the fitted values with resid_vs_fitted, as in 3.1.

logy = np.log(y)
# TODO: fit_log = ... (Part 2's fit, with logy in place of y)
# TODO: print the coefficient table of fit_log
fig, ax = plt.subplots()
# TODO: residuals vs fitted of fit_log (resid_vs_fitted, as in 3.1)
plt.show()

🧠 Explain. Verdict on (A3) now: did the log transformation fix the funnel ?

✍️ Your answer:

…

3.3 The tails

fig, axes = plt.subplots(1, 2, figsize=(10, 3.8))
qq(fit_raw, axes[0], "raw price model")
qq(fit_log, axes[1], "log price model")
plt.tight_layout(); plt.show()

🧠 Explain. Verdict on (A4) for each model, with the plot as evidence. Which side of the log model’s QQ plot breaks, and how many points are responsible?

✍️ Your answer:

…

3.4 Without (A4): the asymptotic table

Section 3.3 said (A4) cannot reasonably be assumed for either model. Every p-value and confidence interval read so far rests on the law $T_j = \hat\beta_j / \widehat{\mathrm{SE}}(\hat\beta_j) \sim t(n-p-1)$, which is exact under (A1)–(A4) and lost without (A4). The lecture’s answer: for large $n$ the central limit theorem makes $T_j$ approximately Gaussian, so the asymptotic p-values stay valid. statsmodels switches to that law with use_t=False: the column t becomes z, and the Gaussian quantile replaces the $t$ quantile in the p-values and the confidence intervals.

🔮 Predict. The log model has about 1 450 rows and 12 coefficients. Before running: will the SEs change? How much will the p-values and the intervals change?

✍️ Your answer:

…

fit_log_z = sm.OLS(logy, X).fit(use_t=False)
print(fit_log_z.summary().tables[1])
print(f"97.5 % quantile: t({int(fit_log.df_resid)}) = {stats.t.ppf(0.975, fit_log.df_resid):.4f}"
      f"   Gaussian = {stats.norm.ppf(0.975):.4f}")

🧠 Explain. Compare with the table of 3.2, column by column: what changed, what did not, and why is the difference this size?

✍️ Your answer:

…

3.5 Beyond homoscedasticity: robust standard errors

Here the log repaired the funnel. But what if it does not? We need to make the variance estimator robust to heteroscedasticity. There are several Heteroscedasticity-Consistent (HC) estimators. We will use HC3, a simple and popular one.

The idea is to stop estimating a single $\sigma^2$ and to estimate $\operatorname{Var}(\hat\beta)$ directly, with an estimator that converges to the true $\operatorname{Var}(\hat\beta)$ whatever the spread does.

🤖 With AI. Find the API rather than being told it. Ask an assistant: “In Python statsmodels, how do I fit an OLS model with HC3 heteroscedasticity-robust standard errors?” Run what it gives you below, then compare the table with the one from Part 2: which columns changed, and which did not?

✍️ Your answer:

…

# TODO: fit_raw_hc = ... (the raw model, HC3 robust covariance)
# TODO: print its coefficient table and compare with Part 2

3.6 The bootstrap

We end this lab with a more general method for uncertainty estimation that works with any model: the bootstrap.

A standard error answers: if I could collect another 1 456 houses from the same market, how much would $\hat\beta$ move? We cannot collect them. So we draw a new sample of 1 456 rows from the rows we have, with replacement, refit, and repeat. The spread of those refits estimates the spread we could not observe.

🔮 Predict. Draw 1 456 rows with replacement from 1 456 rows. Roughly what fraction of the original houses ends up in the draw at least once?

✍️ Your answer:

…

✍️ Your turn. Draw one resample of row positions (integers) and count how many distinct houses it contains. Hint: you can use rng.integers for generating the row positions and np.unique to count the distinct houses.

rng = np.random.default_rng(0)
# TODO: the array idx to be defined as n row positions, drawn with replacement
# TODO: how many distinct houses is that?

One resample is one imaginary alternative dataset. Refit on it and you get one alternative $\hat\beta$. Do it a thousand times and the standard deviation of those thousand estimates is the standard error, measured rather than derived.

✍️ Your turn. Write the loop: at each of the B rounds, draw a fresh resample, refit the raw model on it, and keep the coefficients in row b of boot.

rng = np.random.default_rng(0)
B = 1000

boot = np.zeros((B, X.shape[1]))                                
# TODO: loop over the B rounds
    # TODO: idx, a fresh resample
    # TODO: boot[b], the coefficients refit on the resample

✍️ Your turn. The bootstrap standard error of a coefficient is the standard deviation of its B estimates. Below you have the code to compute the SE and store it in a pandas Series, all you have to do is decide which axis ?

# TODO: the int axis to be defined as the axis along which to compute the SE
boot_se = pd.Series(boot.std(axis=axis, ddof=1), index=X.columns) 

compare = pd.DataFrame({"classical": fit_raw.bse, "HC3": fit_raw_hc.bse,
                        "bootstrap": boot_se})
print(compare.drop("const").round(3))