スクラッチからのクロスバリデーションとn=100での驚き

Textbooks say LOOCV has the lowest bias but highest variance compared to 10 and 5-fold. Coded a K-Fold CV from scratch for learning to test that on simulated data 🔍📊 — and at n=1000 it holds up. At n=100? Not so much. 🤔 The above image was generated via chatGPT. Uploaded all the text of this blog post and asked it to generate a cartoon. Very impressive! It used to be spelling error and gibberish of text in the past, but now cohesive words on image. Just wow. Motivations Crossvalidation is such a crucial step in Machine Learning (and traditional methods) that nowadays is incorporated in easy to use sklearn or tidymodels without us needing to build one from scratch. As with my other learning experience, the best way to learn the concept (other than learning the concept 🤣) is to code it from the ground up and see how it works! In K-Fold CV, the training data is split into K chunks; the model is trained K times, each time holding out a different chunk. Performance is averaged across all K folds, giving a more stable estimate. A special case is Leave-One-Out CV (LOOCV), where each individual observation serves as its own validation set. It’s thorough but computationally expensive. I was told that, bias LOOCV < 10-fold < 5-fold; whereas variance LOOCV > 10-fold > 5-fold. Is that true? Also, what’s with the repeats, does that really reduce variance? Let’s check them out. Objectives: Simulate data with a known data-generating process Implement K-Fold cross-validation from scratch Assessing RMSE Compare candidate models using CV RMSE Verify the best model on a held-out test set Opportunities For Improvement Lessons Learnt Simulate Data library(tidyverse) set.seed(1) n <- 1000 x <- rnorm(n) w <- rnorm(n) y <- 0.5*x^2 + -0.5*w + 0.3*w*x + rnorm(n) df <- tibble(x,y,w) idx <- sample(1:n, size=0.8*n) train <- df[idx, ] test <- df[-idx, ] The above code simulates a dataset with 1000 observations, where the response variable y is generated based on a known data-generating process involving predictors x and w. The dataset is then split into a training set (80%) and a test set (20%). Let’s visualize. df |> mutate(w_cut = cut_interval(w, n=5)) |> ggplot(aes(x=x, y=y, color=w_cut, group=w_cut)) + geom_point(alpha=0.5) + theme_bw() + geom_smooth(method = "gam", se=F) Wow, very interesting visualization where the relationships are definitely not linear here. It’s some form of interaction between x and w. Let’s see if we can recover the underlying data-generating process using K-Fold Cross-Validation. K-Fold Cross-Validation From Scratch folds <- 5 segment_portion <- nrow(train)/folds formula_list <- list(as.formula("y~x"),as.formula("y~I(x^2)"),as.formula("y~I(x^2)+w+w:x"),as.formula("y~I(x^3)+w+w:x"), as.formula("y~w:x"),as.formula("y~w"),as.formula("y~x+w+x:w"),as.formula("y~I(x^2)+w:x"), as.formula("y~I(x^2)+w")) cv_log <- tibble() for (formula in formula_list) { print(formula) predict_log <- y_log <- vector(mode="numeric",length=segment_portion*folds) start <- 1 end <- segment_portion for (fold in 1:folds) { val_i <- train[start:end,] train_i <- train[-c(start:end),] model_i <- lm(formula,train_i) predict_i <- predict(model_i, val_i) predict_log[start:end] <- predict_i y_log[start:end] <- val_i$y start <- end + 1 end <- start + segment_portion - 1 } val_df <- tibble(predict=predict_log,y=y_log) |> mutate(formula=deparse(formula)) cv_log <- cv_log |> bind_rows(val_df) } ## y ~ x ## y ~ I(x^2) ## y ~ I(x^2) + w + w:x ## y ~ I(x^3) + w + w:x ## y ~ w:x ## y ~ w ## y ~ x + w + x:w ## y ~ I(x^2) + w:x ## y ~ I(x^2) + w Alright, what we’ve done above is a manual implementation of K-Fold Cross-Validation. We loop through each formula in our list, and for each formula, we split the training data into 5 folds. For each fold, we train the model on the other 4 folds and validate it on the current fold. We store the predictions and actual values for later evaluation. We basically want to see which formula has the lowest RMSE across the folds. Let’s calculate that next. From the DGP formula, we know that the best model should be y~I(x^2)+w+w:x. Let’s see if we can recover that using K-Fold CV. Assessing RMSE cv_log |> group_by(formula) |> summarize(rmse = sqrt(mean((y-predict)^2))) |> arrange(rmse) |> mutate(rmse = format(rmse, digits = 8)) ## # A tibble: 9 × 2 ## formula rmse ## <chr> <chr> ## 1 y ~ I(x^2) + w + w:x 1.0397338 ## 2 y ~ I(x^2) + w 1.1008693 ## 3 y ~ I(x^2) + w:x 1.1607225 ## 4 y ~ I(x^2) 1.2144857 ## 5 y ~ x + w + x:w 1.2826723 ## 6 y ~ I(x^3) + w + w:x 1.2912046 ## 7 y ~ w 1.3578304 ## 8 y ~ w:x 1.3732478 ## 9 y ~ x 1.4451807 Here our loss function is RMSE since y is a continuous data and we’re trying to predict that. The formula with the lowest RMSE is indeed y~I(x^2)+w+w:x, which matches the underlying data-generating process. OK at least, right now we are able to recover the underlying DGP using 5-Fold Cross-Validation. But is there a difference between 5 fold, 10 fold, or even LOOCV? If there is a difference, how do we even assess that? In the past we were able to assess bias and variance based on a true ATE, but what on earth is a true RMSE !?! To check whether the textbook claim (bias LOOCV < 10-fold < 5-fold; variance LOOCV > 10-fold > 5-fold) holds up, we ran a small simulation with help from Claude Sonnet 5. Since we control the data-generating process, we can compare the “correct formula” (assuming the correct formula has the lowest RMSE as above) with 500 different simulated dataset against a “true” RMSE estimated from a large test set (n=10000) — large enough, by the law of large numbers, to treat as ground truth. Averaging across simulations gives bias (how far off CV runs from the true error) and variance (how much CV’s estimate swings from sample to sample) for each method. Is this legit? 🤔 If the textbook claim is correct, we should be able to observe bias LOOCV < 10-fold < 5-fold; variance LOOCV > 10-fold > 5-fold. Let’s see if we can observe that in the simulation below. Compare Candidate Models # set.seed(1) # k-fold CV RMSE for a given formula and dataset (k = n gives LOOCV) cv_rmse <- function(data, formula, k) { n <- nrow(data) folds <- sample(rep(1:k, length.out = n)) preds <- numeric(n) for (i in 1:k) { train_i <- data[folds != i, ] val_i <- data[folds == i, ] model_i <- lm(formula, train_i) preds[folds == i] <- predict(model_i, val_i) } sqrt(mean((data$y - preds)^2)) } # "true" RMSE: fit on train, evaluate on a large fresh draw from the DGP true_rmse <- function(train, formula, n_test = 10000) { x <- rnorm(n_test); w <- rnorm(n_test) y <- 0.5*x^2 - 0.5*w + 0.3*w*x + rnorm(n_test) test <- tibble(x, y, w) model <- lm(formula, train) sqrt(mean((test$y - predict(model, test))^2)) } formula_true <- as.formula("y ~ I(x^2) + w + w:x") n_sim <- 500 n_train <- 100 results <- vector("list", n_sim) for (s in 1:n_sim) { x <- rnorm(n_train); w <- rnorm(n_train) y <- 0.5*x^2 - 0.5*w + 0.3*w*x + rnorm(n_train) train_s <- tibble(x, y, w) results[[s]] <- tibble( sim = s, true_err = true_rmse(train_s, formula_true), loocv = cv_rmse(train_s, formula_true, k = n_train), cv5 = cv_rmse(train_s, formula_true, k = 5), cv10 = cv_rmse(train_s, formula_true, k = 10) ) } sim_df <- bind_rows(results) sim_long <- sim_df |> pivot_longer(cols = c(loocv, cv5, cv10), names_to = "method", values_to = "cv_estimate") sim_long |> group_by(method) |> summarize( mean_cv_estimate = mean(cv_estimate), mean_true_error = mean(true_err), bias = mean(cv_estimate - true_err), variance = var(cv_estimate), .groups = "drop" ) |> arrange(bias) |> mutate(variance = format(variance, digit = 8)) ## # A tibble: 3 × 5 ## method mean_cv_estimate mean_true_error bias variance ## <chr> <dbl> <dbl> <dbl> <chr> ## 1 loocv 1.00 1.00 0.000703 0.00053227867 ## 2 cv10 1.00 1.00 0.000945 0.00053192050 ## 3 cv5 1.00 1.00 0.00124 0.00053176441 Wow, looking at the results, we can see that the textbook claim holds up. LOOCV has the lowest bias, but the highest variance. 10-fold CV is in between, and 5