k-fold cross validation with modelr and broom

一时失言乱红尘 2022-07-14 14:09 139阅读 0赞

## Assumed knowledge: K-fold Cross validation ##

This post assumes you know what k-fold cross validation is. If you want to brush up, here’s a [fantastic tutorial from Stanford University professors Trevor Hastie and Rob Tibshirani][].

## [ ][Link 1]Creating folds ##

Before worrying about models, we can generate K folds using`crossv_kfold` from the [modelr][] package. Let’s practice with the `mtcars`data to keep things simple.

library(modelr)
    set.seed(1)  # Run to replicate this post
    folds <- crossv_kfold(mtcars, k = 5)
    folds
    #> # A tibble: 5 × 3
    #> train test .id
    #> <list> <list> <chr>
    #> 1 <S3: resample> <S3: resample> 1
    #> 2 <S3: resample> <S3: resample> 2
    #> 3 <S3: resample> <S3: resample> 3
    #> 4 <S3: resample> <S3: resample> 4
    #> 5 <S3: resample> <S3: resample> 5

This function takes a data frame and randomly partitions it’s rows (1 to 32 for`mtcars`) into `k` roughly equal groups. We’ve partitioned the row numbers into `k = 5` groups. The results are returned as a tibble (data frame) like the one above.

Each cell in the `test` column contains a `resample` object, which is an efficient way of referencing a subset of rows in a data frame (`?resample` 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 `1 / k` of the total data set (7 of the 32 rows).

folds$test[[1]]
    #> <resample [7 x 11]> 5, 9, 17, 20, 27, 28, 29

Each cell in `train` also contains a `resample` object, but referencing the rows in all other partitions. For example, the first `train` object references all rows except *5, 9, 17, 20, 27, 28, 29*:

folds$train[[1]]
    #> <resample [25 x 11]> 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, ...

We can now run a model on the data referenced by each `train` object, and validate the model results on each corresponding partition in `test`.

## [ ][Link 2]Fitting models to training data ##

Say we’re interested in predicting Miles Per Gallon (`mpg`) with all other variables. With the whole data set, we’d do this via:

lm(mpg ~ ., data = mtcars)

Instead, we want to run this model on each set of training data (data referenced in each `train` cell). We can do this as follows:

library(dplyr)
    library(purrr)
    
    folds <- folds %>% mutate(model = map(train, ~ lm(mpg ~ ., data = .)))
    folds
    #> # A tibble: 5 × 4
    #> train test .id model
    #> <list> <list> <chr> <list>
    #> 1 <S3: resample> <S3: resample> 1 <S3: lm>
    #> 2 <S3: resample> <S3: resample> 2 <S3: lm>
    #> 3 <S3: resample> <S3: resample> 3 <S3: lm>
    #> 4 <S3: resample> <S3: resample> 4 <S3: lm>
    #> 5 <S3: resample> <S3: resample> 5 <S3: lm>

*  `folds %>% mutate(model = ...)` is adding a new `model` column to the folds tibble.
 *  `map(train, ...)` is applying a function to each of the cells in `train`
 *  `~ lm(...)` is the regression model applied to each `train` cell.
 *  `data = .` specifies that the data for the regression model will be the data referenced by each `train` object.

The result is a new `model` column containing fitted regression models based on each of the `train` data (i.e., the whole data set excluding each partition).

For example, the model fitted to our first set of training data is:

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

## [ ][Link 3]Predicting the test data ##

The next step is to use each model for predicting the outcome variable in the corresponding `test` data. There are many ways to achieve this. One general approach might be:

folds %>% mutate(predicted = map2(model, test, <FUNCTION_TO_PREDICT_TEST_DATA> ))

`map2(model, test, ...)` iterates through each model and set of `test`data in parallel. By referencing these in the function for predicting the test data, this would add a `predicted` column with the predicted results.

For many common models, an elegant alternative is to use `augment` from[broom][]. For regression, `augment` 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 `augment` as follows:

library(broom)
    
    folds %>% mutate(predicted = map2(model, test, ~ augment(.x, newdata = .y)))
    #> # A tibble: 5 × 5
    #> train test .id model predicted
    #> <list> <list> <chr> <list> <list>
    #> 1 <S3: resample> <S3: resample> 1 <S3: lm> <data.frame [7 × 13]>
    #> 2 <S3: resample> <S3: resample> 2 <S3: lm> <data.frame [7 × 13]>
    #> 3 <S3: resample> <S3: resample> 3 <S3: lm> <data.frame [6 × 13]>
    #> 4 <S3: resample> <S3: resample> 4 <S3: lm> <data.frame [6 × 13]>
    #> 5 <S3: resample> <S3: resample> 5 <S3: lm> <data.frame [6 × 13]>

To extract the relevant information from these `predicted` results, we’ll`unnest` the data frames thanks to the [tidyr][] package:

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
    #> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <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 <dbl>,
    #> # .se.fit <dbl>

This was to show you the intermediate steps. In practice we can skip the`mutate` step:

predicted <- 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
    #> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <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 <dbl>,
    #> # .se.fit <dbl>

We now have a tibble of the `test` data for each fold (`.id` = fold number) and the corresponding `.fitted`, or predicted values for the outcome variable (`mpg`) in each case.

## [ ][Link 4]Validating the model ##

We can compute and examine the residuals:

# Compute the residuals
    predicted <- 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 = "loess") +
        theme_minimal()

