Cross-Validation From Scratch and a Surprise at n=100

By Ken Koon Wong | August 5, 2026 Tags: #R #crossvalidation #resampling #bias #variance

According to standard textbooks, Leave-One-Out Cross-Validation (LOOCV) is characterized by having the lowest bias but the highest variance when compared to 10-fold or 5-fold CV. To put this theory to the test, I decided to program a K-Fold CV algorithm from the ground up using simulated data. At $n=1000$, the theory seems to hold.

Fun Side Note: I fed the text of this post into an AI image generator. While it used to produce gibberish and spelling errors, it now generates images with cohesive, readable text!

---

🎯 Motivations

In the modern era of Machine Learning, cross-validation is so fundamental that it's baked into libraries like sklearn or tidymodels. We rarely need to build it ourselves. However, I believe the most effective way to truly grasp a concept (besides just learning the concept 🤣) is to implement it from scratch.

Understanding K-Fold CV

In a K-Fold CV setup, the training data is partitioned into $K$ equal segments (folds). The model is trained $K$ times; in each iteration, a different fold is held out for validation while the remaining $K-1$ folds are used for training. The final performance metric is the average across all $K$ iterations, which provides a more stable estimate of model performance.

LOOCV is simply a special case where $K = n$ (every single observation acts as its own validation set).

The Textbook Hypothesis

I wanted to verify the following claims regarding the relationship between the number of folds, bias, and variance:

• Bias: $\text{LOOCV} < \text{10-fold} < \text{5-fold}$
• Variance: $\text{LOOCV} > \text{10-fold} > \text{5-fold}$

Additionally, I questioned whether repeating the process actually succeeds in reducing variance.

---

🛠️ Objectives

• Simulate a dataset with a known Data Generating Process (DGP).
• Build a K-Fold CV function from scratch.
• Use RMSE (Root Mean Square Error) to assess performance.
• Compare various candidate models via CV RMSE.
• Validate the winning model on a separate test set.

---

📊 Simulating the Data

Using the tidyverse library, I generated a dataset where the response variable $y$ is derived from a specific non-linear formula.

The Data Generating Process (DGP): $y = 0.5x^2 - 0.5w + 0.3wx + \epsilon$

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)

# Split: 80% Training, 20% Testing
idx <- sample(1:n, size = 0.8 * n)
train <- df[idx, ]
test <- df[-idx, ]

To visualize the relationship, I categorized $w$ into intervals:

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 = FALSE)

The visualization confirms that the relationship is decidedly non-linear.

---

💻 K-Fold Cross-Validation Implementation

Here is the logic for the manual CV process:

The Code

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 <- numeric(segment_portion * folds)
  y_log <- numeric(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)
}

---

📈 Assessing RMSE

Since we are predicting a continuous variable, the Root Mean Square Error (RMSE) serves as our loss function.

cv_log |> 
  group_by(formula) |> 
  summarize(rmse = sqrt(mean((y - predict)^2))) |> 
  arrange(rmse) |> 
  mutate(rmse = format(rmse, digits = 8))

Results Table

The model with the lowest RMSE is y ~ I(x^2) + w + w:x, which perfectly matches the original DGP. 5-Fold CV successfully recovered the underlying process.

---

🔍 The Deep Dive: Bias vs. Variance

While 5-fold worked, does the choice of $K$ actually matter? In previous studies, we could measure bias/variance using a true Average Treatment Effect (ATE), but finding a "true RMSE" is trickier.

To test the textbook claim, I collaborated with Claude Sonnet 3.5 to run a larger simulation:

1. Ground Truth: I created a massive test set ($n=10,000$). By the Law of Large Numbers, the RMSE on this set is treated as the "True RMSE."
2. Simulation: I ran the CV process across 500 different simulated datasets.
3. Metrics:
   • Bias: The average difference between the CV estimate and the True RMSE.
   • Variance: The spread of CV estimates across the 500 samples.

The Goal: Verify if $\text{Bias}_{\text{LOOCV}} < \text{Bias}_{\text{10-fold}} < \text{Bias}_{\text{5-fold}}$ and $\text{Var}_{\text{LOOCV}} > \text{Var}_{\text{10-fold}} > \text{Var}_{\text{5-fold}}$.