<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/k-fold-cross-validation-with-modelr-and-broom/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom"> blogR</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><a target="_blank" rel="nofollow" href="https://twitter.com/drsimonj">@drsimonj</a> here to discuss how to conduct k-fold cross validation, with an emphasis on evaluating models supported by <a target="_blank" rel="nofollow" href="http://varianceexplained.org/">David Robinson’s</a> <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/broom/index.html">broom</a> package. Full credit also goes to David, as this is a slightly more detailed version of <a target="_blank" rel="nofollow" href="http://rpubs.com/dgrtwo/cv-modelr">his past post</a>, which I read some time ago and felt like unpacking.</p>
<h2><a target="_blank" rel="nofollow" href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#assumed-knowledge-kfold-cross-validation_2"> </a>Assumed knowledge: K-fold Cross validation</h2>
<p>This post assumes you know what k-fold cross validation is. If you want to brush up, here’s a <a target="_blank" rel="nofollow" href="https://www.youtube.com/watch?v=nZAM5OXrktY">fantastic tutorial from Stanford University professors Trevor Hastie and Rob Tibshirani</a>.</p>
<h2><a target="_blank" rel="nofollow" href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#creating-folds_2"> </a>Creating folds</h2>
<p>Before worrying about models, we can generate K folds using <code>crossv_kfold</code> from the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/modelr/index.html">modelr</a> package. Let’s practice with the <code>mtcars</code> data to keep things simple.</p>
<pre><code>library(modelr)
set.seed(1)  # Run to replicate this post
folds &lt;- crossv_kfold(mtcars, k = 5)
folds
#> # A tibble: 5 × 3
#>            train           test   .id
#>           &lt;list>         &lt;list> &lt;chr>
#> 1 &lt;S3: resample> &lt;S3: resample>     1
#> 2 &lt;S3: resample> &lt;S3: resample>     2
#> 3 &lt;S3: resample> &lt;S3: resample>     3
#> 4 &lt;S3: resample> &lt;S3: resample>     4
#> 5 &lt;S3: resample> &lt;S3: resample>     5
</code></pre>
<p>This function takes a data frame and randomly partitions it’s rows (1 to 32 for <code>mtcars</code>) into <code>k</code> roughly equal groups. We’ve partitioned the row numbers into <code>k = 5</code> groups. The results are returned as a tibble (data frame) like the one above.</p>
<p>Each cell in the <code>test</code> column contains a <code>resample</code> object, which is an efficient way of referencing a subset of rows in a data frame (<code>?resample</code> to learn more). Think of each cell as a reference to the rows of the data frame belonging to each partition. For example, the following tells us that the first partition of the data references rows <em>5, 9, 17, 20, 27, 28, 29</em>, which accounts for roughly <code>1 / k</code> of the total data set (7 of the 32 rows).</p>
<pre><code>folds$test[[1]]
#> &lt;resample [7 x 11]> 5, 9, 17, 20, 27, 28, 29
</code></pre>
<p>Each cell in <code>train</code> also contains a <code>resample</code> object, but referencing the rows in all other partitions. For example, the first <code>train</code> object references all rows <strong>except</strong> <em>5, 9, 17, 20, 27, 28, 29</em>:</p>
<pre><code>folds$train[[1]]
#> &lt;resample [25 x 11]> 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, ...
</code></pre>
<p>We can now run a model on the data referenced by each <code>train</code> object, and validate the model results on each corresponding partition in <code>test</code>.</p>
<h2><a target="_blank" rel="nofollow" href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#fitting-models-to-training-data_2"> </a>Fitting models to training data</h2>
<p>Say we’re interested in predicting Miles Per Gallon (<code>mpg</code>) with all other variables. With the whole data set, we’d do this via:</p>
<pre><code>lm(mpg ~ ., data = mtcars)
</code></pre>
<p>Instead, we want to run this model on each set of training data (data referenced in each <code>train</code> cell). We can do this as follows:</p>
<pre><code>library(dplyr)
library(purrr)

