<div>
<div>
<p>When estimating a model, the quality of the model fit will always be higher in-sample than out-of-sample. A model will always fit the data that it is trained on, warts and all, and may use those warts and statistical noise to make predictions. As a result, a model that performs very well on a data set may perform poorly when used more generally.</p>
<p>Since we often can’t fit a model, put it in use and determine how well it performs (and probably wouldn’t want to even if we could), we need some methods for predicting out-of-sample error when only using the data we have.</p>
<p>The simplest method is to randomly divide the data into a training set and a test set. However, this provides only one test set. If you train a model, test it and then train a new model and re-test it using the same test set and see an improvement, it may have come from the model actually being better or it may come from you over-fitting the test set.</p>
<p>One possible way to address this problem is to use something like k-fold cross validation. The general idea is is to use some portion of your data (typically 90%) to fit a model and then test that model on the withheld portion (typically 10%) to assess the quality of the model fit. You would repeat this process over and over again until each point has been used to estimate the model <span>\(k-1\)</span> times and been included in the test set 1 time.</p>
<p>This can be done with a loop or a member of the <code>apply</code> family, but I’m going to use <a target="_blank" rel="nofollow" href="https://github.com/hadley/purrr"><code>purrr</code></a> to write functional code to fit the models and <a target="_blank" rel="nofollow" href="https://github.com/hadley/modelr"><code>modelr</code></a> to evaluate the models.</p>
<p>My “preamble” looks like this:</p>
<div>
<h2>Data and the Model</h2>
<p>For this example, I’m going to use the <code>cars</code> dataset from the <code>caret</code> package. It looks like this:</p>
<pre><code>data(&quot;cars&quot;, package = &quot;caret&quot;)
glimpse(cars)</code></pre>
<pre><code>## Observations: 804
## Variables: 18
## $ Price       &lt;dbl> 22661.05, 21725.01, 29142.71, 30731.94, 33358.77, ...
## $ Mileage     &lt;int> 20105, 13457, 31655, 22479, 17590, 23635, 17381, 2...
## $ Cylinder    &lt;int> 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,...
## $ Doors       &lt;int> 4, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 4, 4,...
## $ Cruise      &lt;int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ Sound       &lt;int> 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1,...
## $ Leather     &lt;int> 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1,...
## $ Buick       &lt;int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ Cadillac    &lt;int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ Chevy       &lt;int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ Pontiac     &lt;int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ Saab        &lt;int> 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ Saturn      &lt;int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ convertible &lt;int> 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ coupe       &lt;int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ hatchback   &lt;int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ sedan       &lt;int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ wagon       &lt;int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...</code></pre>
<p>We want to build a model to predict the price of the car. We have the mileage, number of cylinders, size, some features (cruise control, sound system, leather seats), brand and then the style of car. The data already has dummy variables constructed but we need to drop one in order to estimate the model (e.g., set it as the reference level). We’ll drop the brand Chevy indicator (setting Chevy to the reference brand) and set the reference style to sedan. We’ll also drop the doors variable as that is perfectly co-linear with the coupe indicator and styles likely contain more information than the number of doors.</p>
<p>Our final model will look something like this:</p>
<pre><code>select(cars, -Doors, -Chevy, -sedan) %>% 
  lm(Price ~ Milage + Cylinder + Cruise + Sound + Leather + Buick + Cadillac + 
  Pontiac + Saab + Saturn + convertible + coupe + hatchback + wagon, data = .)</code></pre>
<p>The <code>modelr</code> package has a useful tool for making the cross-validation folds. <code>crossv_kfold</code> will divide the data into <span>\(k\)</span> (by default 10) folds and returns a tibble with the list-columns of <code>train</code> (the training data for the fold), <code>test</code> (the test set for the fold) and <code>.id</code> (the fold number). For more details on the idea of list-columns, see the <a target="_blank" rel="nofollow" href="http://r4ds.had.co.nz/many-models.html">Many Models Chapter in R for Data Science</a> by Garrett Grolemund and Hadley Wickham.</p>
<pre><code>cars_cv &lt;- cars %>%
  select(-Doors, -Chevy, -sedan) %>%
  crossv_kfold(10)
