---
title: Cross-Validation (Lecture 10)
author: 36-462/662, Spring 2022
date: 17 February 222
output: slidy_presentation
bibliography: locusts.bib
---

\[
\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)}
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]}
\newcommand{\Risk}{r}
\newcommand{\EmpRisk}{\hat{r}}
\newcommand{\Loss}{\ell}
\newcommand{\OptimalStrategy}{\sigma}
\DeclareMathOperator*{\argmin}{argmin}
\newcommand{\ModelClass}{S}
\newcommand{\OptimalModel}{s^*}
\DeclareMathOperator{\tr}{tr}
\newcommand{\optimand}{\theta}
\newcommand{\optimum}{\optimand^*}
\newcommand{\OptimalParameter}{\optimand^*}
\newcommand{\OptDim}{{p}}
\newcommand{\ERM}{\hat{\optimand}}
\newcommand{\ObjFunc}{{M}}
\newcommand{\Hessian}{\mathbf{k}}
\]


```{r, include=FALSE}
# General set-up options
library(knitr)
opts_chunk$set(size="small",background="white", highlight=FALSE,
               cache=TRUE, autodep=TRUE,
               tidy=TRUE, tidy.opts=list(comment=FALSE), warning=FALSE, message=FALSE,
               echo=FALSE)
set.seed(2022-02-17)
```

## Reminders {.smaller}

- We want to learn rules (models, strategies) which will act well in the future
- $\Rightarrow$ We'd _like_ to minimize risk, $\Risk(\optimand) = \Expect{\Loss(Y, s(X;\optimand))} =$ expected loss on new data
- What we _can_ do is minimize empirical risk, $\EmpRisk(\optimand) = n^{-1}\sum_{i=1}^{n}{\Loss(y_i, s(x_i; \optimand))}$
- $\optimum minimizes $\Risk(\optimand)$, $\ERM$ minimizes $\EmpRisk$
- Since $\Risk \neq \EmpRisk$, $\optimum \neq \ERM$
- Since $\EmpRisk \rightarrow \Risk$, $\ERM \rightarrow \optimum$
- Since $\ERM \rightarrow \optimum$, $\Risk(\ERM) \rightarrow \Risk(\optimum)$


## Agenda for today: more direct estimates of risk

- We got mathematical results last time about the risk of estimated model
- We'll mostly look today at procedures which require less theory
    + How do we know those procedures work?  More theory!


## Asymptotic results on the risk of our trained model

