<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/cross-validation-for-predictive-analytics-using-r/"></div></div><div><a data-url="http://www.r-bloggers.com/cross-validation-for-predictive-analytics-using-r/" data-counturl="http://www.r-bloggers.com/cross-validation-for-predictive-analytics-using-r/" data-count="horizontal" data-text="Cross-Validation for Predictive Analytics Using R" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/"> MilanoR</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<div>
<h3>Introduction</h3>
<p>Since ancient times, humankind has always avidly sought a way to predict the future. One of the most widely known examples of this kind of activity in the past is the <a target="_blank" href="https://en.wikipedia.org/wiki/Pythia">Oracle of Delphi</a>, who dispensed previews of the future to her petitioners in the form of divine inspired prophecies<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fn1"><sup>1</sup></a>. In the modern days, the desire to know the future is still of interest to many of us, even if my feeling is that the increasing rapidity of technology innovations we observe everyday has somewhat lessened this instinct: things that few years ago seemed futuristic are now available to the great mass (e.g. the World Wide Web).</p>
<p>Among the many areas of the human being where predictions are highly needed there is <em>business decision making</em>. The tools for formulating predictions about quantities of interest are commonly known as <strong>predictive analytics</strong>, which is itself an essential part of <a target="_blank" href="http://courses.csail.mit.edu/18.337/2015/docs/50YearsDataScience.pdf">data science</a>. At the heart of any prediction there is always a model, which typically depends on some unknown <em>structural parameters</em> (e.g. the coefficients of a regression model) as well as on one or more <em>tuning parameters</em> (e.g. the number of basis functions in a smoothing spline or the degree of a polynomial). The former are commonly estimated using a sample of data, while the latter have to be chosen to guarantee that the model itself provides predictions which are accurate enough. Tuning parameters usually regulate the model complexity and hence are a key ingredient for any predictive task. In this blog entry we focus on the most common strategy for eliciting reasonable values for the tuning parameters, the cross-validation approach.</p>
</div>
<div>
<h3>The Bias-Variance Dilemma</h3>
<p>The reason why one should care about the choice of the tuning parameter values is because these are intimately linked with the accuracy of the predictions returned by the model. What an analyst typically wants is a model that is able to predict well samples that have not been used for estimating the structural parameters (the so called <em>training sample</em>). In other words, a predictive model is considered good when it is capable of predicting previously unseen samples with high accuracy. The accuracy of a model’s predictions is usually gauged using a <strong>loss function</strong>. Popular choices for the loss functions are the mean-squared error for continuous outcomes, or the 0-1 loss for a categorical outcome<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fn2"><sup>2</sup></a>.</p>
<p>At this point, it is important to distinguish between different prediction error concepts:</p>
<ul>
<li>the <strong>training error</strong>, which is the average loss over the <strong>training sample</strong>,</li>
<li>the <strong>test error</strong>, the prediction error over an independent <strong>test sample</strong>.</li>
</ul>
<p>The training error gets smaller as long as the predicted responses are close to the observed responses, and will get larger if for some of the observations, the predicted and observed responses differ substantially. The training error is calculated using the training sample used to fit the model. Clearly, we shouldn’t care too much about the model’s predictive accuracy on the training data. On the contrary, we would like to assess the model’s ability to predict observations never seen during estimation. The test error provides a measure of this ability. In general, one should select the model corresponding to the lowest test error.</p>
<p>The <code>R</code> code below implements these idea via simulated data. In particular, I simulate 100 training sets each of size 50 from a polynomial regression model, and for each I fit a sequence of cubic spline models with degrees of freedom from 1 to 30.</p>
<pre># Generate the training and test samples
seed &lt;- 1809
set.seed(seed)

gen_data &lt;- function(n, beta, sigma_eps) {
    eps &lt;- rnorm(n, 0, sigma_eps)
    x &lt;- sort(runif(n, 0, 100))
    X &lt;- cbind(1, poly(x, degree = (length(beta) - 1), raw = TRUE))
    y &lt;- as.numeric(X %*% beta + eps)
    
    return(data.frame(x = x, y = y))
}

# Fit the models
require(splines)

n_rep &lt;- 100
n_df &lt;- 30
df &lt;- 1:n_df
beta &lt;- c(5, -0.1, 0.004, -3e-05)
n_train &lt;- 50
n_test &lt;- 10000
sigma_eps &lt;- 0.5

xy &lt;- res &lt;- list()
xy_test &lt;- gen_data(n_test, beta, sigma_eps)
for (i in 1:n_rep) {
    xy[[i]] &lt;- gen_data(n_train, beta, sigma_eps)
    x &lt;- xy[[i]][, &quot;x&quot;]
    y &lt;- xy[[i]][, &quot;y&quot;]
    res[[i]] &lt;- apply(t(df), 2, function(degf) lm(y ~ ns(x, df = degf)))
}</pre>
<p>The next plot shows the first simulated training sample together with three fitted models corresponding to cubic splines with 1 (green line), 4 (orange line) and 25 (blue line) degrees of freedom respectively. These numbers have been chosen to show the full set of possibilities one may encounter in practice, i.e., either a model with low variability but high bias (degrees of freedom = 1), or a model with high variability but low bias (degrees of freedom = 25), or a model which tries to find a compromise between bias and variance (degrees of freedom = 4).</p>
<pre># Plot the data
x &lt;- xy[[1]]$x
X &lt;- cbind(1, poly(x, degree = (length(beta) - 1), raw = TRUE))
y &lt;- xy[[1]]$y
plot(y ~ x, col = &quot;gray&quot;, lwd = 2)
lines(x, X %*% beta, lwd = 3, col = &quot;black&quot;)
lines(x, fitted(res[[1]][[1]]), lwd = 3, col = &quot;palegreen3&quot;)
lines(x, fitted(res[[1]][[4]]), lwd = 3, col = &quot;darkorange&quot;)
lines(x, fitted(res[[1]][[25]]), lwd = 3, col = &quot;steelblue&quot;)
legend(x = &quot;topleft&quot;, legend = c(&quot;True function&quot;, &quot;Linear fit (df = 1)&quot;, &quot;Best model (df = 4)&quot;, 
    &quot;Overfitted model (df = 25)&quot;), lwd = rep(3, 4), col = c(&quot;black&quot;, &quot;palegreen3&quot;, 
    &quot;darkorange&quot;, &quot;steelblue&quot;), text.width = 32, cex = 0.85)</pre>
