Worked Example: House Price Prediction

Step 1 — Load and inspect the raw file

The dataset has 1,200 rows and 13 columns. Before touching anything, look at it.

import pandas as pd
import numpy as np

df = pd.read_csv("housing.csv")

print(df.shape)              # (1200, 13)
print(df.dtypes)
print(df.isna().sum())
print(df.describe())

(1200, 13)

floor_area_m2            float64
bedrooms                 float64
bathrooms                float64
property_type             object
age_years                float64
distance_to_centre_km    float64
epc_rating                object
has_garden                object
postcode                  object
days_on_market           float64
survey_valuation         float64
agent_id                  object
sale_price               float64

Missing values per column:
floor_area_m2             14
bedrooms                   8
bathrooms                 11
age_years                  6
distance_to_centre_km      0
epc_rating                22
has_garden                 4
... (others complete)

The file is a mix of numeric and object columns, with modest missingness in a handful of numeric and categorical columns. Nothing catastrophic. The next step is the important one and it does not involve the data at all.

Step 2 — Audit the features for leakage before exploring

From Section 3.2, three columns are removed on principle:

# These columns are recorded after the sale or are a consequence of it.
LEAKY = ["days_on_market", "survey_valuation", "agent_id"]

df = df.drop(columns=LEAKY)

TARGET = "sale_price"
y = df[TARGET]
X_raw = df.drop(columns=[TARGET])

Step 3 — Split before any preprocessing

from sklearn.model_selection import train_test_split

# 60% train, 20% validation, 20% test
X_train, X_temp, y_train, y_temp = train_test_split(
    X_raw, y, test_size=0.40, random_state=42
)
X_val, X_test, y_val, y_test = train_test_split(
    X_temp, y_temp, test_size=0.50, random_state=42
)

print(f"train {len(X_train):>5}   val {len(X_val):>5}   test {len(X_test):>5}")

train   720   val   240   test   240

Every subsequent decision — imputation values, scaling statistics, feature selection, model choice — is made using the training and validation sets. The test set is locked away and will be used once, at Step 10.

Step 4 — Explore the training set only

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

print(X_train["property_type"].value_counts())
print(X_train["epc_rating"].value_counts())
print(X_train["has_garden"].value_counts())

# Implausible values?
print("floor area <= 0 :", (X_train["floor_area_m2"] <= 0).sum())
print("age < 0         :", (X_train["age_years"] < 0).sum())
print("bedrooms > 10   :", (X_train["bedrooms"] > 10).sum())
print("n postcodes     :", X_train["postcode"].nunique())

count      720.00
mean    284,318.42
std      94,207.16
min      58,000.00
25%     218,500.00
50%     276,000.00
75%     340,750.00
max     612,000.00
skewness: 0.94

property_type   detached 241 · semi 268 · terraced 174 · flat 37
epc_rating      E 41 · D 138 · C 291 · B 217 · A 33
has_garden      True 519 · False 201

floor area <= 0 : 1
age < 0         : 2
bedrooms > 10   : 1
n postcodes     : 386

Three things to act on.

• The target is right-skewed (0.94). A log transform will help. We will fit both and compare on validation rather than assuming.
• Four implausible numeric values. One non-positive floor area, two negative ages, one property with 11 bedrooms. Each needs investigation before a decision.
• 386 postcodes in 720 training rows. One-hot encoding would create 386 features averaging 1.9 observations each. That is a guaranteed overfit. We will group into broad regions instead.

Step 5 — Handle the implausible values

# Inspect the suspicious rows before deciding
print(X_train.loc[X_train["floor_area_m2"] <= 0, ["floor_area_m2", "bedrooms", "bathrooms"]])
print(X_train.loc[X_train["age_years"] < 0, ["age_years", "property_type"]])
print(X_train.loc[X_train["bedrooms"] > 10, ["bedrooms", "floor_area_m2", "property_type"]])

      floor_area_m2  bedrooms  bathrooms
1183           -1.0      3.0        2.0     ← impossible, treat as missing

      age_years  property_type
402        -5.0      terraced             ← impossible, treat as missing
917        -2.0      semi

      bedrooms  floor_area_m2  property_type
664       11.0          312.0      detached   ← large but plausible for a large house

The decisions:

• A floor area of −1 is a data-entry error. Treat it as missing and let the imputation handle it.
• Negative ages are impossible. Treat as missing.
• 11 bedrooms in a 312 m² detached house is unusual but entirely possible. Keep it. Deleting real extreme values because they are inconvenient is one of the most common mistakes in applied regression.

for frame in (X_train, X_val, X_test):
    frame.loc[frame["floor_area_m2"] <= 0, "floor_area_m2"] = np.nan
    frame.loc[frame["age_years"] < 0, "age_years"] = np.nan

