Exercise: choose fairly, judge once

Session 9 — Applied Data Science 2026

Author
Affiliation

Claudius Gräbner-Radkowitsch

EUF and JKU

Published

23 11 2026

NoteThe task

Build a model that predicts the price of a house in Ames from a fitting sample of 300 sales — a different 300 from the live coding — and say, honestly, how wrong it will be on houses it has never seen.

Alone or in pairs. Two parts: Part A and Part B.

ImportantStart here

One new package:

install.packages("glmnet")

Download ames.csv from the session page into data/, create session09.qmd, and start with this setup chunk. It does the data preparation from the live coding — read it, you will be asked about it. The only difference: your 300 sales are drawn with set.seed(8).

```{r}
#| label: setup
#| include: false
library(tidyverse)
library(broom)
library(glmnet)
library(here)

i_am("session09.qmd")

ames <- read_csv(here("data", "ames.csv"))

# Three parts: holdout (used once, at the end), fitting sample, rest unused
set.seed(2026)
holdout_rows <- sample(nrow(ames), round(0.2 * nrow(ames)))
rest_rows <- setdiff(seq_len(nrow(ames)), holdout_rows)
set.seed(8)
fit_rows <- sample(rest_rows, 300)

# A missing lot frontage is stored as 0: fill in the median of the sales
# outside the holdout, and keep a flag
fill_value <- ames |>
  slice(rest_rows) |>
  filter(Lot_Frontage > 0) |>
  pull(Lot_Frontage) |>
  median()

houses <- ames |>
  mutate(
    across(where(is.character),          # categories in a fixed order
           \(x) factor(x, levels = sort(unique(x), method = "radix"))),
    flag_lot_frontage = as.integer(Lot_Frontage == 0),
    Lot_Frontage = if_else(Lot_Frontage == 0, fill_value, Lot_Frontage),
    ln_price = log(Sale_Price)
  ) |>
  select(-Sale_Price)

# One matrix, built once for all houses
X <- model.matrix(ln_price ~ ., data = houses)[, -1]
colnames(X) <- make.names(colnames(X))
# drop columns fewer than 10 houses outside the holdout have
X <- X[, colSums(X[rest_rows, ] != 0) >= 10]

dat <- as_tibble(X) |>
  mutate(ln_price = houses$ln_price)

holdout  <- dat[holdout_rows, ]
sales300 <- dat[fit_rows, ]

# A fold number from 1 to 5 for every sale of the fitting sample
set.seed(1)
sales300 <- sales300 |>
  mutate(fold = sample(rep(1:5, length.out = n())))

rmse <- function(actual, predicted) sqrt(mean((actual - predicted)^2))
```

Part A — Cross-validation

A1 · The flag and the fill value

  1. In two sentences: what does flag_lot_frontage record, and why not simply leave the zeros in Lot_Frontage? Use ames and a group_by() with summarise() to show that the houses with a zero differ from the others.
  2. fill_value is computed from rest_rows only. Why not from all houses?

A2 · Cross-validate three models

Write the function from the live coding:

cv_rmse <- function(formula) {
  purrr::map_dbl(1:5, \(k) {
    fit_data  <- sales300 |> filter(fold != k) |> select(-fold)
    test_data <- sales300 |> filter(fold == k)
    fit <- lm(formula, data = fit_data)
    rmse(pull(test_data, ln_price),
         predict(fit, newdata = test_data))
  })
}

Run it for these three models and report, for each, the mean of the five test errors and the smallest and largest of them:

Model Formula
m1 ln_price ~ Gr_Liv_Area
m2 ln_price ~ Gr_Liv_Area + Year_Built + Total_Bsmt_SF + Garage_Cars
m3 ln_price ~ .

Then fit m3 on all 300 sales and compute its error on those 300 sales. Compare it with its cross-validated error. What has happened? How many of its coefficients are NA — and why can R not estimate them?

Leave out the fold number: lm(ln_price ~ ., data = select(sales300, -fold)). sum(is.na(coef(...))) counts the NAs.

A3 · Your own model

Build an m2b: m2 plus up to six columns you think should matter for the price of a house. colnames(sales300) lists them all — dummies included, for example NeighborhoodNorthridge_Heights or Central_AirY. Does cross-validation agree that your additions help?

A4 · One round is not enough

Suppose you had only looked at round 1 — fold 1 as the test set. Would you have ranked m1, m2 and m3 the same way as the mean of the five rounds? Is that true for every round?

Part B — Lasso and ridge

B1 · Let the lasso choose

X_300     <- sales300 |> select(-ln_price, -fold) |> as.matrix()
X_holdout <- holdout |> select(-ln_price) |> as.matrix()
y_300     <- pull(sales300, ln_price)

lasso_cv <- cv.glmnet(X_300, y_300, alpha = 1, foldid = sales300$fold)