</div>
<p><a target="_blank" rel="attachment wp-att-2268" href="http://i2.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_1-e1462284035741.png"><img sizes="(max-width: 571px) 100vw, 571px" srcset="http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_1-e1462284035741-300x191.png 300w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_1-e1462284035741-768x488.png 768w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_1-e1462284035741-1024x651.png 1024w, http://i2.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_1-e1462284035741.png?resize=456%2C290 1310w" alt="cross-validation_1" class="wp-image-2268" src="http://i2.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_1-e1462284035741.png?resize=456%2C290" data-recalc-dims="1"></a></p>
<div>
<p>Then, for each training sample and fitted model, I compute the corresponding test error using a large test sample generated from the same (known!) population. These are represented in the following plot together with their <em>averages</em>, which are shown using thicker lines<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fn3"><sup>3</sup></a>. The solid points represent the three models illustrated in the previous diagram.</p>
<pre># Compute the training and test errors for each model
pred &lt;- list()
mse &lt;- te &lt;- matrix(NA, nrow = n_df, ncol = n_rep)
for (i in 1:n_rep) {
    mse[, i] &lt;- sapply(res[[i]], function(obj) deviance(obj)/nobs(obj))
    pred[[i]] &lt;- mapply(function(obj, degf) predict(obj, data.frame(x = xy_test$x)), 
        res[[i]], df)
    te[, i] &lt;- sapply(as.list(data.frame(pred[[i]])), function(y_hat) mean((xy_test$y - 
        y_hat)^2))
}

# Compute the average training and test errors
av_mse &lt;- rowMeans(mse)
av_te &lt;- rowMeans(te)

# Plot the errors
plot(df, av_mse, type = &quot;l&quot;, lwd = 2, col = gray(0.4), ylab = &quot;Prediction error&quot;, 
    xlab = &quot;Flexibilty (spline's degrees of freedom [log scaled])&quot;, ylim = c(0, 
        1), log = &quot;x&quot;)
abline(h = sigma_eps, lty = 2, lwd = 0.5)
for (i in 1:n_rep) {
    lines(df, te[, i], col = &quot;lightpink&quot;)
}
for (i in 1:n_rep) {
    lines(df, mse[, i], col = gray(0.8))
}
lines(df, av_mse, lwd = 2, col = gray(0.4))
lines(df, av_te, lwd = 2, col = &quot;darkred&quot;)
points(df[1], av_mse[1], col = &quot;palegreen3&quot;, pch = 17, cex = 1.5)
points(df[1], av_te[1], col = &quot;palegreen3&quot;, pch = 17, cex = 1.5)
points(df[which.min(av_te)], av_mse[which.min(av_te)], col = &quot;darkorange&quot;, pch = 16, 
    cex = 1.5)
points(df[which.min(av_te)], av_te[which.min(av_te)], col = &quot;darkorange&quot;, pch = 16, 
    cex = 1.5)
points(df[25], av_mse[25], col = &quot;steelblue&quot;, pch = 15, cex = 1.5)
points(df[25], av_te[25], col = &quot;steelblue&quot;, pch = 15, cex = 1.5)
legend(x = &quot;top&quot;, legend = c(&quot;Training error&quot;, &quot;Test error&quot;), lwd = rep(2, 2), 
    col = c(gray(0.4), &quot;darkred&quot;), text.width = 0.3, cex = 0.85)</pre>
<p><a target="_blank" rel="attachment wp-att-2269" href="http://i2.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_2-e1462283947527.png"><img sizes="(max-width: 561px) 100vw, 561px" srcset="http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_2-e1462283947527-300x194.png 300w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_2-e1462283947527-768x495.png 768w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_2-e1462283947527-1024x661.png 1024w, http://i2.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_2-e1462283947527.png?resize=456%2C294 1333w" alt="cross-validation_2" class="wp-image-2269" src="http://i2.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_2-e1462283947527.png?resize=456%2C294" data-recalc-dims="1"></a></p>
<p>One can see that the training errors decrease monotonically as the model gets more complicated (and less smooth). On the other side, even if the test error initially decreases, from a certain flexibility level on it starts increasing again. The change point occurs in correspondence of the orange model, that is, the model that provides a good compromise between bias and variance. The reason why the test error starts increasing for degrees of freedom larger than 3 or 4 is the so called <strong>overfitting</strong> problem. Overfitting is the tendency of a model to adapt too well to the training data, at the expense of generalization to previously unseen data points. In other words, an overfitted model fits the noise in the data rather than the actual underlying relationships among the variables. Overfitting usually occurs when a model is unnecessarily complex.</p>
<p>It is possible to show that the (expected) test error for a given observation in the test set can be decomposed into the sum of three components, namely</p>
<div><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;block&quot;&gt;&lt;mtext&gt;Expected Test Error&lt;/mtext&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mtext&gt;Irreducible Noise&lt;/mtext&gt;&lt;mo&gt;+&lt;/mo&gt;&lt;mo stretchy=&quot;false&quot;&gt;(&lt;/mo&gt;&lt;mtext&gt;Model Bias&lt;/mtext&gt;&lt;msup&gt;&lt;mo stretchy=&quot;false&quot;&gt;)&lt;/mo&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/msup&gt;&lt;mo&gt;+&lt;/mo&gt;&lt;mtext&gt;Model Variance&lt;/mtext&gt;&lt;mo&gt;,&lt;/mo&gt;&lt;/math&gt;" tabindex="0"><span><span><span>Expected Test Error </span><span>= </span><span>Irreducible Noise </span><span>+ </span><span>(</span><span>Model Bias</span><span><span>)^</span><span>2 </span></span><span>+ </span><span>Model Variance</span></span></span></span></div>
<div></div>
<div>which is known as the <strong>bias-variance decomposition</strong>. The first term is the data generating process variance. This term is unavoidable because we live in a noisy stochastic world, where even the best ideal model has non-zero error. The second term originates from the difficulty to catch the correct functional form of the relationship that links the dependent and independent variables (sometimes it is also called the <em>approximation bias</em>). The last term is due to the fact that we estimate our models using only a limited amount of data. Fortunately, this terms gets closer and closer to zero as long as we collect more and more training data. Typically, the more complex (i.e., flexible) we make the model, the lower the bias but the higher the variance. This general phenomenon is known as the<strong>bias-variance trade-off</strong>, and the challenge is to find a model which provides a good compromise between these two issues.</div>
<p>Clearly, the situation illustrated above is only ideal, because in practice:</p>
<ul>
<li>We do not know the true model that generates the data. Indeed, our models are typically more or less mis-specified.</li>
<li>We do only have a limited amount of data.</li>
</ul>
<p>One way to overcome these hurdles and approximate the search for the optimal model is to use the cross-validation approach.</p>
</div>
<div>
<h3>A Solution: Cross-Validation</h3>
<p>In essence, all these ideas bring us to the conclusion that it is not advisable to compare the predictive accuracy of a set of models using the same observations used for estimating the models. Therefore, for assessing the models’ predictive performance we should use an independent set of data (the test sample). Then, the model showing the lowest error on the test sample (i.e., the lowest test error) is identified as the best.</p>
<p>Unfortunately, in many cases it is not possible to draw a (possibly large) independent set of observations for testing the models’ performance, because collecting data is typically an expensive activity. The immediate reaction to this statement is that we can solve this issue by <em>splitting</em> the available data in two sets, one of which will be used for training while the other is used for testing. The split is usually performed randomly to guarantee that the two parts have the same distribution<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fn4"><sup>4</sup></a>.</p>
<p>Even if data splitting provides an unbiased estimate of the test error, it is often quite noisy. A possible solution<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fn5"><sup>5</sup></a> is to use <strong>cross-validation</strong> (CV). In its basic version, the so called <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span></span></span><span>k</span></span></span>-fold cross-validation, the samples are randomly partitioned into <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span></span></span><span>k</span></span></span> sets (called <em>folds</em>) of roughly equal size. A model is fit using all the samples except the first subset. Then, the prediction error of the fitted model is calculated using the first held-out samples. The same operation is repeated for each fold and the model’s performance is calculated by averaging the errors across the different test sets. <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span></span></span><span>k</span></span></span> is usually fixed at 5 or 10 . Cross-validation provides an estimate of the test error for each model<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fn6"><sup>6</sup></a>. Cross-validation is one of the most widely-used method for model selection, and for choosing tuning parameter values.</p>
<p>The code below illustrates <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span></span></span><span>k</span></span></span>-fold cross-validation using the same simulated data as above but not pretending to know the data generating process. In particular, I generate 100 observations and choose <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mn&gt;10&lt;/mn&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span><span>=</span><span>10</span></span></span><span>k=10</span></span></span>. Together with the training error curve, in the plot I report both the CV and test error curves. Additionally, I provide also the standard error bars, which are the standard errors of the individual prediction error for each of the <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mn&gt;10&lt;/mn&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span><span>=</span><span>10</span></span></span><span>k=10</span></span></span> parts.</p>
<pre>set.seed(seed)