- Two asymptotic approximations:
\begin{eqnarray}
\Expect{\Risk(\ERM)} & \approx & \Risk(\optimum) + \frac{1}{2n}\tr{\mathbf{j}\Hessian}\\
\Expect{\Risk(\ERM)} & \approx & \Expect{\EmpRisk(\ERM)} + \frac{1}{n}\tr{\mathbf{j} \Hessian}
\end{eqnarray}
    + See [lecture 9](https://www.stat.cmu.edu/~cshalizi/dm/22/lectures/09/lecture-09.html#(16)) for definitions & estimation methods for the matrices
- This lets us estimate risk from the training data:
\[
\EmpRisk(\ERM) + \frac{1}{n}\tr{\mathbf{j}\Hessian}
\]
is an unbiased estimate of $\Expect{\Risk(\ERM)}
    + And in fact $\Risk(\ERM) \rightarrow \Expect{\Risk(\ERM)}$
- The quantity $n^{-1}\tr{\mathbf{j}\Hessian}$ is therefore called the **optimism**
- Notice: When $\mathrm{dim}(\optimand) = \OptDim$, $\Hessian$ is $[\OptDim \times \OptDim]$, $\mathbf{j}$ is also $[\OptDim\times \OptDim]$, so $\tr{\mathbf{j}\Hessian}$ is a sum of $\OptDim$ terms
- $\therefore$ the optimism is $O(p/n)$, at least roughly speaking




## Why aren't we happy with using the optimism to estimate true risk?

- We relied on a lot of calculus
- The formulas only apply when the model has continuous parameters
    + Think of adjusting regression slopes (continuous) vs. choosing the shape of a tree or neural network, or the number of nearest neighbors to use
       * Sometimes we try to make dis-continuous problems continuous, like deciding if the slope on some interaction term should be $\neq 0$
- The formulas only apply when the _loss function_ is continuous
    + log loss is continuous, 0-1 loss is not
- The formulas require us to take derivatives
    + Some of us hate calculus and wish it would go away
    + Others just don't like the work
	

## "Why think, when you can do the experiment?"

- If we want to see how well our model will generalize to new data, _have it generalize to new data_

1. Wait for new data to appear, or go out to gather it
    + Actual science (not data science); slow, expensive, painful, does not play to our strengths as statisticians
2. Get someone else's data, see if the model generalizes
    + Who is this person, who magically has exactly the type of data we need?
3. See how well our model generalizes from one part of our data to another


## Out-of-sample prediction evaluation, in general

- Divide the data into a training and testing part
    + Or a training set and a validation set, etc.
- Estimate the model using _only_ the training set
    + If we're comparing multiple models, estimate all of them on the training set
- Find the model(s)'s average loss on the training set
- Optional: consider multiple training/testing splits, and average loss over testing sets

**Make sure the training and testing data are independent**









## We have a data set


```{r}
# Load a data set
library(gamair)
data(chicago)
# Create an additional column which gives the date
day.zero <- as.Date("1993-12-31")
chicago$date <- chicago$time + day.zero
```

Deaths per day in Chicago, from all causes except accidents, plus measurements of air pollution and temperature.

```{r}
kable(head(chicago))
```

## We fit models and want to know how well they do

```{r, echo=TRUE, eval=FALSE}
chicago.lm <- lm(death ~ pm10median + o3median + so2median + tmpd,
                 death=chicago)
```

- ... this could go on through many models, many models
- We want the best model
- "Best" = predicts best on new data
  + We don't want to presume any of them are true


## Sample splitting

- To see how well a model predicts new data, _use it to predict new data_
- Divide the data at random into two halves, **training** and **testing**
- Fit all the models on the training set
- Evaluate the predictions of the _fitted_ models on the testing set
- Pick the best


## Sample splitting {.smaller}

```{r, echo=TRUE}
chicago$split <- sample(c("train", "test"), size=nrow(chicago), replace=TRUE)
head(chicago$split)
```

```{r}
kable(head(chicago))
```


## Sample splitting {.smaller}

```{r, echo=TRUE}
kable(head(chicago[chicago$split=="train",]))
```

```{r, echo=TRUE}
kable(head(chicago[chicago$split=="test",]))
```

## Sample splitting {.smaller}

- Re-number the data so points $1, 2, \ldots n/2$ are in the training set
  and $n/2+1, \ldots n$ are in the testing set
- We fit our model models $m_1, m_2, \ldots m_r$ using the training set,
  to get $\hat{m}_1, \hat{m}_2, \ldots \hat{m}_r$
- Our sample-splititng scores are then
\[
SSS_j \equiv \frac{1}{n/2}\sum_{i=n/2+1}^{n}{\ell(y_i, \hat{m}_j(x_i))}
\]
- Because we made the split at random, $\hat{m}_j$ is independent of the testing data, so $SSS_j$ is an unbiased estimate of $\mathbb{E}\left[\ell(Y, \hat{m}_j(X))\right]$.
- Pick the model with the best sample-splitting score


## Drawbacks of sample splitting

- _One_ random split is noisy
  + Often not such a big issue
- Which model predicts best _should_ change with $n$
    + $\Risk(\ERM) = \Risk(\optimum) + \frac{1}{2n}\tr{\mathbf{j}\Hessian^{-1}}$ so $\Risk(\ERM) = \Risk(\optimum) + O(\OptDim/n)$
	+ Say model 2 is more complicated than model 1 ($\OptDim_2 > \OptDim_1$) but a better approximation to the truth, $\Risk(\optimum_1) > \Risk(\optimum_2)$
	+ For small enough $n$, $\Risk(\ERM_1) < \Risk(\ERM_2)$, but for for large enough $n$, $\Risk(\ERM_2) < \Risk(\ERM_1)$
- The optimal model at sample size $n/2$ isn't the same as the optimal model with sample size $n$
- Obvious solution to the sample-size-bias issue is to make the training set bigger than the testing set
    + ... but now the score on the testing is noisier


## Enter cross-validation {.smaller}

- Divide the data into $k$ equally-sized chunks
  + Say $k=3$ for demo purposes

```{r, echo=TRUE}
head(rep(1:3, length.out=nrow(chicago)))
chicago$fold <- sample(rep(1:3, length.out=nrow(chicago)))
head(chicago$fold)
```

```{r, echo=TRUE}
kable(head(chicago[chicago$fold==1,]))
```

```{r, echo=TRUE}
kable(head(chicago[chicago$fold==2,]))
```

```{r, echo=TRUE}
kable(head(chicago[chicago$fold==3,]))
```


## How to actually do cross-validation

- For each fold $f$ in $1:k$,
  + Set the data points assigned to fold $f$ aside as the testing set
  + Estimate the model(s) on the rest of the data, getting $\hat{m}_j^{(-f)}$
  + Evaluate the models on the testing set,
  \[
  CV^{(f)}_j = \frac{1}{n/k}\sum_{i \in \mathrm{fold}(f)}{\Loss(y_i, \hat{m}_j^{(-f)}(x_i))}
  \]
    * Again, an unbiased estimate of $\Expect{\Loss(Y, \hat{m}_{j}^{(-f)}(X))}$
- Average over folds,
\[
CV_j = \frac{1}{k}\sum_{f=1}^{k}{CV^{(f)}_j}
\]
- Take the model with the best CV score



## Why should this work?

- Sample size isn't $n/2$ but $\frac{k-1}{k}n$, which is a lot closer to $n$
  + The sample-size-bias is attenuated
- _Every_ data point gets used to evaluate the models, and every data point
  gets used exactly the same number of times
- We're averaging over multiple random splits which reduces some of the randomness
  + But not all; any two of the training sets have $n\frac{k-2}{k}$ of the data points in common, which creates correlations between $\hat{m}_j^{(-f_1)}$ and $\hat{m}_j^{(-f_2)}$





## How do we know this works?

- Could do some theory, but it's surprisingly hard...
- Instead, let's simulate!

- "Synthetic" (i.e., fake) data example, so we know what the real answer is, and can tell which methods do and do not work

```{r}
# 20 standard-Gaussian X's
x = rnorm(20)
# Quadratic Y's
y = 7*x^2 - 0.5*x + rnorm(20)
```

```{r}
# Initial plot of training data plus true regression curve
plot(x,y,
      main="Noisy quadratic truth, plus polynomial fits")
curve(7*x^2-0.5*x, col="grey", add=TRUE)

# Fit polynomials and add them to the plot
# Make a list of model formulae for polynomials of degree 0 to 9
poly.formulae <- c("y~1", paste("y ~ poly(x,", 1:9, ")",
                                sep=""))
  # because "y ~ poly(x,0)" doesn't work
poly.formulae <- sapply(poly.formulae, as.formula)
# Get evenly spaced points for pretty plotting models
df.plot <- data.frame(x=seq(min(x), max(x), length.out=200))
# Fit polynomials of order 0 to 9
  # Store them in a list
fitted.models <- list(length=length(poly.formulae))
for (model_index in 1:length(poly.formulae)) {
  fm <- lm(formula=poly.formulae[[model_index]])
  lines(df.plot$x, predict(fm,newdata=df.plot),
        lty=model_index)
  fitted.models[[model_index]] <- fm
}
legend("top", legend=paste("degree",0:9), lty=1:10, ncol=2, cex=0.5)
```

- Visually, it looks like more complicated models fit better
  + Here "more complicated" = "more parameters"


## In-sample error just goes down with the order of the polynomial

- Quantitatively, more complicated models _do_ have lower mean-squared error
  + They have to, because degree-$d$ polynomials are special cases of degree-$(d+1)$ polynomials, and we fit by minimizing MSE

```{r}
mse.q <- sapply(fitted.models, function(mdl)
  { mean(residuals(mdl)^2) })
plot(0:9,mse.q,type="b",xlab="polynomial degree",
     ylab="mean squared error",
     log="y",
     main="In-sample MSE for polynomials")
```


## The more complicated models don't generalize well

- Add a lot of new data points from the _same_ distribution (blue dots)

```{r}
# Plot the old curves with testing data
x.new = rnorm(2e4) #20k points
y.new = 7*x.new^2 - 0.5*x.new + rnorm(length(x.new))
# 20k data points is very "busy", 2k gets the message across
plot(x.new[1:2e3], y.new[1:2e3],
     xlab="x", ylab="y", pch=24, cex=0.1, col="blue",
     xlim=range(x.new),ylim=2*range(y.new),
     main="Noisy quadratic curve plus polynomials fit to the first 20 data points")
curve(7*x^2-0.5*x, col="grey", add=TRUE)
# Use predict() on new data to see what the models predict
newdf <- data.frame(x=seq(from=2*min(x.new), to=2*max(x.new), length.out=200))
for (model_index in 1:length(fitted.models)) {
  fm <- fitted.models[[model_index]]
  lines(newdf$x, predict(fm, newdf), lty=model_index)
}
points(x,y)
legend("top", legend=paste("degree",0:9), lty=1:10, ncol=2, cex=0.5)
```


## More complicated models generalize _very badly_

- Use all the new data points to get the risk (= expected loss) of the
polynomials we fit to those first 20 data points

```{r}
# Calculate an already-fitted model's MSE on new data
  # Inputs: the model object
  # Outputs: the MSE
  # ATTN: Hard-wired to use "x.new" and "y.new" --- can you improve this?
gmse <- function(mdl) {
  predictions <- predict(mdl, data.frame(x=x.new))
  residuals <- y.new - predictions
  return(mean(residuals^2))
}
gmse.q <- sapply(fitted.models, gmse)
plot(0:9, mse.q, type="b",
     xlab="polynomial degree",
     ylab="mean squared error",
     log="y",
     ylim=c(min(mse.q), max(gmse.q)),
     main="Out-of-sample MSEs of the polynomials fit to the first 20 data points")
lines(0:9,gmse.q,lty=2,col="blue")
points(0:9,gmse.q,pch=24,col="blue")
```

- Note: logarithmic scale of vertical axis!
- ATTN: Can you add a good legend?


## $R^2$ and $R^2_{adjusted}$ are useless for goodness-of-fit

```{r}
# Extract R^2 and adjusted R^2 from a fitted model
  # Inputs: the lm model object
  # Output: vector of length 2, R^2 first
extract.rsqd <- function(mdl) {
  c(summary(mdl)$r.squared, summary(mdl)$adj.r.squared)
}
rsqd.q <- sapply(fitted.models, extract.rsqd)
plot(0:9, rsqd.q[1,],type="b",
     xlab="polynomial degree",
     ylab=expression(R^2),
     ylim=c(0,1),
     main="R^2 and adjusted R^2 for polynomials fit to first 20 data points")
lines(0:9, rsqd.q[2,], type="b", lty="dashed")
legend("bottomright",
       legend=c(expression(R^2), expression(R[adj]^2)),
       lty=c("solid","dashed"))
```


## k-fold CV for linear models

```{r, echo=TRUE}
# General function to do k-fold CV for a bunch of linear models
  # Inputs: dataframe to fit all models on, list or vector of model formulae,
    # number of folds of cross-validation
  # Output: vector of cross-validated MSEs for the models
cv.lm <- function(data, formulae, nfolds=5) {
  # Strip data of NA rows
    # ATTN: Better to check whether NAs are in variables used by the models
  data <- na.omit(data)
  # Make sure the formulae have type "formula"
  formulae <- sapply(formulae, as.formula)
  # Extract the name of the response variable from each formula
    # ATTN: CV doesn't make a lot of sense unless these are all the same!
  responses <- sapply(formulae, response.name)
  names(responses) <- as.character(formulae)
  n <- nrow(data)
  # Assign each data point to a fold, at random
    # see ?sample for the effect of sample(x) on a vector x
  fold.labels <- sample(rep(1:nfolds, length.out=n))
  mses <- matrix(NA, nrow=nfolds, ncol=length(formulae))
  colnames <- as.character(formulae)
  # ATTN: Can you replace the double for() loop below by defining a new
  # function and then calling outer()?
  for (fold in 1:nfolds) {
    test.rows <- which(fold.labels == fold)
    train <- data[-test.rows,]
    test <- data[test.rows,]
    for (form in 1:length(formulae)) {
       # Fit the model on the training data
       current.model <- lm(formula=formulae[[form]],
                           data=train)
       # Generate predictions on the testing data
       predictions <- predict(current.model, newdata=test)
       # Get the responses on the testing data
       test.responses <- test[,responses[form]]
       test.errors <- test.responses - predictions
       mses[fold, form] <- mean(test.errors^2)
    }
  }
  return(colMeans(mses))
}

# Extract the name of the response variable from a regression formula
  # Presumes response is the left-most variable
  # ATTN: Can you write a more robust version using terms()?
  # Inputs: regression formula
  # Outputs: name of the response variable
response.name <- function(formula) {
  var.names <- all.vars(formula)
  return(var.names[1])
}
```

## How well does cross-validation work?

- Remember that our original data had 20 points, and we're seeing how  we generalize to 20,000 points from the same distribution

```{r}
# Make a little data frame out of our little data
little.df <- data.frame(x=x, y=y)
# CV for the polynomials (defaults to five-fold)
cv.q <- cv.lm(little.df, poly.formulae)
plot(0:9, mse.q, type="b",
     xlab="polynomial degree",
     ylab="mean squared error", log="y",
     ylim=c(min(mse.q), max(gmse.q)))
lines(0:9, gmse.q,lty=2, col="blue",type="b",pch=2)
lines(0:9, cv.q,lty=3, col="red",type="b", pch=3)
legend("topleft",
       legend=c("In-sample", "Generalization", "CV"),
       col=c("black","blue","red"), lty=1:3, pch=1:3)
```


## Leave-one-out-CV (LOOCV)

- Take the limit, $k=n$
- Each data point gets "left out" _once_, to be the testing set
- Fit the models with the remaining $n-1$ data points as the training set
- See how well all models predict the left-out point
- Average over all $n$ points


## Short-cut formula for LOOCV

- Ordinarily LOOCV is the most time-consuming sort of CV
  + Need to re-estimate the model $n$ times
- **IF** the model is a linear smoother, $\hat{m}(x_0) = \sum_{j=1}^{n}{w(x_0, x_j) y_j}$, then there is a short-cut [@Wahba-spline-models]:
\[
LOOCV = \frac{1}{n}\sum_{i=1}^{n}{\left(\frac{y_i - \hat{m}(x_i)}{1-w_{ii}}\right)^2}
\]
- We estimate the model _once_, and get the diagonal of its smoother/hat matrix $w_{ii} = w(x_i, x_i)$
  + We often have to work this out anyway


## $k$-fold CV vs. LOOCV

- LOOCV is unbiased as an estimate of risk at sample size $n-1$
    + So a pretty good estimate of the risk at sample size $n$
- LOOCV tends to have high variance
- If there is a true model, LOOCV tends to pick a model which is _bigger_ than the true model
- $k$-fold CV tends to be better at selecting the true model
    + It's estimating risk at the sample size $\frac{k-1}{k}n$, which puts more complex models at more of a disadvantage (compared to LOOCV)
- On the other hand, LOOCV is better _predictively_
  + That is, the risk of the model picked by LOOCV tends to be smaller than the risk of the model picked by $k$-fold
  + There are theoretical reasons to think we can't get _both_ predictive optimality _and_ consistent selection of the true model [@YYang-conflict-between-identification-and-estimation; @Claeskens-Hjort-model-selection]
- Advice:
  + In prediction tasks, _if time allows_, use LOOCV
  + Otherwise, 5- or 10- fold CV


## Reminder: _Best_ CV score is optimistically biased

- The CV score of our _best_ model is going to be optimistic about how well
  that model will generalize
  + For basically the same reason that in-sample error is optimistic
    + "The winner's curse" (see backup)
- The easiest way to get an unbiased estimate of how well the winner will
  perform is to evaluate it on held-out data, which we never use for anything
  else
- Alternately, @Tibshirani-Tibshirani-bias-correction-for-cv suggest a bias-correction formula
  

## Summing up

- In-sample error is going to be optimistically biased, often badly so
  + $R^2$ and $R^2_{adjusted}$ are useless
- The "optimism" gives a calculus-based way of estimating true risk (and so of selecting models)
    + Asymptotic, needs calculus, so may not apply, may need too much thinking
- Data splitting gives a direct estimate of true risk, but at a smaller sample size
     + Which we know from the theory leading to the optimism formula can really matter!
- Use cross-validation
  + Typically 5- or 10- fold
  + Leave-one-out if you can use the short-cut formula or $n$ is really small
- Remember the winner's CV score will still be an optimistic estimate  of true risk
  + Use a hold-out set if you really need to know the risk of the _selected_ model


## Backup: Further reading

- Data splitting and cross-validation go back in statistical _practice_ (as opposed to theory) for many decades
  + @Stone1974 and @Geisser-predictive-sample-reuse seem to be the first formal publications
  + @Geisser-Eddy-predictive-approach made the generality of the method clear
  + Another notable early contributor to cross-validation as a general method was Grace Wahba; her book @Wahba-spline-models covers this (with references to her original papers from the 1970s)
- Computer scientists adopted cross-validation from statisticians around 1990; this is was one of the steps that led to the modern discipline of "machine learning"
- @Arlot-Celisse-on-cv is a good recent review of cross-validation methods in practice
  + Cross-validation for _dependent_ data raises its own issues [@Burman-Chow-Nolan-cross-validation-for-dependent; @Racine-consistent-cv-for-dependent-data]
- Important references on the _theory_ of cross-validation include:
  + @Gyorfi-Kohler-Krzyzak-Walk-nonparametric-regrssion, chs. 7 and 8
  + @van-der-Laan-Dudoit-unified-cv-methodology
  + @van-der-Vaart-Dudoit-van-der-Laan-cv-oracle
  + These are all much more mathematically involved than this course...
  
  
## Backup: The Winner's Curse

- Enpirical risk = true risk plus noise, $\EmpRisk(\optimand) = \Risk(\optimand) + \gamma(\optimand)$
- $\Expect{\gamma(\optimand)} = 0$
- When we pick $\ERM = \argmin_{\optimand}{\EmpRisk(\optimand)}$, we
    + tend to pick a value of $\optimand$ with $\Risk$ close to the minimum
    + tend to pick a $\optimand$ where $\gamma(\optimand) < 0$
- If data splitting, CV, etc., give us an unbiased estimate of true risk, but we use those estimates to pick one model out of many, the score of the _selected_ model is a biased estimate of _its_ risk
    + Same argument as above!
- Analogy: the **winner's curse** in auctions [@Thaler-winners-curse]
    + Bids in auctions $=$ true economic value of the object _plus_ noise in estimating its value
    + The highest bidder tends to be someone for whom the object is very valuable
         * And has that much money...
    + The highest bidder _also_ tends to be someone with a lot of noise in their estimate of the value
     + If you win an auction, that is evidence that you over-paid, in and of itself
 - Estimation and model selection are like the winner's curse
     + Even if we use an unbiased score to select a parameter value or a model, winning the competition is evidence of optimism
 
 

## References
