Examples
Polynomial regression for theAuto data
The data for this example are drawn from theISLR2package for R, associated withJames et al.(2021). The presentation here is close (though not identical) tothat in the original source(James et al., 2021,secs. 5.1, 5.3), and it demonstrates the use of thecv() function in thecv package.1
TheAuto dataset contains information about 392cars:
data("Auto",package="ISLR2")head(Auto)#> mpg cylinders displacement horsepower weight acceleration year origin#> 1 18 8 307 130 3504 12.0 70 1#> 2 15 8 350 165 3693 11.5 70 1#> 3 18 8 318 150 3436 11.0 70 1#> 4 16 8 304 150 3433 12.0 70 1#> 5 17 8 302 140 3449 10.5 70 1#> 6 15 8 429 198 4341 10.0 70 1#> name#> 1 chevrolet chevelle malibu#> 2 buick skylark 320#> 3 plymouth satellite#> 4 amc rebel sst#> 5 ford torino#> 6 ford galaxie 500dim(Auto)#> [1] 392 9With the exception oforigin (which we don’t use here),these variables are largely self-explanatory, except possibly for unitsof measurement: for details seehelp("Auto", package="ISLR2").
We’ll focus here on the relationship ofmpg (miles pergallon) tohorsepower, as displayed in the followingscatterplot:
mpg vshorsepower for theAutodata
The relationship between the two variables is monotone, decreasing,and nonlinear. FollowingJames et al.(2021), we’ll consider approximating the relationship by apolynomial regression, with the degree of the polynomial\(p\) ranging from 1 (a linear regression) to10.2Polynomial fits for\(p = 1\) to\(5\) are shown in the following figure:
plot(mpg~ horsepower,data = Auto)horsepower<-with(Auto,seq(min(horsepower),max(horsepower),length =1000))for (pin1:5) { m<-lm(mpg~poly(horsepower, p),data = Auto) mpg<-predict(m,newdata =data.frame(horsepower = horsepower))lines(horsepower, mpg,col = p+1,lty = p,lwd =2)}legend("topright",legend =1:5,col =2:6,lty =1:5,lwd =2,title ="Degree",inset =0.02)mpg vshorsepower for theAutodata
The linear fit is clearly inappropriate; the fits for\(p = 2\) (quadratic) through\(4\) are very similar; and the fit for\(p = 5\) may over-fit the data by chasingone or two relatively highmpg values at the right (but seethe CV results reported below).
The following graph shows two measures of estimated (squared) erroras a function of polynomial-regression degree: The mean-squared error(“MSE”), defined as\(\mathsf{MSE} =\frac{1}{n}\sum_{i=1}^n (y_i - \widehat{y}_i)^2\), and the usualresidual variance, defined as\(\widehat{\sigma}^2 = \frac{1}{n - p - 1}\sum_{i=1}^n (y_i - \widehat{y}_i)^2\). The former necessarilydeclines with\(p\) (or, more strictly,can’t increase with\(p\)), while thelatter gets slightly larger for the largest values of\(p\), with the “best” value, by a smallmargin, for\(p = 7\).
library("cv")# for mse() and other functionsvar<- mse<-numeric(10)for (pin1:10) { m<-lm(mpg~poly(horsepower, p),data = Auto) mse[p]<-mse(Auto$mpg,fitted(m)) var[p]<-summary(m)$sigma^2}plot(c(1,10),range(mse, var),type ="n",xlab ="Degree of polynomial, p",ylab ="Estimated Squared Error")lines(1:10, mse,lwd =2,lty =1,col =2,pch =16,type ="b")lines(1:10, var,lwd =2,lty =2,col =3,pch =17,type ="b")legend("topright",inset =0.02,legend =c(expression(hat(sigma)^2),"MSE"),lwd =2,lty =2:1,col =3:2,pch =17:16)Estimated squared error as a function of polynomial degree,\(p\)
The code for this graph uses themse() function from thecv package to compute the MSE for each fit.
Usingcv()
The genericcv() function has an"lm"method, which by default performs\(k =10\)-fold CV:
m.auto<-lm(mpg~poly(horsepower,2),data = Auto)summary(m.auto)#>#> Call:#> lm(formula = mpg ~ poly(horsepower, 2), data = Auto)#>#> Residuals:#> Min 1Q Median 3Q Max#> -14.714 -2.594 -0.086 2.287 15.896#>#> Coefficients:#> Estimate Std. Error t value Pr(>|t|)#> (Intercept) 23.446 0.221 106.1 <2e-16 ***#> poly(horsepower, 2)1 -120.138 4.374 -27.5 <2e-16 ***#> poly(horsepower, 2)2 44.090 4.374 10.1 <2e-16 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#>#> Residual standard error: 4.37 on 389 degrees of freedom#> Multiple R-squared: 0.688, Adjusted R-squared: 0.686#> F-statistic: 428 on 2 and 389 DF, p-value: <2e-16cv(m.auto)#> R RNG seed set to 860198#> cross-validation criterion (mse) = 19.605The"lm" method by default usesmse() asthe CV criterion and the Woodbury matrix identity(Hager, 1989) to update the regression with eachfold deleted without having literally to refit the model. (Computationaldetails are discussed in a separate vignette.) Theprint()method for"cv" objects simply print the cross-validationcriterion, here the MSE. We can usesummary() to obtainmore information (as we’ll do routinely in the sequel):
summary(cv.m.auto<-cv(m.auto))#> R RNG seed set to 393841#> 10-Fold Cross Validation#> method: Woodbury#> criterion: mse#> cross-validation criterion = 19.185#> bias-adjusted cross-validation criterion = 19.175#> full-sample criterion = 18.985The summary reports the CV estimate of MSE, a biased-adjustedestimate of the MSE (the bias adjustment is explained in the finalsection), and the MSE is also computed for the original, full-sampleregression. Because the division of the data into 10 folds is random,cv() explicitly (randomly) generates and saves a seed forR’s pseudo-random number generator, to make the results replicable. Theuser can also specify the seed directly via theseedargument tocv().
Theplot() method for"cv" objects graphsthe CV criterion (here MSE) by fold or the coefficients estimates witheach fold deleted:
CV criterion (MSE) for cases in each fold.
Regression coefficients with each fold omitted.
To perform LOO CV, we can set thek argument tocv() to the number of cases in the data, herek=392, or, more conveniently, tok="loo" ork="n":
summary(cv(m.auto,k ="loo"))#> n-Fold Cross Validation#> method: hatvalues#> criterion: mse#> cross-validation criterion = 19.248For LOO CV of a linear model,cv() by default uses thehatvalues from the model fit to the full data for the LOO updates, andreports only the CV estimate of MSE. Alternative methods are to use theWoodbury matrix identity or the “naive” approach of literally refittingthe model with each case omitted. All three methods produce exactresults for a linear model (within the precision of floating-pointcomputations):
summary(cv(m.auto,k ="loo",method ="naive"))#> n-Fold Cross Validation#> method: naive#> criterion: mse#> cross-validation criterion = 19.248#> bias-adjusted cross-validation criterion = 19.248#> full-sample criterion = 18.985summary(cv(m.auto,k ="loo",method ="Woodbury"))#> n-Fold Cross Validation#> method: Woodbury#> criterion: mse#> cross-validation criterion = 19.248#> bias-adjusted cross-validation criterion = 19.248#> full-sample criterion = 18.985The"naive" and"Woodbury" methods alsoreturn the bias-adjusted estimate of MSE and the full-sample MSE, butbias isn’t an issue for LOO CV.
Comparing competing models
Thecv() function also has a method that can be appliedto a list of regression models for the same data, composed using themodels() function. For\(k\)-fold CV, the same folds are used forthe competing models, which reduces random error in their comparison.This result can also be obtained by specifying a common seed for R’srandom-number generator while applyingcv() separately toeach model, but employing a list of models is more convenient for both\(k\)-fold and LOO CV (where there isno random component to the composition of the\(n\) folds).
We illustrate with the polynomial regression models of varying degreefor theAuto data (discussed previously), beginning byfitting and saving the 10 models:
for (pin1:10) { command<-paste0("m.", p,"<- lm(mpg ~ poly(horsepower, ", p,"), data=Auto)")eval(parse(text = command))}objects(pattern ="m\\.[0-9]")#> [1] "m.1" "m.10" "m.2" "m.3" "m.4" "m.5" "m.6" "m.7" "m.8" "m.9"summary(m.2)# for example, the quadratic fit#>#> Call:#> lm(formula = mpg ~ poly(horsepower, 2), data = Auto)#>#> Residuals:#> Min 1Q Median 3Q Max#> -14.714 -2.594 -0.086 2.287 15.896#>#> Coefficients:#> Estimate Std. Error t value Pr(>|t|)#> (Intercept) 23.446 0.221 106.1 <2e-16 ***#> poly(horsepower, 2)1 -120.138 4.374 -27.5 <2e-16 ***#> poly(horsepower, 2)2 44.090 4.374 10.1 <2e-16 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#>#> Residual standard error: 4.37 on 389 degrees of freedom#> Multiple R-squared: 0.688, Adjusted R-squared: 0.686#> F-statistic: 428 on 2 and 389 DF, p-value: <2e-16The convoluted code within the loop to produce the 10 models insuresthat the model formulas are of the form, e.g.,mpg ~ poly(horsepower, 2) rather thanmpg ~ poly(horsepower, p), which may not work properly (seebelow).
We then applycv() to the list of 10 models (thedata argument is required):
# 10-fold CVcv.auto.10<-cv(models(m.1, m.2, m.3, m.4, m.5, m.6, m.7, m.8, m.9, m.10),data = Auto,seed =2120)cv.auto.10[1:2]# for the linear and quadratic models#> Model model.1:#> cross-validation criterion = 24.246#>#> Model model.2:#> cross-validation criterion = 19.346summary(cv.auto.10)#>#> Model model.1:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 24.246#> bias-adjusted cross-validation criterion = 24.23#> full-sample criterion = 23.944#>#> Model model.2:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 19.346#> bias-adjusted cross-validation criterion = 19.327#> full-sample criterion = 18.985#>#> Model model.3:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 19.396#> bias-adjusted cross-validation criterion = 19.373#> full-sample criterion = 18.945#>#> Model model.4:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 19.444#> bias-adjusted cross-validation criterion = 19.414#> full-sample criterion = 18.876#>#> Model model.5:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 19.014#> bias-adjusted cross-validation criterion = 18.983#> full-sample criterion = 18.427#>#> Model model.6:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 18.954#> bias-adjusted cross-validation criterion = 18.916#> full-sample criterion = 18.241#>#> Model model.7:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 18.898#> bias-adjusted cross-validation criterion = 18.854#> full-sample criterion = 18.078#>#> Model model.8:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 19.126#> bias-adjusted cross-validation criterion = 19.068#> full-sample criterion = 18.066#>#> Model model.9:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 19.342#> bias-adjusted cross-validation criterion = 19.269#> full-sample criterion = 18.027#>#> Model model.10:#> 10-Fold Cross Validation#> method: Woodbury#> cross-validation criterion = 20.012#> bias-adjusted cross-validation criterion = 19.882#> full-sample criterion = 18.01# LOO CVcv.auto.loo<-cv(models(m.1, m.2, m.3, m.4, m.5, m.6, m.7, m.8, m.9, m.10),data = Auto,k ="loo")cv.auto.loo[1:2]# linear and quadratic models#> Model model.1:#> cross-validation criterion = 24.232#>#> Model model.2:#> cross-validation criterion = 19.248summary(cv.auto.loo)#>#> Model model.1:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 24.232#>#> Model model.2:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 19.248#>#> Model model.3:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 19.335#>#> Model model.4:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 19.424#>#> Model model.5:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 19.033#>#> Model model.6:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 18.979#>#> Model model.7:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 18.833#>#> Model model.8:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 18.961#>#> Model model.9:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 19.069#>#> Model model.10:#> n-Fold Cross Validation#> method: hatvalues#> cross-validation criterion = 19.491Because we didn’t supply names for the models in the calls to themodels() function, the namesmodel.1,model.2, etc., are generated by the function.
Finally, we extract and graph the adjusted MSEs for\(10\)-fold CV and the MSEs for LOO CV (seethe section below on manipulating"cv" objects:
cv.mse.10<-as.data.frame(cv.auto.10,rows="cv",columns="criteria" )$adjusted.criterioncv.mse.loo<-as.data.frame(cv.auto.loo,rows="cv",columns="criteria" )$criterionplot(c(1,10),range(cv.mse.10, cv.mse.loo),type ="n",xlab ="Degree of polynomial, p",ylab ="Cross-Validated MSE")lines(1:10, cv.mse.10,lwd =2,lty =1,col =2,pch =16,type ="b")lines(1:10, cv.mse.loo,lwd =2,lty =2,col =3,pch =17,type ="b")legend("topright",inset =0.02,legend =c("10-Fold CV","LOO CV"),lwd =2,lty =2:1,col =3:2,pch =17:16)Cross-validated 10-fold and LOO MSE as a function of polynomial degree,\(p\)
Alternatively, we can use theplot() method for"cvModList" objects to compare the models, though withseparate graphs for 10-fold and LOO CV:
plot(cv.auto.10,main="Polynomial Regressions, 10-Fold CV",axis.args=list(labels=1:10),xlab="Degree of Polynomial, p")plot(cv.auto.loo,main="Polynomial Regressions, LOO CV",axis.args=list(labels=1:10),xlab="Degree of Polynomial, p")Cross-validated 10-fold and LOO MSE as a function of polynomial degree,\(p\)
In this example, 10-fold and LOO CV produce generally similarresults, and also results that are similar to those produced by theestimated error variance\(\widehat{\sigma}^2\) for each model,reported above (except for the highest-degree polynomials, where the CVresults more clearly suggest over-fitting).
Let’s try fitting the polynomial regressions using a loop with themodel formulampg ~ poly(horsepower, p) (fitting 3 ratherthan 10 polynomial regressions for brevity):
mods<-vector(3,mode="list")for (pin1:3){ mods[[p]]<-lm(mpg~poly(horsepower, p),data = Auto)}cv(models(mods),data = Auto,seed =2120,method ="naive")#> Warning in checkFormula(model, colnames(data)): the following variable is#> in the model formula but not in the data set:#> p#> expect errors or incorrect results#> Warning in checkFormula(model, colnames(data)): the following variable is#> in the model formula but not in the data set:#> p#> expect errors or incorrect results#> Warning in checkFormula(model, colnames(data)): the following variable is#> in the model formula but not in the data set:#> p#> expect errors or incorrect results#> Model model.1:#> cross-validation criterion = 19.396#>#> Model model.2:#> cross-validation criterion = 19.396#>#> Model model.3:#> cross-validation criterion = 19.396p#> [1] 3Notice that all the MSE values (thecross-validation criterion) are the same, because the"naive" method (as opposed to the default"Woodbury" method) usesupdate() to refit themodels with each fold omitted, and when the models are refit, the value3 is used forp for all models. Joshua PhilippEntrop brought this problem to our attention
Logistic regression for theMroz data
TheMroz data set from thecarDatapackage(associated with Fox & Weisberg,2019) has been used by several authors to illustrate binarylogistic regression; see, in particularFox &Weisberg (2019). The data were originally drawn from the U.S.Panel Study of Income Dynamics and pertain to married women. Here are afew cases in the data set:
data("Mroz",package ="carData")head(Mroz,3)#> lfp k5 k618 age wc hc lwg inc#> 1 yes 1 0 32 no no 1.2102 10.91#> 2 yes 0 2 30 no no 0.3285 19.50#> 3 yes 1 3 35 no no 1.5141 12.04tail(Mroz,3)#> lfp k5 k618 age wc hc lwg inc#> 751 no 0 0 43 no no 0.88814 9.952#> 752 no 0 0 60 no no 1.22497 24.984#> 753 no 0 3 39 no no 0.85321 28.363The response variable in the logistic regression islfp,labor-force participation, a factor coded"yes" or"no". The remaining variables are predictors:
k5, number of children 5 years old of younger in thewoman’s household;k618, number of children between 6 and 18 yearsold;age, in years;wc, wife’s college attendance,"yes"or"no";hc, husband’s college attendance;lwg, the woman’s log wage rate if she is employed, orherimputed wage rate, if she is not(avariable that Fox & Weisberg, 2019 show is problematicallydefined); andinc, family income, in $1000s, exclusive of wife’sincome.
We use theglm() function to fit a binary logisticregression to theMroz data:
m.mroz<-glm(lfp~ .,data = Mroz,family = binomial)summary(m.mroz)#>#> Call:#> glm(formula = lfp ~ ., family = binomial, data = Mroz)#>#> Coefficients:#> Estimate Std. Error z value Pr(>|z|)#> (Intercept) 3.18214 0.64438 4.94 7.9e-07 ***#> k5 -1.46291 0.19700 -7.43 1.1e-13 ***#> k618 -0.06457 0.06800 -0.95 0.34234#> age -0.06287 0.01278 -4.92 8.7e-07 ***#> wcyes 0.80727 0.22998 3.51 0.00045 ***#> hcyes 0.11173 0.20604 0.54 0.58762#> lwg 0.60469 0.15082 4.01 6.1e-05 ***#> inc -0.03445 0.00821 -4.20 2.7e-05 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#>#> (Dispersion parameter for binomial family taken to be 1)#>#> Null deviance: 1029.75 on 752 degrees of freedom#> Residual deviance: 905.27 on 745 degrees of freedom#> AIC: 921.3#>#> Number of Fisher Scoring iterations: 4BayesRule(ifelse(Mroz$lfp=="yes",1,0),fitted(m.mroz,type ="response"))#> [1] 0.30677#> attr(,"casewise loss")#> [1] "y != round(yhat)"In addition to the usually summary output for a GLM, we show theresult of applying theBayesRule() function from thecv package to predictions derived from the fittedmodel. Bayes rule, which predicts a “success” in a binary regressionmodel when the fitted probability of success [i.e.,\(\phi = \Pr(y = 1)\)] is\(\widehat{\phi} \ge .5\) and a “failure” if\(\widehat{\phi} \lt .5\).3 The firstargument toBayesRule() is the binary {0, 1} response, andthe second argument is the predicted probability of success.BayesRule() returns the proportion of predictions that arein error, as appropriate for a “cost” function.
The value returned byBayesRule() is associated with an“attribute” named"casewise loss" and set to"y != round(yhat)", signifying that the Bayes rule CVcriterion is computed as the mean of casewise values, here 0 if theprediction for a case matches the observed value and 1 if it does not(signifying a prediction error). Themse() function fornumeric responses is also calculated as a casewise average. Some othercriteria, such as the median absolute error, computed by themedAbsErr() function in thecv package,aren’t averages of casewise components. The distinction is importantbecause, to our knowledge, the statistical theory of cross-validation,for example, inDavison & Hinkley(1997),Bates, Hastie, & Tibshirani(2023), andArlot & Celisse(2010), is developed for CV criteria like MSE that are means ofcasewise components. As a consequence, we limit computation of biasadjustment and confidence intervals (see below) to criteria that arecasewise averages.
In this example, the fitted logistic regression incorrectly predicts31% of the responses; we expect this estimate to be optimistic giventhat the model is used to “predict” the data to which it is fit.
The"glm" method forcv() is largelysimilar to the"lm" method, although the default algorithm,selected explicitly bymethod="exact", refits the modelwith each fold removed (and is thus equivalent tomethod="naive" for"lm" models). Forgeneralized linear models,method="Woodbury" or (for LOOCV)method="hatvalues" provide approximate results (see thecomputational and technical vignette for details):
summary(cv(m.mroz,criterion = BayesRule,seed =248))#> R RNG seed set to 248#> 10-Fold Cross Validation#> method: exact#> criterion: BayesRule#> cross-validation criterion = 0.32404#> bias-adjusted cross-validation criterion = 0.31952#> 95% CI for bias-adjusted CV criterion = (0.28607, 0.35297)#> full-sample criterion = 0.30677summary(cv(m.mroz,criterion = BayesRule,seed =248,method ="Woodbury"))#> R RNG seed set to 248#> 10-Fold Cross Validation#> method: Woodbury#> criterion: BayesRule#> cross-validation criterion = 0.32404#> bias-adjusted cross-validation criterion = 0.31926#> 95% CI for bias-adjusted CV criterion = (0.28581, 0.35271)#> full-sample criterion = 0.30677To ensure that the two methods use the same 10 folds, we specify theseed for R’s random-number generator explicitly; here, and as is commonin our experience, the"exact" and"Woodbury"algorithms produce nearly identical results. The CV estimates ofprediction error are slightly higher than the estimate based on all ofthe cases.
The printed output includes a 95% confidence interval for thebias-adjusted Bayes rule CV criterion.Bates etal. (2023) show that these confidence intervals are unreliablefor models fit to small samples, and by defaultcv()computes them only when the sample size is 400 or larger and when the CVcriterion employed is an average of casewise components, as is the casefor Bayes rule. See the final section of the vignette for details of thecomputation of confidence intervals for bias-adjusted CV criteria.
Here are results of applying LOO CV to the Mroz model, using both theexact and the approximate methods:
summary(cv(m.mroz,k ="loo",criterion = BayesRule))#> n-Fold Cross Validation#> method: exact#> criterion: BayesRule#> cross-validation criterion = 0.32005#> bias-adjusted cross-validation criterion = 0.3183#> 95% CI for bias-adjusted CV criterion = (0.28496, 0.35164)#> full-sample criterion = 0.30677summary(cv(m.mroz,k ="loo",criterion = BayesRule,method ="Woodbury"))#> n-Fold Cross Validation#> method: Woodbury#> criterion: BayesRule#> cross-validation criterion = 0.32005#> bias-adjusted cross-validation criterion = 0.3183#> 95% CI for bias-adjusted CV criterion = (0.28496, 0.35164)#> full-sample criterion = 0.30677summary(cv(m.mroz,k ="loo",criterion = BayesRule,method ="hatvalues"))#> n-Fold Cross Validation#> method: hatvalues#> criterion: BayesRule#> cross-validation criterion = 0.32005To the number of decimal digits shown, the three methods produceidentical results for this example.