folds &lt;- folds %>% mutate(model = map(train, ~ lm(mpg ~ ., data = .)))
folds
#> # A tibble: 5 × 4
#>            train           test   .id    model
#>           &lt;list>         &lt;list> &lt;chr>   &lt;list>
#> 1 &lt;S3: resample> &lt;S3: resample>     1 &lt;S3: lm>
#> 2 &lt;S3: resample> &lt;S3: resample>     2 &lt;S3: lm>
#> 3 &lt;S3: resample> &lt;S3: resample>     3 &lt;S3: lm>
#> 4 &lt;S3: resample> &lt;S3: resample>     4 &lt;S3: lm>
#> 5 &lt;S3: resample> &lt;S3: resample>     5 &lt;S3: lm>
</code></pre>
<ul>
<li> <code>folds %>% mutate(model = ...)</code> is adding a new <code>model</code> column to the folds tibble.</li>
<li> <code>map(train, ...)</code> is applying a function to each of the cells in <code>train</code></li>
<li> <code>~ lm(...)</code> is the regression model applied to each <code>train</code> cell.</li>
<li> <code>data = .</code> specifies that the data for the regression model will be the data referenced by each <code>train</code> object.</li>
</ul>
<p>The result is a new <code>model</code> column containing fitted regression models based on each of the <code>train</code> data (i.e., the whole data set excluding each partition).</p>
<p>For example, the model fitted to our first set of training data is:</p>
<pre><code>folds$model[[1]] %>% summary()
#> 
#> Call:
#> lm(formula = mpg ~ ., data = .)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.6540 -0.9116  0.0439  0.9520  4.2811 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)  
#> (Intercept) -44.243933  31.884363  -1.388   0.1869  
#> cyl           0.844966   1.064141   0.794   0.4404  
#> disp          0.016800   0.015984   1.051   0.3110  
#> hp            0.004685   0.022741   0.206   0.8398  
#> drat          3.950410   1.989177   1.986   0.0670 .
#> wt           -4.487007   2.016341  -2.225   0.0430 *
#> qsec          2.327131   1.243095   1.872   0.0822 .
#> vs           -3.963492   3.217176  -1.232   0.2382  
#> am           -0.550804   2.333252  -0.236   0.8168  
#> gear          5.476604   2.648708   2.068   0.0577 .
#> carb         -1.595979   1.104272  -1.445   0.1704  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2.159 on 14 degrees of freedom
#> Multiple R-squared:  0.9092, Adjusted R-squared:  0.8444 
#> F-statistic: 14.02 on 10 and 14 DF,  p-value: 1.205e-05
</code></pre>
<h2><a target="_blank" rel="nofollow" href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#predicting-the-test-data_2"> </a>Predicting the test data</h2>
<p>The next step is to use each model for predicting the outcome variable in the corresponding <code>test</code> data. There are many ways to achieve this. One general approach might be:</p>
<pre><code>folds %>% mutate(predicted = map2(model, test, &lt;FUNCTION_TO_PREDICT_TEST_DATA> ))
</code></pre>
<p><code>map2(model, test, ...)</code> iterates through each model and set of <code>test</code> data in parallel. By referencing these in the function for predicting the test data, this would add a <code>predicted</code> column with the predicted results.</p>
<p>For many common models, an elegant alternative is to use <code>augment</code> from <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/broom/index.html">broom</a>. For regression, <code>augment</code> will take a fitted model and a new data frame, and return a data frame of the predicted results, which is what we want! Following above, we can use <code>augment</code> as follows:</p>
<pre><code>library(broom)

folds %>% mutate(predicted = map2(model, test, ~ augment(.x, newdata = .y)))
#> # A tibble: 5 × 5
#>            train           test   .id    model             predicted
#>           &lt;list>         &lt;list> &lt;chr>   &lt;list>                &lt;list>
#> 1 &lt;S3: resample> &lt;S3: resample>     1 &lt;S3: lm> &lt;data.frame [7 × 13]>
#> 2 &lt;S3: resample> &lt;S3: resample>     2 &lt;S3: lm> &lt;data.frame [7 × 13]>
#> 3 &lt;S3: resample> &lt;S3: resample>     3 &lt;S3: lm> &lt;data.frame [6 × 13]>
#> 4 &lt;S3: resample> &lt;S3: resample>     4 &lt;S3: lm> &lt;data.frame [6 × 13]>
#> 5 &lt;S3: resample> &lt;S3: resample>     5 &lt;S3: lm> &lt;data.frame [6 × 13]>
</code></pre>
<p>To extract the relevant information from these <code>predicted</code> results, we’ll <code>unnest</code> the data frames thanks to the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/tidyr/index.html">tidyr</a> package:</p>
<pre><code>library(tidyr)