cars_cv</code></pre>
<pre><code>## # A tibble: 10 × 3
##             train           test   .id
##            &lt;list>         &lt;list> &lt;chr>
## 1  &lt;S3: resample> &lt;S3: resample>    01
## 2  &lt;S3: resample> &lt;S3: resample>    02
## 3  &lt;S3: resample> &lt;S3: resample>    03
## 4  &lt;S3: resample> &lt;S3: resample>    04
## 5  &lt;S3: resample> &lt;S3: resample>    05
## 6  &lt;S3: resample> &lt;S3: resample>    06
## 7  &lt;S3: resample> &lt;S3: resample>    07
## 8  &lt;S3: resample> &lt;S3: resample>    08
## 9  &lt;S3: resample> &lt;S3: resample>    09
## 10 &lt;S3: resample> &lt;S3: resample>    10</code></pre>
<p>You might notice that instead of the elements being <code>data.frame</code> or <code>tibble</code>, they are <code>resample</code> objects. <code>resample</code> objects are cool because they are lazy. Instead of making 10 copies of <code>cars</code>, <code>crossv_kfold()</code> returned what are essentially pointers to the original data. So instead of storing the data 10 times, we store it once plus a little bit for the pointers. Indeed, we can see that in the object size:</p>
<pre><code>pryr::object_size(cars)</code></pre>
<pre><code>## 63.6 kB</code></pre>
<pre><code>pryr::object_size(cars_cv)</code></pre>
<pre><code>## 93.7 kB</code></pre>
<p>Awesome, right?</p>
<p>Now it is time to fit the model. It all can be down within the very tidy tidyverse (see what I did there?). The code below uses <code>purrr:map()</code> to fit the model for each one of the 10 folds and returns the RMSE on the test set.</p>
<pre><code>cars_cv %>% 
  mutate(train = map(train, as_tibble)) %>%
  mutate(model = map(train, ~ lm(Price ~ Mileage + Cylinder + Cruise + Sound + 
                                   Leather + Buick + Cadillac + Pontiac + 
                                   Saab + Saturn + convertible + coupe + 
                                   hatchback + wagon, 
                                 data = .))) %>%
  mutate(rmse = map2_dbl(model, test, rmse)) %>%
  select(.id, rmse)</code></pre>
<pre><code>## # A tibble: 10 × 2
##      .id     rmse
##    &lt;chr>    &lt;dbl>
## 1     01 2372.046
## 2     02 2844.509
## 3     03 3398.572
## 4     04 3165.065
## 5     05 2838.691
## 6     06 2675.153
## 7     07 3720.843
## 8     08 2692.567
## 9     09 2362.007
## 10    10 3016.648</code></pre>
<p>It is reasonably fast too.</p>
<pre><code>ten_fold &lt;- function(x) {
  x %>% 
    mutate(model = map(train, ~ lm(Price ~ Mileage + Cylinder + Cruise + Sound + 
                                   Leather + Buick + Cadillac + Pontiac + 
                                   Saab + Saturn + convertible + coupe + 
                                   hatchback + wagon, 
                                 data = .))) %>%
  mutate(rmse = map2_dbl(model, test, rmse)) %>%
  select(.id, rmse)
}

microbenchmark::microbenchmark(ten_fold(cars_cv))</code></pre>
<pre><code>## Unit: milliseconds
##               expr      min       lq     mean   median       uq      max
##  ten_fold(cars_cv) 45.47155 50.18167 61.01467 53.27679 60.67367 139.0201
##  neval
##    100</code></pre>
<p><code>microbenchmark</code> puts it a median time of about 50 ms. We can probably make it even faster though. We don’t actually need to store the model and we can do it all in one step with a function.</p>
<pre><code>fitter &lt;- function(training_data, testing_data) {
  lm(Price ~ Mileage + Cylinder + Cruise + Sound + Leather + Buick + Cadillac + 
       Pontiac + Saab + Saturn + convertible + coupe + hatchback + wagon, 
     data = training_data) %>% 
    rmse(testing_data)
}

ten_fold_2 &lt;- function(x) {
  x %>% 
    mutate(rmse = map2_dbl(train, test, fitter)) %>%
    select(.id, rmse)
}

microbenchmark::microbenchmark(ten_fold_2(cars_cv))</code></pre>
<pre><code>## Unit: milliseconds
##                 expr     min       lq     mean   median       uq      max
##  ten_fold_2(cars_cv) 46.1202 49.16919 52.21734 51.25691 53.72752 100.7599
##  neval
##    100</code></pre>
<p>Hey, that saved 2 or 3 milliseconds on the median. Totally worth it. (If your model fit takes a long time, it is probably worth it).</p>
<p>It is possible to parallelize the model fitting using <a target="_blank" rel="nofollow" href="https://github.com/hadley/multidplyr/"><code>multidplyr</code></a> as the backend for <code>dplyr</code>; however, <code>multidplyr</code> doesn’t seem to work well with the <code>resample</code> objects. The easiest solution is to store nine copies of the data to allow for the cross-validation to work most simply - but I’m more likely to hurt for shortage of memory than shortage of CPU cycles/time when dealing with health care administrative data. Alternatively, you could use <code>do()</code> and <code>multidplyr</code> to filter by fold ID, fit a model, compute RMSE and proceed from there, but that solution is much less simple.</p>
<p>Packages like <code>caret</code> do provide powerful tools for cross-validation; however, I like the ability to keep my entire workflow smooth and within the tidyverse. I found that, in general, I prefer to use the <code>modelr</code> tools and do the model fitting the “old fashioned way.”</p>
</div>
</div>

<p></p>
<p><br>
</p>
<p></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="https://www.r-bloggers.com/easy-cross-validation-in-r-with-modelr/"></div></div><div><a data-url="https://www.r-bloggers.com/easy-cross-validation-in-r-with-modelr/" data-counturl="https://www.r-bloggers.com/easy-cross-validation-in-r-with-modelr/" data-count="vertical" data-text="Easy Cross Validation in R with `modelr`" 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://jacobsimmering.com/2016/11/11/CrossValidationInR/"> jacobsimmering.com</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://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="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/nDV3thG0KRo"></a></div></div>