ggplot(data = tidy(lasso_cv), mapping = aes(x = log(lambda), y = sqrt(estimate))) +
  geom_ribbon(mapping = aes(ymin = sqrt(conf.low), ymax = sqrt(conf.high)),
              fill = "grey85") +
  geom_line() +
  geom_vline(xintercept = log(c(lasso_cv$lambda.min, lasso_cv$lambda.1se)),
             linetype = "dashed") +
  labs(x = "log(λ)  →  stricter", y = "Cross-validated RMSE (log price)")
  1. How many predictors does the lasso keep at lambda.min, and how many at lambda.1se?
  2. The chart shows the two as dashed lines, and a grey band of ± one standard error around the curve. Explain in one or two sentences how lambda.1se is chosen, and why one might prefer it.
  3. What is the cross-validated RMSE at lambda.min? Compare it with your models from A2 and A3. Why is that comparison fair?

coef(lasso_cv, s = "lambda.1se") gives all coefficients; count the ones that are not zero with sum(... != 0) - 1 (minus the intercept). tidy(lasso_cv) and lasso_cv$cvm hold squared errors, hence the sqrt(): the CV RMSE at lambda.min is sqrt(min(lasso_cv$cvm)).

A quick look without ggplot: plot(lasso_cv) draws the same curve in base R — on the squared-error scale, with error bars instead of a band.

B2 · Read the survivors

Rank the survivors at lambda.1se by their effect per standard deviation:

sd_300 <- apply(X_300, 2, sd)

tidy(lasso_cv$glmnet.fit) |>
  filter(lambda == lasso_cv$lambda.1se, term != "(Intercept)") |>
  mutate(per_sd = estimate * sd_300[term]) |>
  arrange(desc(abs(per_sd)))
  1. What does a per_sd of 0.05 mean, in words? Why rank by per_sd and not by estimate?
  2. Which of your A3 choices did the lasso also keep? Pick one survivor you did not think of and write one sentence on why it might predict the price — and one sentence on why that does not mean it is worth that much.

B3 · Ridge

Fit the same with alpha = 0. What is its cross-validated RMSE? How many of its coefficients are exactly zero — and why those? Which would you rather explain to someone selling a house, lasso or ridge?

B4 · The verdict, once

Now — and only now — use the holdout. Compute the holdout RMSE of m2, m3 (both fitted on sales300), the lasso at lambda.min and ridge at lambda.min. Write one paragraph for someone who wants to sell a house: which model, how many predictors, and how wrong it typically is — as a factor (exp(rmse)) and in dollars for a $180,000 house.

For glmnet: predict(lasso_cv, X_holdout, s = "lambda.min"). For lm: predict(fit, newdata = holdout).

If you finish early

a) All the sales

Replace the 300 sales by all 2,344 that are not in the holdout: sales_all <- dat[rest_rows, ], give them new folds, and run B1, B3 and B4 again. How far apart are OLS with every column and the lasso now? Why?

b) Optional: Airbnb, where the penalty does not pay

The same tools on different data: airbnb-hackney.csv (on the session page) has 4,478 usable Airbnb listings in Hackney, London, priced in US dollars per night, with 62 predictors. This chunk prepares it — it sets the holdout aside first, then fills in missing review scores with the median of the other listings and keeps a flag, just like Lot_Frontage:

raw <- read_csv(here("data", "airbnb-hackney.csv"))

usable <- raw |>
  filter(!is.na(price), !is.na(n_beds), !is.na(n_bathrooms))

# The holdout first, as for Ames: the rare amenities and the medians below
# are computed from the other 80% only
set.seed(2026)
hold_bnb <- sample(nrow(usable), round(0.2 * nrow(usable)))
usable_work <- usable[-hold_bnb, ]

amenity_counts <- usable_work |> select(starts_with("d_")) |> colSums()
rare_amenities <- names(amenity_counts)[amenity_counts < 20]

gappy <- c("n_review_scores_rating", "n_reviews_per_month", "n_days_since",
           "p_host_response_rate", "usd_cleaning_fee")
fill <- usable_work |>
  summarise(across(all_of(gappy), \(x) median(x, na.rm = TRUE)))

airbnb <- usable |>
  select(-any_of(rare_amenities)) |>
  select(-usd_price_day, -neighbourhood_cleansed, -cancellation_policy,
         -room_type, -property_type, -f_neighbourhood_cleansed) |>
  mutate(
    across(all_of(gappy), \(x) as.integer(is.na(x)), .names = "flag_{.col}"),
    across(all_of(gappy), \(x) replace_na(x, fill[[cur_column()]])),
    ln_price = log(price)
  ) |>
  select(-price, -flag_n_days_since) |>
  mutate(across(where(is.character), factor))
  1. Why does a listing without a review score need a flag? Compare the median price of listings with and without one.
  2. Use hold_bnb from the chunk as the holdout, fit OLS with all predictors on the rest, and the lasso and ridge with cv.glmnet() (build the matrix with model.matrix(ln_price ~ ., data = airbnb)[, -1]). Compare the three on the holdout.
  3. Why do lasso and ridge not help here? Count observations per predictor, and compare with the 300 Ames sales.

c) The tidymodels way

The live coding document has an appendix that does the lasso with the tidymodels framework, including the median imputation done properly inside each fold. Install tidymodels and run it.