folds %>%
  mutate(predicted = map2(model, test, ~ augment(.x, newdata = .y))) %>% 
  unnest(predicted)
#> # A tibble: 32 × 14
#>      .id   mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
#>    &lt;chr> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl>
#> 1      1  18.7     8 360.0   175  3.15 3.440 17.02     0     0     3     2
#> 2      1  22.8     4 140.8    95  3.92 3.150 22.90     1     0     4     2
#> 3      1  14.7     8 440.0   230  3.23 5.345 17.42     0     0     3     4
#> 4      1  33.9     4  71.1    65  4.22 1.835 19.90     1     1     4     1
#> 5      1  26.0     4 120.3    91  4.43 2.140 16.70     0     1     5     2
#> 6      1  30.4     4  95.1   113  3.77 1.513 16.90     1     1     5     2
#> 7      1  15.8     8 351.0   264  4.22 3.170 14.50     0     1     5     4
#> 8      2  21.0     6 160.0   110  3.90 2.875 17.02     0     1     4     4
#> 9      2  21.4     6 258.0   110  3.08 3.215 19.44     1     0     3     1
#> 10     2  24.4     4 146.7    62  3.69 3.190 20.00     1     0     4     2
#> # ... with 22 more rows, and 2 more variables: .fitted &lt;dbl>,
#> #   .se.fit &lt;dbl>
</code></pre>
<p>This was to show you the intermediate steps. In practice we can skip the <code>mutate</code> step:</p>
<pre><code>predicted &lt;- folds %>% unnest(map2(model, test, ~ augment(.x, newdata = .y)))
predicted
#> # A tibble: 32 × 14
#>      .id   mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
#>    &lt;chr> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl>
#> 1      1  18.7     8 360.0   175  3.15 3.440 17.02     0     0     3     2
#> 2      1  22.8     4 140.8    95  3.92 3.150 22.90     1     0     4     2
#> 3      1  14.7     8 440.0   230  3.23 5.345 17.42     0     0     3     4
#> 4      1  33.9     4  71.1    65  4.22 1.835 19.90     1     1     4     1
#> 5      1  26.0     4 120.3    91  4.43 2.140 16.70     0     1     5     2
#> 6      1  30.4     4  95.1   113  3.77 1.513 16.90     1     1     5     2
#> 7      1  15.8     8 351.0   264  4.22 3.170 14.50     0     1     5     4
#> 8      2  21.0     6 160.0   110  3.90 2.875 17.02     0     1     4     4
#> 9      2  21.4     6 258.0   110  3.08 3.215 19.44     1     0     3     1
#> 10     2  24.4     4 146.7    62  3.69 3.190 20.00     1     0     4     2
#> # ... with 22 more rows, and 2 more variables: .fitted &lt;dbl>,
#> #   .se.fit &lt;dbl>
</code></pre>
<p>We now have a tibble of the <code>test</code> data for each fold (<code>.id</code> = fold number) and the corresponding <code>.fitted</code>, or predicted values for the outcome variable (<code>mpg</code>) in each case.</p>
<h2><a target="_blank" rel="nofollow" href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#validating-the-model_2"> </a>Validating the model</h2>
<p>We can compute and examine the residuals:</p>
<pre><code># Compute the residuals
predicted &lt;- predicted %>% 
  mutate(residual = .fitted - mpg)

# Plot actual v residual values
library(ggplot2)
predicted %>%
  ggplot(aes(mpg, residual)) +
    geom_hline(yintercept = 0) +
    geom_point() +
    stat_smooth(method = &quot;loess&quot;) +
    theme_minimal()
</code></pre>
<p><a target="_blank" rel="nofollow" href="https://i1.wp.com/svbtleusercontent.com/h7h2tbbfsdlua.png?ssl=1"><img alt="unnamed-chunk-12-1.png" src="https://i2.wp.com/svbtleusercontent.com/h7h2tbbfsdlua_small.png?w=456&ssl=1" data-recalc-dims="1"></a></p>
<p>It looks like our models could be overestimating <code>mpg</code> around 20-30 and underestimating higher <code>mpg</code>. But there are clearly fewer data points, making prediction difficult.</p>
<p>We can also use these values to calculate the overall proportion of variance accounted for by each model:</p>
<pre><code>rs &lt;- predicted %>%
  group_by(.id) %>% 
  summarise(
    sst = sum((mpg - mean(mpg)) ^ 2), # Sum of Squares Total
    sse = sum(residual ^ 2),          # Sum of Squares Residual/Error
    r.squared = 1 - sse / sst         # Proportion of variance accounted for
    )