n_train &lt;- 100
xy &lt;- gen_data(n_train, beta, sigma_eps)
x &lt;- xy$x
y &lt;- xy$y

fitted_models &lt;- apply(t(df), 2, function(degf) lm(y ~ ns(x, df = degf)))
mse &lt;- sapply(fitted_models, function(obj) deviance(obj)/nobs(obj))

n_test &lt;- 10000
xy_test &lt;- gen_data(n_test, beta, sigma_eps)
pred &lt;- mapply(function(obj, degf) predict(obj, data.frame(x = xy_test$x)), 
    fitted_models, df)
te &lt;- sapply(as.list(data.frame(pred)), function(y_hat) mean((xy_test$y - y_hat)^2))

n_folds &lt;- 10
folds_i &lt;- sample(rep(1:n_folds, length.out = n_train))
cv_tmp &lt;- matrix(NA, nrow = n_folds, ncol = length(df))
for (k in 1:n_folds) {
    test_i &lt;- which(folds_i == k)
    train_xy &lt;- xy[-test_i, ]
    test_xy &lt;- xy[test_i, ]
    x &lt;- train_xy$x
    y &lt;- train_xy$y
    fitted_models &lt;- apply(t(df), 2, function(degf) lm(y ~ ns(x, df = degf)))
    x &lt;- test_xy$x
    y &lt;- test_xy$y
    pred &lt;- mapply(function(obj, degf) predict(obj, data.frame(ns(x, df = degf))), 
        fitted_models, df)
    cv_tmp[k, ] &lt;- sapply(as.list(data.frame(pred)), function(y_hat) mean((y - 
        y_hat)^2))
}
cv &lt;- colMeans(cv_tmp)

require(Hmisc)

plot(df, mse, type = &quot;l&quot;, lwd = 2, col = gray(0.4), ylab = &quot;Prediction error&quot;, 
    xlab = &quot;Flexibilty (spline's degrees of freedom [log scaled])&quot;, main = paste0(n_folds, 
        &quot;-fold Cross-Validation&quot;), ylim = c(0.1, 0.8), log = &quot;x&quot;)
lines(df, te, lwd = 2, col = &quot;darkred&quot;, lty = 2)
cv_sd &lt;- apply(cv_tmp, 2, sd)/sqrt(n_folds)
errbar(df, cv, cv + cv_sd, cv - cv_sd, add = TRUE, col = &quot;steelblue2&quot;, pch = 19, 
    lwd = 0.5)
lines(df, cv, lwd = 2, col = &quot;steelblue2&quot;)
points(df, cv, col = &quot;steelblue2&quot;, pch = 19)
legend(x = &quot;topright&quot;, legend = c(&quot;Training error&quot;, &quot;Test error&quot;, &quot;Cross-validation error&quot;), 
    lty = c(1, 2, 1), lwd = rep(2, 3), col = c(gray(0.4), &quot;darkred&quot;, &quot;steelblue2&quot;), 
    text.width = 0.4, cex = 0.85)</pre>
