Cross validation on a single model (not model comparison)

Question

I understand the method of cross validation to be to leave out some part of a dataset (whether that be one data point at a time = LOO, or subsets = K fold), and train the model on some data, test the model's predictive accuracy with the remaining data, and repeat.

This 'method' then should tell you how well a model predicts out of sample, yet I only seem to see folks use it to compare models (ask which model does a better job predicting out of sample), by comparing these 'relative' model scores such as ELPD, LOOIC, WAIC (https://cran.r-project.org/web/packages/loo/vignettes/loo2-example.html, https://avehtari.github.io/modelselection/CV-FAQ.html).

It seems like one way to see if the model does a decent job at predicting, is to compare the model scores of a model with half the data to that of the other half (e.g. in Rloo_compare(loo(firsthalfmodel),loo(secondhalfmodel))), but that seems like cross validation within cross validation, since functions like loo are supposed to be doing cross validation themselves.

Is there some way that I can make a statement about a single model without comparing it to another with LOO or K-fold CV?

If I can assess one model, ELPD is often an output from LOO (see example R code below), but its interpretation doesn't make sense to me outside of a model comparison example.

From: https://avehtari.github.io/modelselection/CV-FAQ.html

"ELPD: The theoretical expected log pointwise predictive density for a new observations"

So this somehow tells me how predictive my model is, but I don't understand the implications of the numbers that come from an output, and I cannot seem to find this information anywhere - aside from a model comparison context.

An R example:

library(rstanarm)
set.seed(707)

dat<-data.frame(x = rnorm(1000),
                y = 0.5 + x*.2
)

mod1<-stan_glm(y ~ x,data=dat)
loo(mod1)


Computed from 4000 by 1000 log-likelihood matrix

         Estimate   SE
elpd_loo    172.7 22.8
p_loo         3.2  0.2
looic      -345.4 45.6
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

From: https://cran.r-project.org/web/packages/loo/vignettes/loo2-example.html

"If we had a well-specified model we would expect the estimated effective number of parameters (p_loo) to be smaller than or similar to the total number of parameters in the model."

Here p_loo is over 3, which is more than 3 times the total number of parameters in the model (only x). I am guessing that this would indicate model misspecification, how much to worry however seems elusive.

Can someone give, in layperson's terms, what these other numbers are referring to - in a non-model comparison context. What can I say about this model, given this output? For example, is the model a good fit to the data? Does it do a good job predicting? How would one go about making a statement about how well this model performed? Is there any overfitting going on? Or perhaps this method doesn't answer any of these questions?

I have skimmed this resource: https://arxiv.org/pdf/1507.04544.pdf which is dense, yet it all seems to be over my head, so I am looking for an answer that you might give to your grandparent I suppose.

Answer 1

Cross validation provides a point estimator of the recognition (or error) rate, and thus does say something about a single model. The question remains, however, how good this estimator is or, more precisely, how a confidence interval can be estimated for a recognition rate estimated by cross validation.

You can combine $n$-fold cross validation with the jackknife estimator for the variance, which cyclically omits one sample i, estimates the observable on the remining samples as $\theta_{(i)}$, and then computes the standard deviation as $$\sigma_{JK}(\hat{\theta}) = \sqrt{\frac{n-1}{n}\sum_{i=1}^n (\theta_{(i)}-\theta_{(.)})^2} \quad\mbox{ with } \quad \theta_{(.)}=\frac{1}{n}\sum_{i=1}^n\theta_{(i)}$$ Don't get confused that samples are cyclically ommitted twice: once in the jackknife procedure and then inside it in the LOO error rate estimation.

Out of curiosity, I have tried this out on the Iris dataset with a multivariate Gaussian Bayes classifier (the R function is called "qda", but mathematically this is equivalent):

library(MASS)

n <- nrow(iris)
rate <- rep(0,n) # memory preallocation

# compute LOO recognition rate for each left out sample
for (i in 1:n) {
    qda.result <- qda(iris[-i,-5], iris[-i,5], prior=rep(1/3,3), CV=TRUE)
    rate[i] <- mean(qda.result$class == iris$Species[-i])
}

# compute mean and jackknife variance
rate.m <- mean(rate)
sigma.jk <- sqrt( ((n-1)/n) * sum((rate - rate.m)^2) )
cat(sprintf("Jackknife LOO recognition rate: %f +/- %f\n", rate.m, sigma.jk))

This yields:

Jackknife LOO recognition rate: 0.973199 +/- 0.011573

I do not know, however, how good the coverage probability of confidence intervals based on $\sigma_{JK}$ is in this particular case, because leave-one-out has been used both for computing each recognition rate estimator rate[i] and for estimating the variance therefrom. I have a gut feeling that this looses some "degrees of freedom", and the variance estimator might be somewhat too small.

Maybe someone knows theoretical results about the statistical properties of this approach?

Answer 2

 This 'method' then should tell you how well a model predicts out of sample

Yes, and that is its original purpose, and also part of why it is called cross validation.

Depending on the figure of merit you use, that is in itself a proper statement about about a certain kind of generalization error.

In my field (chemometrics), e.g. $RMSE_{CV}$ (root mean squared error estimated via cross validation) is widely used as estimate of the model's performance.

---

IMHO, there are a few common but rather unfortunate misunderstandings around cross validation and machine learning.

One of them is "cross validation provides model optimzation" - which is IMHO better explained by saying that

1. cross validation is a scheme to generate tests to measure generalization performance of a given model, and
2. that information can then be used as target function for model optimization.

(Or alternatively, as generalization performance estimate)