rs
#> # A tibble: 5 × 4
#>     .id      sst       sse r.squared
#>   &lt;chr>    &lt;dbl>     &lt;dbl>     &lt;dbl>
#> 1     1 321.5886 249.51370 0.2241214
#> 2     2 127.4371  31.86994 0.7499164
#> 3     3 202.6600  42.19842 0.7917773
#> 4     4 108.4733  50.79684 0.5317113
#> 5     5 277.3283  59.55946 0.7852385

# Overall
mean(rs$r.squared)
#> [1] 0.616553
</code></pre>
<p>So, across the folds, the regression models have accounted for an average of 61.66% of the variance of <code>mpg</code> in new, test data.</p>
<p>Plotting these results:</p>
<pre><code>rs %>% 
  ggplot(aes(r.squared, fill  = .id)) +
    geom_histogram() +
    geom_vline(aes(xintercept = mean(r.squared)))  # Overall mean
</code></pre>
<p><a target="_blank" rel="nofollow" href="https://i2.wp.com/svbtleusercontent.com/hmnlslyagqzqa.png?ssl=1"><img alt="unnamed-chunk-14-1.png" src="https://i2.wp.com/svbtleusercontent.com/hmnlslyagqzqa_small.png?w=456&ssl=1" data-recalc-dims="1"></a></p>
<p>Although the model performed well on average, it performed pretty poorly when fold 1 was used as test data.</p>
<h2><a target="_blank" rel="nofollow" href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#all-at-once_2"> </a>All at once</h2>
<p>With this new knowledge, we can do something similar to the <code>k = 20</code> case shown in <a target="_blank" rel="nofollow" href="http://rpubs.com/dgrtwo/cv-modelr">David’s post</a>. See that you can understand what’s going on here:</p>
<pre><code>set.seed(1)
# Select four variables from the mpg data set in ggplot2
ggplot2::mpg %>% select(year, cyl, drv, hwy) %>% 
  # Create 20 folds (5% of the data in each partition)
  crossv_kfold(k = 20) %>%
  # Fit a model to training data
  mutate(model = map(train, ~ lm(hwy ~ ., data = .))) %>%
  # Unnest predicted values on test data
  unnest(map2(model, test, ~ augment(.x, newdata = .y))) %>% 
  # Compute R-squared values for each partition
  group_by(.id) %>%
  summarise(
    sst = sum((hwy - mean(hwy)) ^ 2),
    sse = sum((hwy - .fitted) ^ 2),
    r.squared = 1 - sse / sst
  ) %>% 
  # Plot
  ggplot(aes(r.squared)) +
    geom_density() +
    geom_vline(aes(xintercept = mean(r.squared))) +
    theme_minimal()
</code></pre>
<p><a target="_blank" rel="nofollow" href="https://i0.wp.com/svbtleusercontent.com/dwq9es20tiia.png?ssl=1"><img alt="unnamed-chunk-15-1.png" src="https://i1.wp.com/svbtleusercontent.com/dwq9es20tiia_small.png?w=456&ssl=1" data-recalc-dims="1"></a></p>
<h2><a target="_blank" rel="nofollow" href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#sign-off_2"> </a>Sign off</h2>
<p>Thanks for reading and I hope this was useful for you.</p>
<p>For updates of recent blog posts, follow <a target="_blank" rel="nofollow" href="https://twitter.com/drsimonj">@drsimonj</a> on Twitter, or email me at <a target="_blank" rel="nofollow" href="mailto:drsimonjackson@gmail.com">drsimonjackson@gmail.com</a> to get in touch.</p>
<p>If you’d like the code that produced this blog, check out the <a target="_blank" rel="nofollow" href="https://github.com/drsimonj/blogR">blogR GitHub repository</a>.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom"> blogR</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/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/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/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/VxpOYmkCbr8"></a></div>