<p><a 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 href="http://varianceexplained.org/">David Robinson’s</a> <a 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 href="http://rpubs.com/dgrtwo/cv-modelr">his past post</a>, which I read some time ago and felt like unpacking.</p>
<h2><a href="https://drsimonj.svbtle.com/#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 href="https://www.youtube.com/watch?v=nZAM5OXrktY">fantastic tutorial from Stanford University professors Trevor Hastie and Rob Tibshirani</a>.</p>
<h2><a href="https://drsimonj.svbtle.com/#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 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 5, 9, 17, 20, 27, 28, 29, 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> 5, 9, 17, 20, 27, 28, 29:</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 href="https://drsimonj.svbtle.com/#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 href="https://drsimonj.svbtle.com/#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 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 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 href="https://drsimonj.svbtle.com/#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 href="https://svbtleusercontent.com/h7h2tbbfsdlua.png"></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 href="https://svbtleusercontent.com/hmnlslyagqzqa.png"></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 href="https://drsimonj.svbtle.com/#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 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 href="https://svbtleusercontent.com/dwq9es20tiia.png"></a></p>
<h2><a href="https://drsimonj.svbtle.com/#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 href="https://twitter.com/drsimonj">@drsimonj</a> on Twitter, or email me at <a 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 href="https://github.com/drsimonj/blogR">blogR GitHub repository</a>.</p>