[![unnamed-chunk-12-1.png][]][unnamed-chunk-12-1.png 1]

It looks like our models could be overestimating `mpg` around 20-30 and underestimating higher `mpg`. But there are clearly fewer data points, making prediction difficult.

We can also use these values to calculate the overall proportion of variance accounted for by each model:

rs <- 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
    #> <chr> <dbl> <dbl> <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

So, across the folds, the regression models have accounted for an average of 61.66% of the variance of `mpg` in new, test data.

Plotting these results:

rs %>% 
      ggplot(aes(r.squared, fill  = .id)) +
        geom_histogram() +
        geom_vline(aes(xintercept = mean(r.squared)))  # Overall mean

[![unnamed-chunk-14-1.png][]][unnamed-chunk-14-1.png 1]

Although the model performed well on average, it performed pretty poorly when fold 1 was used as test data.

## [ ][Link 5]All at once ##

With this new knowledge, we can do something similar to the `k = 20` case shown in [David’s post][David_s post]. See that you can understand what’s going on here:

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()

[![unnamed-chunk-15-1.png][]][unnamed-chunk-15-1.png 1]

## [ ][Link 6]Sign off ##

Thanks for reading and I hope this was useful for you.

For updates of recent blog posts, follow [@drsimonj][drsimonj] on Twitter, or email me at [drsimonjackson@gmail.com][drsimonjackson_gmail.com] to get in touch.

If you’d like the code that produced this blog, check out the [blogR GitHub repository][].

[fantastic tutorial from Stanford University professors Trevor Hastie and Rob Tibshirani]: https://www.youtube.com/watch?v=nZAM5OXrktY
[Link 1]: https://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#creating-folds_2
[modelr]: https://cran.r-project.org/web/packages/modelr/index.html
[Link 2]: https://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#fitting-models-to-training-data_2
[Link 3]: https://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#predicting-the-test-data_2
[broom]: https://cran.r-project.org/web/packages/broom/index.html
[tidyr]: https://cran.r-project.org/web/packages/tidyr/index.html
[Link 4]: https://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#validating-the-model_2
[unnamed-chunk-12-1.png]: /images/20220714/911dceeb5555400abe7d8b645fb08b6f.png
[unnamed-chunk-12-1.png 1]: https://svbtleusercontent.com/h7h2tbbfsdlua.png
[unnamed-chunk-14-1.png]: /images/20220714/413757ee990a4b239f4005f106af832b.png
[unnamed-chunk-14-1.png 1]: https://svbtleusercontent.com/hmnlslyagqzqa.png
[Link 5]: https://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#all-at-once_2
[David_s post]: http://rpubs.com/dgrtwo/cv-modelr
[unnamed-chunk-15-1.png]: /images/20220714/8f183367a00342259148bd8924bb9805.png
[unnamed-chunk-15-1.png 1]: https://svbtleusercontent.com/dwq9es20tiia.png
[Link 6]: https://drsimonj.svbtle.com/k-fold-cross-validation-with-modelr-and-broom#sign-off_2
[drsimonj]: https://twitter.com/drsimonj
[drsimonjackson_gmail.com]: mailto:drsimonjackson@gmail.com
[blogR GitHub repository]: https://github.com/drsimonj/blogR