Cross-Validation From Scratch and a Surprise at n=100 | Everyday Is A School Day
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. 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)<br>n 1000<br>x rnorm(n)<br>w rnorm(n)<br>y 0.5*x^2 + -0.5*w + 0.3*w*x + rnorm(n)<br>df tibble(x,y,w)<br>idx sample(1:n, size=0.8*n)<br>train df[idx, ]<br>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 |><br>mutate(w_cut = cut_interval(w, n=5)) |><br>ggplot(aes(x=x, y=y, color=w_cut, group=w_cut)) +<br>geom_point(alpha=0.5) +<br>theme_bw() +<br>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<br>segment_portion nrow(train)/folds<br>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"),<br>as.formula("y~w:x"),as.formula("y~w"),as.formula("y~x+w+x:w"),as.formula("y~I(x^2)+w:x"),<br>as.formula("y~I(x^2)+w"))
cv_log tibble()
for (formula in formula_list) {<br>print(formula)<br>predict_log y_log vector(mode="numeric",length=segment_portion*folds)<br>start 1<br>end segment_portion
for (fold in 1:folds) {<br>val_i train[start:end,]<br>train_i train[-c(start:end),]<br>model_i lm(formula,train_i)<br>predict_i predict(model_i, val_i)<br>predict_log[start:end] predict_i<br>y_log[start:end] val_i$y<br>start end + 1<br>end start + segment_portion - 1
val_df tibble(predict=predict_log,y=y_log) |><br>mutate(formula=deparse(formula))<br>cv_log cv_log |><br>bind_rows(val_df)
## y ~ x<br>## y ~ I(x^2)<br>## y ~ I(x^2) + w + w:x<br>## y ~ I(x^3) + w + w:x<br>## y ~ w:x<br>## y ~ w<br>## y ~ x + w + x:w<br>## y ~ I(x^2) + w:x<br>## y ~ I(x^2) + w<br>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 |><br>group_by(formula) |><br>summarize(rmse = sqrt(mean((y-predict)^2))) |><br>arrange(rmse) |><br>mutate(rmse = format(rmse, digits = 8))
## # A tibble: 9 × 2<br>## formula rmse<br>##<br>## 1 y ~ I(x^2) + w + w:x 1.0397338<br>## 2 y ~ I(x^2) + w 1.1008693<br>## 3 y ~ I(x^2) + w:x 1.1607225<br>## 4 y ~ I(x^2) 1.2144857<br>## 5 y ~ x + w + x:w 1.2826723<br>## 6 y ~ I(x^3) + w + w:x 1.2912046<br>## 7 y ~ w 1.3578304<br>## 8 y ~ w:x 1.3732478<br>## 9 y ~ x 1.4451807<br>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...