<p><a target="_blank" rel="attachment wp-att-2270" href="http://i0.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_3.png"><img sizes="(max-width: 562px) 100vw, 562px" srcset="http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_3-300x214.png 300w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_3-768x549.png 768w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_3-1024x731.png 1024w, http://i0.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_3.png?resize=456%2C325 1344w" alt="cross-validation_3" class="wp-image-2270" src="http://i0.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_3.png?resize=456%2C325" data-recalc-dims="1"></a></p>
<p>Often a “one-standard error” rule is used with cross-validation, according to which one should choose the most parsimonious model whose error is no more than one standard error above the error of the best model. In the example above, the best model (that for which the CV error is minimized) uses 3 degrees of freedom, which also satisfies the requirement of the one-standard error rule.</p>
<p>The case where <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mi&gt;n&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span><span>=</span><span>n</span></span></span><span>k=n</span></span></span> corresponds to the so called <strong>leave-one-out cross-validation</strong> (LOOCV) method. In this case the test set contains a single observation. The advantages of LOOCV are: 1) it doesn’t require random numbers to select the observations to test, meaning that it doesn’t produce different results when applied repeatedly, and 2) it has far less bias than <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span></span></span><span>k</span></span></span>-fold CV because it employs larger training sets containing <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;n&lt;/mi&gt;&lt;mo&gt;&amp;#x2212;&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/math&gt;" tabindex="0"><span><span><span>n</span><span>−</span><span>1</span></span></span><span>n−1</span></span></span> observations each. On the other side, LOOCV presents also some drawbacks: 1) it is potentially quite intense computationally, and 2) due to the fact that any two training sets share <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;n&lt;/mi&gt;&lt;mo&gt;&amp;#x2212;&lt;/mo&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/math&gt;" tabindex="0"><span><span><span>n</span><span>−</span><span>2</span></span></span><span>n−2</span></span></span> points, the models fit to those training sets tend to be strongly correlated with each other.</p>
<p>The code below implements LOOCV using the same example I discussed so far. The next plot shows that most of the times LOOCV does not provide dramatically different results with respect to CV.</p>
<pre>require(splines)

loocv_tmp &lt;- matrix(NA, nrow = n_train, ncol = length(df))
for (k in 1:n_train) {
  train_xy &lt;- xy[-k, ]
  test_xy &lt;- xy[k, ]
  x &lt;- train_xy$x
  y &lt;- train_xy$y
  fitted_models &lt;- apply(t(df), 2, function(degf) lm(y ~ ns(x, df = degf)))
  pred &lt;- mapply(function(obj, degf) predict(obj, data.frame(x = test_xy$x)),
                 fitted_models, df)
  loocv_tmp[k, ] &lt;- (test_xy$y - pred)^2
}
loocv &lt;- colMeans(loocv_tmp)