Step 6 — Group the postcodes

The postcode column has too many levels to one-hot encode. We will keep the outward code — the first part, such as "SW1" — which captures broad geography while cutting the cardinality substantially.

def reduce_postcode(series):
    return series.str.split().str[0]

for frame in (X_train, X_val, X_test):
    frame["postcode_area"] = reduce_postcode(frame["postcode"])
    frame.drop(columns=["postcode"], inplace=True)

print(X_train["postcode_area"].nunique())   # 24

24

24 areas is manageable. One-hot encoding produces 24 binary columns, each with around 30 training observations on average — enough to estimate a coefficient from.

Step 7 — Build the pipeline

The pipeline encapsulates every transformation so that the correct fit/transform discipline is enforced automatically.

from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder, OrdinalEncoder
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LinearRegression

numeric_cols = [
    "floor_area_m2", "bedrooms", "bathrooms",
    "age_years", "distance_to_centre_km",
]

categorical_cols = ["property_type", "postcode_area", "has_garden"]

ordinal_cols = ["epc_rating"]     # E < D < C < B < A

numeric_pipe = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("scale",  StandardScaler()),
])

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

ordinal_pipe = Pipeline([
    ("impute", SimpleImputer(strategy="most_frequent")),
    ("ordinal", OrdinalEncoder(
        categories=[["E", "D", "C", "B", "A"]]
    )),
])

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

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

Three details worth noting.

• Median imputation rather than mean, because floor area is mildly skewed and the median is more robust to the very large houses.
• drop="first" in the one-hot encoder. Without it, the dummy variables plus the intercept are perfectly collinear — the dummy variable trap from Tutorial 3. Dropping one level per categorical makes XᵀX invertible.
• Ordered encoding for epc_rating. EPC ratings are ordinal — A is better than B — so treating them as ordered integers preserves that information. One-hot encoding would throw it away.

Step 8 — Fit and evaluate on validation

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

def report(name, y_true, y_pred):
    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)
    ratio = rmse / mae
    print(f"{name:<28} RMSE £{rmse:>9,.0f}   MAE £{mae:>9,.0f}   "
          f"R² {r2:>6.3f}   RMSE/MAE {ratio:.2f}")

# --- baseline ---
baseline_pred = np.full(len(y_val), y_train.mean())
report("mean baseline", y_val, baseline_pred)

# --- model, raw target ---
model.fit(X_train, y_train)
report("linear (raw price)", y_val, model.predict(X_val))

# --- model, log target ---
model_log = Pipeline([
    ("preprocess", preprocess),
    ("regressor",  LinearRegression()),
])

# log transform requires a wrapper so the prediction is exp'd back
from sklearn.compose import TransformedTargetRegressor

model_log = Pipeline([
    ("preprocess", preprocess),
    ("regressor",  TransformedTargetRegressor(
        regressor=LinearRegression(),
        func=np.log, inverse_func=np.exp,
    )),
])

model_log.fit(X_train, y_train)
report("linear (log price)", y_val, model_log.predict(X_val))

mean baseline                RMSE £   94,207   MAE £  71,438   R²  0.000   RMSE/MAE 1.32
linear (raw price)           RMSE £   46,318   MAE £  34,905   R²  0.758   RMSE/MAE 1.33
linear (log price)           RMSE £   44,061   MAE £  33,112   R²  0.781   RMSE/MAE 1.33

The log transform improves the validation RMSE by about 5% and R² by two points. That is a real but modest gain, and it comes with a change in what the model is optimising — it now minimises proportional error rather than absolute error. Given that a seller's intuition about accuracy is proportional ("within about 10%"), the log model is the better choice even on the raw metrics.

Both models comfortably beat the mean baseline, which had an RMSE of £94,207. That is the first question to answer and it is answered clearly.

Step 9 — Diagnose

import matplotlib.pyplot as plt

y_pred_val = model_log.predict(X_val)
resid = y_val - y_pred_val

fig, ax = plt.subplots(1, 2, figsize=(12, 4))

# residual plot against fitted values
ax[0].axhline(0, color="grey", linestyle="--", linewidth=1)
ax[0].scatter(y_pred_val, resid, alpha=0.4, s=12)
ax[0].set_xlabel("Fitted sale price (£)")
ax[0].set_ylabel("Residual (£)")
ax[0].set_title("Residuals vs fitted")

# residual plot against floor area — the most important feature
ax[1].axhline(0, color="grey", linestyle="--", linewidth=1)
ax[1].scatter(X_val["floor_area_m2"], resid, alpha=0.4, s=12)
ax[1].set_xlabel("Floor area (m²)")
ax[1].set_ylabel("Residual (£)")
ax[1].set_title("Residuals vs floor area")

What the plots show:

• Residuals vs fitted. The cloud is broadly shapeless, centred on zero. There is a mild widening at higher fitted values — heteroscedasticity — but no curvature. The model has captured the non-linear structure.
• Residuals vs floor area. No visible curvature, confirming that the relationship between size and price is roughly linear once the other features are included. There is more spread for larger properties.

The widening is worth investigating. It suggests the model's errors are proportionally larger for expensive properties, which is consistent with the raw scale results — fitting on the log target should have reduced it, but it has not eliminated it.

from sklearn.model_selection import learning_curve

train_sizes, train_scores, val_scores = learning_curve(
    model_log, X_train, y_train,
    train_sizes=np.linspace(0.1, 1.0, 8),
    scoring="neg_root_mean_squared_error",
    cv=5, random_state=42,
)

train_rmse = -train_scores.mean(axis=1)
val_rmse   = -val_scores.mean(axis=1)

for size, tr, va in zip(train_sizes, train_rmse, val_rmse):
    print(f"n={size:>4}   train £{tr:>9,.0f}   val £{va:>9,.0f}   gap £{va-tr:>9,.0f}")

n=  72   train £   31,904   val £   56,182   gap £   24,278
n= 144   train £   37,806   val £   51,240   gap £   13,434
n= 216   train £   40,913   val £   48,891   gap £    7,978
n= 288   train £   42,540   val £   47,502   gap £    4,962
n= 360   train £   43,388   val £   46,751   gap £    3,363
n= 432   train £   44,050   val £   46,255   gap £    2,205
n= 504   train £   44,372   val £   45,945   gap £    1,573
n= 720   train £   44,612   val £   44,061   gap £   -551

The training and validation curves converge as the training set grows, and the gap has essentially closed by the full training size. The validation error is still declining slowly, but the flattening is clear.

The diagnosis is that the model is well fitted but on the edge of high bias. More data would help a little — the validation curve has not quite plateaued — but the model has essentially extracted what these features can provide.

Step 10 — Evaluate on the test set, once

y_pred_test = model_log.predict(X_test)

report("test set — log model", y_test, y_pred_test)

# Also report the baseline on the test set for comparison
baseline_test = np.full(len(y_test), y_train.mean())
report("test set — baseline", y_test, baseline_test)

test set — log model       RMSE £   45,120   MAE £  33,806   R²  0.772   RMSE/MAE 1.33
test set — baseline        RMSE £   94,318   MAE £  71,502   R² -0.001   RMSE/MAE 1.32

The test performance is essentially identical to the validation performance. That is the good outcome: it means every decision made during validation generalised. A large gap in either direction would have been a warning — higher test error would mean we overfitted the validation set, and lower test error would suggest the validation split was unlucky.

Step 11 — Read the coefficients

# Extract the fitted coefficients from the pipeline
regressor   = model_log.named_steps["regressor"].regressor_
preprocess  = model_log.named_steps["preprocess"]

feature_names = preprocess.get_feature_names_out()
coefs = pd.DataFrame({
    "feature":    feature_names,
    "coef_log":   regressor.coef_,
})

# Convert log-scale coefficients to a percentage price change
coefs["pct_change_per_unit"] = (np.exp(coefs["coef_log"]) - 1) * 100

print(coefs.sort_values("coef_log", key=abs, ascending=False).to_string(index=False))

                feature  coef_log  pct_change_per_unit
   num__distance_to_centre_km   -0.241                -21.4
            num__floor_area_m2    0.198                 21.9
                num__age_years   -0.071                 -6.9
               ord__epc_rating    0.058                  6.0
             num__bathrooms     0.034                  3.5
              num__bedrooms     0.021                  2.1
     cat__property_type_semi   -0.134                -12.5
cat__property_type_terraced   -0.268                -23.5
     cat__property_type_flat   -0.412                -33.8
        cat__has_garden_True    0.043                  4.4
      cat__postcode_area_SW1    0.187                 20.6
      cat__postcode_area_N1     0.142                 15.3
                                     ...

Three cautions for the report.

• Scaled features give standardised coefficients. Because the numeric features were passed through StandardScaler, the coefficient on floor area is per standard deviation, not per square metre. For the report we would rescale it — divide by the training-set standard deviation — to give a per-unit figure that a reader can act on.
• Bedrooms and bathrooms have small coefficients. That is not because bedrooms do not matter; it is because they are correlated with floor area, and once floor area is in the model they add little. Their coefficients describe a comparison that rarely occurs in the data: two properties with identical floor area, location and age but different bedroom counts.
• None of these are causal. The distance coefficient does not mean "moving a house one kilometre closer to the centre would raise its price by 21%". It means "in this dataset, houses one kilometre closer to the centre sold for about 21% more, holding the other measured features constant."