plot(df, mse, type = &quot;l&quot;, lwd = 2, col = gray(.4), ylab = &quot;Prediction error&quot;,
     xlab = &quot;Flexibilty (spline's degrees of freedom [log scaled])&quot;,
     main = &quot;Leave-One-Out Cross-Validation&quot;, ylim = c(.1, .8), log = &quot;x&quot;)
lines(df, cv, lwd = 2, col = &quot;steelblue2&quot;, lty = 2)
lines(df, loocv, lwd = 2, col = &quot;darkorange&quot;)
legend(x = &quot;topright&quot;, legend = c(&quot;Training error&quot;, &quot;10-fold CV error&quot;, &quot;LOOCV error&quot;),
       lty = c(1, 2, 1), lwd = rep(2, 3), col = c(gray(.4), &quot;steelblue2&quot;, &quot;darkorange&quot;),
       text.width = .3, cex = .85)</pre>

</div>
<div>
<h2><a target="_blank" rel="attachment wp-att-2271" href="http://i1.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_4.png"><img sizes="(max-width: 600px) 100vw, 600px" srcset="http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_4-300x214.png 300w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_4-768x549.png 768w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_4-1024x731.png 1024w, http://i1.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_4.png?resize=456%2C325 1344w" alt="cross-validation_4" class="wp-image-2271" src="http://i1.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_4.png?resize=456%2C325" data-recalc-dims="1"></a></h2>
<h3>Doing Cross-Validation With <code>R</code>: the <code>caret</code> Package</h3>
<p>There are many <code>R</code> packages that provide functions for performing different flavors of CV. In my opinion, one of the best implementation of these ideas is available in the <code>caret</code> package by Max Kuhn <span>(see Kuhn and Johnson 2013)</span><a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fn7"><sup>7</sup></a>. The aim of the <code>caret</code> package (acronym of <strong>c</strong>lassification <strong>a</strong>nd <strong>re</strong>gression <strong>t</strong>raining) is to provide a very general and efficient suite of commands for building and assessing predictive models. It allows to compare the predictive accuracy of a multitude of models (currently more than 200), including the most recent ones from machine learning. The comparison of different models can be done using cross-validation as well as with other approaches. The package also provides many options for data pre-processing. It is not my aim to provide here a thorough presentation of all the package features. Rather, I will focus only on a handful of its functions, those that allow to perform CV. For more details on the other package functions, you can inspect the package documentation and its <a target="_blank" href="http://topepo.github.io/caret/">website</a>. To illustrate these feature I will use some data for a credit scoring application whose data can be found <a target="_blank" href="https://github.com/gastonstat/CreditScoring">here</a>.</p>
<p>Since credit scoring is a classification problem, I will use the number of misclassified observations as the loss measure. The data set contains information about 4,455 individuals for the following variables:</p>
<table>
<thead>
<tr>
<th align="left">Variable</th>
<th align="left">Description</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left"><em>Status</em></td>
<td align="left">credit status</td>
</tr>
<tr>
<td align="left"><em>Seniority</em></td>
<td align="left">job seniority (years)</td>
</tr>
<tr>
<td align="left"><em>Home</em></td>
<td align="left">type of home ownership</td>
</tr>
<tr>
<td align="left"><em>Time</em></td>
<td align="left">time of requested loan</td>
</tr>
<tr>
<td align="left"><em>Age</em></td>
<td align="left">client’s age</td>
</tr>
<tr>
<td align="left"><em>Marital</em></td>
<td align="left">marital status</td>
</tr>
<tr>
<td align="left"><em>Records</em></td>
<td align="left">existence of records</td>
</tr>
<tr>
<td align="left"><em>Job</em></td>
<td align="left">type of job</td>
</tr>
<tr>
<td align="left"><em>Expenses</em></td>
<td align="left">amount of expenses</td>
</tr>
<tr>
<td align="left"><em>Income</em></td>
<td align="left">amount of income</td>
</tr>
<tr>
<td align="left"><em>Assets</em></td>
<td align="left">amount of assets</td>
</tr>
<tr>
<td align="left"><em>Debt</em></td>
<td align="left">amount of debt</td>
</tr>
<tr>
<td align="left"><em>Amount</em></td>
<td align="left">amount requested of loan</td>
</tr>
<tr>
<td align="left"><em>Price</em></td>
<td align="left">price of good</td>
</tr>
</tbody>
</table>
<p>Here I use the “cleaned” version of the data set, where some pre-processing has already been performed (i.e., removal of few observations, imputation of missing values and categorization of continuous predictors). The tidy data are contained in the file <code>CleanCreditScoring.csv</code>.</p>
<pre>require(RCurl)
require(prettyR)

url &lt;- &quot;https://raw.githubusercontent.com/gastonstat/CreditScoring/master/CleanCreditScoring.csv&quot;
cs_data &lt;- getURL(url)
cs_data &lt;- read.csv(textConnection(cs_data))
describe(cs_data)</pre>

<pre>## Description of cs_data</pre>

<pre>## 
##  Numeric 
##              mean  median          var       sd valid.n
## Seniority    7.99    5.00        66.85     8.18    4446
## Time        46.45   48.00       214.56    14.65    4446
## Age         37.08   36.00       120.70    10.99    4446
## Expenses    55.60   51.00       381.06    19.52    4446
## Income     140.63  124.00      6428.50    80.18    4446
## Assets    5354.95 3000.00 133040726.62 11534.33    4446
## Debt       342.26    0.00   1549264.52  1244.69    4446
## Amount    1038.76 1000.00    225385.62   474.75    4446
## Price     1462.48 1400.00    395081.60   628.56    4446
## Finrat      72.62   77.10       415.78    20.39    4446
## Savings      3.86    3.12        13.89     3.73    4446
## 
##  Factor 
##          
## Status       good     bad
##   Count   3197.00 1249.00
##   Percent   71.91   28.09
## Mode good 
##          
## Home        owner   rent parents  other   priv ignore
##   Count   2106.00 973.00  782.00 319.00 246.00  20.00
##   Percent   47.37  21.88   17.59   7.17   5.53   0.45
## Mode owner 
##          
## Marital   married single separated widow divorced
##   Count   3238.00 973.00    130.00 67.00    38.00
##   Percent   72.83  21.88      2.92  1.51     0.85
## Mode married 
##          
## Records   no_rec yes_rec
##   Count   3677.0   769.0
##   Percent   82.7    17.3
## Mode no_rec 
##          
## Job         fixed freelance partime others
##   Count   2803.00   1021.00  451.00 171.00
##   Percent   63.05     22.96   10.14   3.85
## Mode fixed 
##           
## seniorityR sen (-1,1] sen (3,8] sen (14,99] sen (1,3] sen (8,14]
##    Count      1042.00       978      880.00    789.00     757.00
##    Percent      23.44        22       19.79     17.75      17.03
## Mode sen (-1,1] 
##          
## timeR     time (48,99] time (24,36] time (36,48] time (12,24] time (0,12]
##   Count        1949.00       991.00       885.00       441.00      180.00
##   Percent        43.84        22.29        19.91         9.92        4.05
## Mode time (48,99] 
##          
## ageR      age (30,40] age (40,50] age (25,30] age (0,25] age (50,99]
##   Count       1415.00      900.00      781.00     699.00      651.00
##   Percent       31.83       20.24       17.57      15.72       14.64
## Mode age (30,40] 
##          
## expensesR exp (0,40] exp (40,50] exp (50,60] exp (60,80] exp (80,1e+04]
##   Count      1219.00      999.00      979.00      798.00         451.00
##   Percent      27.42       22.47       22.02       17.95          10.14
## Mode exp (0,40] 
##          
## incomeR   inc (80,110] inc (140,190] inc (0,80] inc (110,140]
##   Count         954.00        915.00     886.00        866.00
##   Percent        21.46         20.58      19.93         19.48
##          
## incomeR   inc (190,1e+04]
##   Count            825.00
##   Percent           18.56
## Mode inc (80,110] 
##          
## assetsR   asset (-1,0] asset (3e+03,5e+03] asset (8e+03,1e+06]
##   Count        1626.00              937.00              719.00
##   Percent        36.57               21.08               16.17
##          
## assetsR   asset (0,3e+03] asset (5e+03,8e+03]
##   Count            626.00               538.0
##   Percent           14.08                12.1
## Mode asset (-1,0] 
##          
## debtR     debt (-1,0] debt (500,1.5e+03] debt (2.5e+03,1e+06] debt (0,500]
##   Count       3667.00             230.00               197.00       193.00
##   Percent       82.48               5.17                 4.43         4.34
##          
## debtR     debt (1.5e+03,2.5e+03]
##   Count                   159.00
##   Percent                   3.58
## Mode debt (-1,0] 
##          
## amountR   am (900,1.1e+03] am (1.1e+03,1.4e+03] am (600,900] am (0,600]
##   Count             945.00               925.00       911.00     895.00
##   Percent            21.26                20.81        20.49      20.13
##          
## amountR   am (1.4e+03,1e+05]
##   Count               770.00
##   Percent              17.32
## Mode am (900,1.1e+03] 
##          
## priceR    priz (1.5e+03,1.8e+03] priz (1e+03,1.3e+03] priz (0,1e+03]
##   Count                  1028.00               985.00         821.00
##   Percent                  23.12                22.15          18.47
##          
## priceR    priz (1.8e+03,1e+05] priz (1.3e+03,1.5e+03]
##   Count                 811.00                 801.00
##   Percent                18.24                  18.02
## Mode priz (1.5e+03,1.8e+03] 
##          
## finratR   finr (80,90] finr (90,100] finr (50,70] finr (70,80] finr (0,50]
##   Count         995.00        960.00       954.00       821.00       716.0
##   Percent        22.38         21.59        21.46        18.47        16.1
## Mode finr (80,90] 
##          
## savingsR  sav (2,4] sav (0,2] sav (6,99] sav (4,6] sav (-99,0]
##   Count      1396.0   1111.00      827.0    814.00       298.0
##   Percent      31.4     24.99       18.6     18.31         6.7
## Mode sav (2,4]</pre>
<p>The <code>caret</code> package provides functions for splitting the data as well as functions that automatically do all the job for us, namely functions that create the resampled data sets, fit the models, and evaluate performance.</p>
<p>Among the functions for data splitting I just mention <code>createDataPartition()</code> and <code>createFolds()</code>. The former allows to create one or more test/training random partitions of the data, while the latter randomly splits the data into <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span></span></span><span>k</span></span></span> subsets. In both functions the random sampling is done within the levels of <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;y&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>y</span></span></span><span>y</span></span></span> (when <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;y&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>y</span></span></span><span>y</span></span></span> is categorical) to balance the class distributions within the splits. These functions return vectors of indexes that can then be used to subset the original sample into training and test sets.</p>
<pre>require(caret)

classes &lt;- cs_data[, &quot;Status&quot;]
predictors &lt;- cs_data[, -match(c(&quot;Status&quot;, &quot;Seniority&quot;, &quot;Time&quot;, &quot;Age&quot;, &quot;Expenses&quot;, 
    &quot;Income&quot;, &quot;Assets&quot;, &quot;Debt&quot;, &quot;Amount&quot;, &quot;Price&quot;, &quot;Finrat&quot;, &quot;Savings&quot;), colnames(cs_data))]

train_set &lt;- createDataPartition(classes, p = 0.8, list = FALSE)
str(train_set)</pre>

<pre>##  int [1:3558, 1] 1 2 3 4 5 6 7 8 9 11 ...
##  - attr(*, &quot;dimnames&quot;)=List of 2
##   ..$ : NULL
##   ..$ : chr &quot;Resample1&quot;</pre>

<pre>train_predictors &lt;- predictors[train_set, ]
train_classes &lt;- classes[train_set]
test_predictors &lt;- predictors[-train_set, ]
test_classes &lt;- classes[-train_set]

set.seed(seed)
cv_splits &lt;- createFolds(classes, k = 10, returnTrain = TRUE)
str(cv_splits)</pre>

<pre>## List of 10
##  $ Fold01: int [1:4002] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold02: int [1:4002] 2 3 4 5 6 7 8 9 10 11 ...
##  $ Fold03: int [1:4001] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold04: int [1:4002] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold05: int [1:4001] 1 2 3 4 7 8 10 11 12 13 ...
##  $ Fold06: int [1:4001] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold07: int [1:4001] 1 2 4 5 6 7 9 10 11 12 ...
##  $ Fold08: int [1:4002] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold09: int [1:4001] 1 2 3 5 6 8 9 11 12 13 ...
##  $ Fold10: int [1:4001] 1 3 4 5 6 7 8 9 10 11 ...</pre>
<p>To automatically split the data, fit the models and assess the performance, one can use the <code>train()</code> function in the <code>caret</code> package. The code below shows an example of the <code>train()</code> function on the credit scoring data by modeling the outcome using all the predictors available with a penalized logistic regression. More specifically, I use the <code>glmnet</code> package <span>(Friedman, Hastie, and Tibshirani 2008)</span>, that fits a generalized linear model via penalized maximum likelihood. The algorithm implemented in the package computes the regularization path for the <em>elastic-net</em> penalty over a grid of values for the <em>regularization parameter</em> <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03BB;&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>λ</span></span></span><span>λ</span></span></span>. The tuning parameter <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03BB;&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>λ</span></span></span><span>λ</span></span></span> controls the overall strength of the penalty. A second tuning parameter, called the <em>mixing percentage</em> and denoted with <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03B1;&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>α</span></span></span><span>α</span></span></span>, represents the elastic-net penalty <span>(Zou and Hastie 2005)</span>. This parameter takes value in <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mo stretchy=&quot;false&quot;&gt;[&lt;/mo&gt;&lt;mn&gt;0&lt;/mn&gt;&lt;mo&gt;,&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;mo stretchy=&quot;false&quot;&gt;]&lt;/mo&gt;&lt;/math&gt;" tabindex="0"><span><span><span>[</span><span>0</span><span>,</span><span>1</span><span>]</span></span></span><span>[0,1]</span></span></span> and bridges the gap between the <em>lasso</em> (<span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03B1;&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/math&gt;" tabindex="0"><span><span><span>α</span><span>=</span><span>1</span></span></span><span>α=1</span></span></span>) and the <em>ridge</em> (<span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03B1;&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mn&gt;0&lt;/mn&gt;&lt;/math&gt;" tabindex="0"><span><span><span>α</span><span>=</span><span>0</span></span></span><span>α=0</span></span></span>) approaches.</p>
<p>The <code>train()</code> function requires the model formula together with the indication of the model to fit and the grid of tuning parameter values to use. In the code below this grid is specified through the <code>tuneGrid</code> argument, while <code>trControl</code> provides the method to use for choosing the optimal values of the tuning parameters (in our case, 10-fold cross-validation). Finally, the <code>preProcess</code> argument allows to apply a series of pre-processing operations on the predictors (in our case, centering and scaling the predictor values).</p>
<pre>require(glmnet)</pre>

<pre>## Warning: package 'Matrix' was built under R version 3.2.5</pre>

<pre>set.seed(seed)

cs_data_train &lt;- cs_data[train_set, ]
cs_data_test &lt;- cs_data[-train_set, ]

glmnet_grid &lt;- expand.grid(alpha = c(0,  .1,  .2, .4, .6, .8, 1),
                           lambda = seq(.01, .2, length = 20))
glmnet_ctrl &lt;- trainControl(method = &quot;cv&quot;, number = 10)
glmnet_fit &lt;- train(Status ~ ., data = cs_data_train,
                    method = &quot;glmnet&quot;,
                    preProcess = c(&quot;center&quot;, &quot;scale&quot;),
                    tuneGrid = glmnet_grid,
                    trControl = glmnet_ctrl)
glmnet_fit</pre>

<pre>## glmnet 
## 
## 3558 samples
##   26 predictor
##    2 classes: 'bad', 'good' 
## 
## Pre-processing: centered (68), scaled (68) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 3202, 3203, 3202, 3203, 3202, 3202, ... 
## Resampling results across tuning parameters:
## 
##   alpha  lambda  Accuracy   Kappa       
##   0.0    0.01    0.8021427  0.4613907413
##   0.0    0.02    0.7998916  0.4520486081
##   0.0    0.03    0.7976412  0.4402614685
##   0.0    0.04    0.7987633  0.4407093800
##   0.0    0.05    0.7982015  0.4355350784
##   0.0    0.06    0.7979182  0.4313111542
##   0.0    0.07    0.7953893  0.4205306747
##   0.0    0.08    0.7931413  0.4105376360
##   0.0    0.09    0.7922978  0.4050557210
##   0.0    0.10    0.7892072  0.3920192662
##   0.0    0.11    0.7841454  0.3722554927
##   0.0    0.12    0.7824600  0.3640031420
##   0.0    0.13    0.7807739  0.3557226473
##   0.0    0.14    0.7793694  0.3482341774
##   0.0    0.15    0.7807746  0.3491274474
##   0.0    0.16    0.7810571  0.3472621824
##   0.0    0.17    0.7796511  0.3411817028
##   0.0    0.18    0.7796526  0.3373484610
##   0.0    0.19    0.7807746  0.3374411775
##   0.0    0.20    0.7785267  0.3288713594
##   0.1    0.01    0.8015794  0.4596697325
##   0.1    0.02    0.8010160  0.4506840219
##   0.1    0.03    0.7987672  0.4392696923
##   0.1    0.04    0.7962367  0.4270137579
##   0.1    0.05    0.7953909  0.4179208158
##   0.1    0.06    0.7922986  0.4039787586
##   0.1    0.07    0.7892056  0.3879423318
##   0.1    0.08    0.7880804  0.3782808312
##   0.1    0.09    0.7841454  0.3596294655
##   0.1    0.10    0.7807723  0.3432057135
##   0.1    0.11    0.7748718  0.3176779656
##   0.1    0.12    0.7726223  0.3058039618
##   0.1    0.13    0.7706544  0.2940236563
##   0.1    0.14    0.7672804  0.2743480275
##   0.1    0.15    0.7641905  0.2598076390
##   0.1    0.16    0.7613815  0.2443603691
##   0.1    0.17    0.7585718  0.2288614463
##   0.1    0.18    0.7549153  0.2107358484
##   0.1    0.19    0.7518231  0.1940038500
##   0.1    0.20    0.7495751  0.1809109004
##   0.2    0.01    0.8021396  0.4602812440
##   0.2    0.02    0.7982070  0.4415693991
##   0.2    0.03    0.7951147  0.4239150748
##   0.2    0.04    0.7917424  0.4066236599
##   0.2    0.05    0.7911742  0.3954568138
##   0.2    0.06    0.7875218  0.3766126319
##   0.2    0.07    0.7813388  0.3485492065
##   0.2    0.08    0.7765588  0.3236774842
##   0.2    0.09    0.7723398  0.3024505090
##   0.2    0.10    0.7706583  0.2883765872
##   0.2    0.11    0.7664425  0.2668721852
##   0.2    0.12    0.7627868  0.2474032680
##   0.2    0.13    0.7554779  0.2144108884
##   0.2    0.14    0.7509796  0.1882593847
##   0.2    0.15    0.7445165  0.1570413702
##   0.2    0.16    0.7419869  0.1374198064
##   0.2    0.17    0.7372108  0.1088243186
##   0.2    0.18    0.7346811  0.0884942046
##   0.2    0.19    0.7304645  0.0666304068
##   0.2    0.20    0.7293401  0.0587304848
##   0.4    0.01    0.7998916  0.4500521540
##   0.4    0.02    0.7951147  0.4262079087
##   0.4    0.03    0.7894920  0.3986179726
##   0.4    0.04    0.7833091  0.3669902064
##   0.4    0.05    0.7793694  0.3406658720
##   0.4    0.06    0.7765596  0.3202875649
##   0.4    0.07    0.7726246  0.2966774741
##   0.4    0.08    0.7639128  0.2564191601
##   0.4    0.09    0.7546352  0.2079469868
##   0.4    0.10    0.7456401  0.1554338724
##   0.4    0.11    0.7369299  0.1037171598
##   0.4    0.12    0.7324316  0.0745712422
##   0.4    0.13    0.7296218  0.0585154515
##   0.4    0.14    0.7282173  0.0515741506
##   0.4    0.15    0.7270937  0.0460444276
##   0.4    0.16    0.7237213  0.0275193624
##   0.4    0.17    0.7203474  0.0079998432
##   0.4    0.18    0.7189429  0.0000000000
##   0.4    0.19    0.7189429  0.0000000000
##   0.4    0.20    0.7189429  0.0000000000
##   0.6    0.01    0.7998940  0.4479646834
##   0.6    0.02    0.7900538  0.4056859000
##   0.6    0.03    0.7830274  0.3693670539
##   0.6    0.04    0.7765572  0.3347009487
##   0.6    0.05    0.7757161  0.3177186986
##   0.6    0.06    0.7655990  0.2676239096
##   0.6    0.07    0.7549153  0.2082278135
##   0.6    0.08    0.7422693  0.1330653578
##   0.6    0.09    0.7332758  0.0786721840
##   0.6    0.10    0.7296218  0.0585273347
##   0.6    0.11    0.7279364  0.0501825990
##   0.6    0.12    0.7223176  0.0195861846
##   0.6    0.13    0.7189429  0.0000000000
##   0.6    0.14    0.7189429  0.0000000000
##   0.6    0.15    0.7189429  0.0000000000
##   0.6    0.16    0.7189429  0.0000000000
##   0.6    0.17    0.7189429  0.0000000000
##   0.6    0.18    0.7189429  0.0000000000
##   0.6    0.19    0.7189429  0.0000000000
##   0.6    0.20    0.7189429  0.0000000000
##   0.8    0.01    0.7959582  0.4342802453
##   0.8    0.02    0.7869623  0.3901784671
##   0.8    0.03    0.7776832  0.3436570802
##   0.8    0.04    0.7745925  0.3177549651
##   0.8    0.05    0.7650348  0.2651758625
##   0.8    0.06    0.7518270  0.1887158977
##   0.8    0.07    0.7363681  0.0981410803
##   0.8    0.08    0.7299027  0.0607058454
##   0.8    0.09    0.7270937  0.0452185157
##   0.8    0.10    0.7189429  0.0008725808
##   0.8    0.11    0.7189429  0.0000000000
##   0.8    0.12    0.7189429  0.0000000000
##   0.8    0.13    0.7189429  0.0000000000
##   0.8    0.14    0.7189429  0.0000000000
##   0.8    0.15    0.7189429  0.0000000000
##   0.8    0.16    0.7189429  0.0000000000
##   0.8    0.17    0.7189429  0.0000000000
##   0.8    0.18    0.7189429  0.0000000000
##   0.8    0.19    0.7189429  0.0000000000
##   0.8    0.20    0.7189429  0.0000000000
##   1.0    0.01    0.7920241  0.4209892004
##   1.0    0.02    0.7807794  0.3655088353
##   1.0    0.03    0.7754360  0.3289251638
##   1.0    0.04    0.7686897  0.2864096482
##   1.0    0.05    0.7526713  0.1943914570
##   1.0    0.06    0.7349620  0.0898775655
##   1.0    0.07    0.7290600  0.0557404582
##   1.0    0.08    0.7209100  0.0108497047
##   1.0    0.09    0.7189429  0.0000000000
##   1.0    0.10    0.7189429  0.0000000000
##   1.0    0.11    0.7189429  0.0000000000
##   1.0    0.12    0.7189429  0.0000000000
##   1.0    0.13    0.7189429  0.0000000000
##   1.0    0.14    0.7189429  0.0000000000
##   1.0    0.15    0.7189429  0.0000000000
##   1.0    0.16    0.7189429  0.0000000000
##   1.0    0.17    0.7189429  0.0000000000
##   1.0    0.18    0.7189429  0.0000000000
##   1.0    0.19    0.7189429  0.0000000000
##   1.0    0.20    0.7189429  0.0000000000
## 
## Accuracy was used to select the optimal model using  the largest value.
## The final values used for the model were alpha = 0 and lambda = 0.01.</pre>

<pre>trellis.par.set(caretTheme())
plot(glmnet_fit, scales = list(x = list(log = 2)))</pre>
<p><a target="_blank" rel="attachment wp-att-2272" href="http://i0.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_5.png"><img sizes="(max-width: 565px) 100vw, 565px" srcset="http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_5-300x214.png 300w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_5-768x549.png 768w, http://www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_5-1024x731.png 1024w, http://i0.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_5.png?resize=456%2C325 1344w" alt="cross-validation_5" class="wp-image-2272" src="http://i0.wp.com/www.milanor.net/blog/wp-content/uploads/2016/05/articolo_milanor_5.png?resize=456%2C325" data-recalc-dims="1"></a></p>
<p>The previous plot shows the “accuracy”, that is the percentage of correctly classified observations, for the penalized logistic regression model with each combination of the two tuning parameters <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03B1;&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>α</span></span></span><span>α</span></span></span>and <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03BB;&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>λ</span></span></span><span>λ</span></span></span>. The optimal tuning parameter values are <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03B1;&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;/math&gt;" tabindex="0"><span><span><span>α</span><span>=</span></span></span><span>α=</span></span></span> 0 and <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;&amp;#x03BB;&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;/math&gt;" tabindex="0"><span><span><span>λ</span><span>=</span></span></span><span>λ=</span></span></span> 0.01.</p>
<p>Then, it is possible to predict new samples with the identified optimal model using the <code>predict</code> method:</p>
<pre>pred_classes &lt;- predict(glmnet_fit, newdata = cs_data_test)
table(pred_classes)</pre>

<pre>## pred_classes
##  bad good 
##  172  716</pre>

<pre>pred_probs &lt;- predict(glmnet_fit, newdata = cs_data_test, type = &quot;prob&quot;)
head(pred_probs)</pre>

<pre>##          bad      good
## 1 0.07142151 0.9285785
## 2 0.04231067 0.9576893
## 3 0.03736701 0.9626330
## 4 0.14796622 0.8520338
## 5 0.12416939 0.8758306
## 6 0.42359516 0.5764048</pre>
<p>If you need to deepen your knowledge of predictive analytics, you may find something interesting in the R Course <a target="_blank" href="http://www.quantide.com/website/courses/data-mining-with-r/">Data Mining with R</a>.</p>
<p>Stay tuned for the next article on the <a target="_blank" href="http://www.milanor.net/blog/">MilanoR blog</a>!</p>
</div>
<div>
<h3>References</h3>
<div>
<div>
<p>Efron, B., and R. Tibshirani. 1993. <em>An Introduction to the Bootstrap</em>. CRC Press.</p>
</div>
<div>
<p>Friedman, J., T. Hastie, and R. Tibshirani. 2008. “Regularization Paths for Generalized Linear Models via Coordinate Descent.” <em>Journal of Statistical Software</em> 33 (1): 1–22.</p>
</div>
<div>
<p>Hastie, T., R. Tibshirani, and J. Friedman. 2009. <em>The Elements of Statistical Learning</em>. 2nd ed. Springer.</p>
</div>
<div>
<p>James, G., D. Witten, T. Hastie, and R. Tibshirani. 2013. <em>An Introduction to Statistical Learning</em>. Springer.</p>
</div>
<div>
<p>Kuhn, M., and K. Johnson. 2013. <em>Applied Predictive Modeling</em>. Springer.</p>
</div>
<div>
<p>Zou, H., and T. Hastie. 2005. “Regularization and Variable Selection via the Elastic Net.” <em>Journal of the Royal Statistical Association B</em> 67 (2): 301–20.</p>
</div>
</div>
</div>
<div>
<hr>
<ol>
<li>By the way, it seems that the oracular powers appeared to be associated with hallucinogenic gases that puffed out from the temple floor.<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fnref1">↩</a></li>
<li>You can find a thorough formal illustration of all these concepts in <span>Hastie, Tibshirani, and Friedman (2009)</span>, Chapter 7. A somewhat simpler presentation can be found in <span>James et al. (2013)</span>.<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fnref2">↩</a></li>
<li>More precisely, the light red curves correspond to what is called <strong>conditional test error</strong>, which means that each curve is conditional on the corresponding training sample. The heavier red curve correspond to the <strong>expected test error</strong>. In general, we would like to focus on the conditional test error for the particular training sample we have. However, this curve is very difficult to be estimated and in practice the expected test error is typically targeted. As we will see, cross-validation is a method for estimating the expected test error. For more details see <span>Hastie, Tibshirani, and Friedman (2009)</span>.<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fnref3">↩</a></li>
<li>A variant of the purely random split is to use stratified random sampling in order to create subsets that are balanced with respect to the outcome. This is useful in particular in classification problems when one class has a disproportionately small frequency compared to the others.<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fnref4">↩</a></li>
<li>An alternative approach for the same objective is the <strong>bootstrap</strong>, that won’t be illustrated here (see <span>Efron and Tibshirani (1993)</span>).<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fnref5">↩</a></li>
<li>More precisely, cross-validation provides an estimate of the <em>expected</em> test error.<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fnref6">↩</a></li>
<li>The <code>boot</code> package contains also a nice function called <code>cv.glm</code>, which implements <span><span data-mathml="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;mi&gt;k&lt;/mi&gt;&lt;/math&gt;" tabindex="0"><span><span><span>k</span></span></span><span>k</span></span></span>-fold cross-validation for generalized linear models.<a target="_blank" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/#fnref7">↩</a></li>
</ol>
</div>
</div>
<p>The post <a target="_blank" rel="nofollow" href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/">Cross-Validation for Predictive Analytics Using R</a> appeared first on <a target="_blank" rel="nofollow" href="http://www.milanor.net/blog">MilanoR</a>.</p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/cross-validation-for-predictive-analytics-using-r/"></div></div><div><a data-url="http://www.r-bloggers.com/cross-validation-for-predictive-analytics-using-r/" data-counturl="http://www.r-bloggers.com/cross-validation-for-predictive-analytics-using-r/" data-count="vertical" data-text="Cross-Validation for Predictive Analytics Using R" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/"> MilanoR</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/swE9qztQOZ4"></a></div>