Tidy Resampling Redux with Agricultural Economics Data
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
(No statistical graphs in this one. This is what my dog Artemis looks like when she wants my attention during work hours.)
Mindy L. Mallory (@ace_prof
) wrote a blog post on Machine Learning and Econometrics: Model Selection and Assessment Statistical Learning Style where she has a great description of the variancebias tradeoff, resampling, and model complexity using some data from agricultural economics. I asked if I could take her code and make a more tidy version of the resampling part. Here it is…
Here are the R packages that I’ll be using:
library(tidyverse) library(rlang) library(rsample) library(yardstick) library(purrr)
First, duplicate the process of reading in the data and adding two new columns:
stocks < read.csv('http://blog.mindymallory.com/wpcontent/uploads/2018/02/stocks.csv') %>% as_tibble() %>% mutate( USStockUse = USEndingStocks/USTotalUse, WorldStockUse = ROWEndingStocks/WorldTotalUse ) stocks
## # A tibble: 43 x 8 ## Year USEndingStocks ROWEndingStocks USTotalUse ## <int> <dbl> <int> <dbl> ## 1 1975 633. 36411 5767. ## 2 1976 1136. 39491 5789. ## 3 1977 1436. 40833 6207. ## 4 1978 1710. 47957 6995. ## 5 1979 2034. 59481 7604. ## 6 1980 1392. 67180 7282. ## 7 1981 2537. 62725 6975. ## 8 1982 3523. 60273 7249. ## 9 1983 1006. 63421 6693. ## 10 1984 1648. 76287 7032. ## # ... with 33 more rows, and 4 more variables: ## # WorldTotalUse <int>, PriceRecievedFarmers <dbl>, ## # USStockUse <dbl>, WorldStockUse <dbl>
The blog post has an excellent description of crossvalidation and looked at five different models that encoded the US and World stockuse predictors. Either a log or inversetransformation was applied and then polynomial basis functions were used on these features to demonstrate overfitting.
The blog post has some for
loops to do the resampling and I volunteered to show how to do it with some tidy modeling packages.
Tidy CrossValidation
First, let’s take the easy part. Instead of using for
loops, we can use the new infrastructure in the tidyverse to resample the model. The rsample
package has some functions for different types of resampling and we will use the same procedure as the original post:
set.seed(918) resamp_info < vfold_cv(stocks, v = 5) resamp_info
## # 5fold crossvalidation ## # A tibble: 5 x 2 ## splits id ## <list> <chr> ## 1 <S3: rsplit> Fold1 ## 2 <S3: rsplit> Fold2 ## 3 <S3: rsplit> Fold3 ## 4 <S3: rsplit> Fold4 ## 5 <S3: rsplit> Fold5
The first column in the tibble is a set of "rsplit
" objects that define how the data are split for each fold of crossvalidation. Each one fully and efficiently encapsulates everything what is needed to get the two divisions of the original data. In rsample
, to avoid naming confusion, we label the two resulting data sets as:

The analysis data are those that we selected in the resample. For a bootstrap, this is the sample with replacement. For 5fold crossvalidation, this is the 80% of the data. These data are often used to fit a model or calculate a statistic in traditional bootstrapping.

The assessment data are usually the section of the original data not covered by the analysis set. Again, in 5fold CV, this is the 20% held out. These data are often used to evaluate the performance of a model that was fit to the analysis data.
To get these partitions for the first split there are functions analysis
and assessment
that return the appropriate data frames when given an rsplit
:
# printing just shows the #rows per analysis/assessment/overall resamp_info$splits[[1]]
## <34/9/43>
# data used for modeling: analysis(resamp_info$splits[[1]])
## # A tibble: 34 x 8 ## Year USEndingStocks ROWEndingStocks USTotalUse ## <int> <dbl> <int> <dbl> ## 1 1975 633. 36411 5767. ## 2 1976 1136. 39491 5789. ## 3 1977 1436. 40833 6207. ## 4 1979 2034. 59481 7604. ## 5 1980 1392. 67180 7282. ## 6 1981 2537. 62725 6975. ## 7 1982 3523. 60273 7249. ## 8 1983 1006. 63421 6693. ## 9 1984 1648. 76287 7032. ## 10 1985 4040. 75069 6494. ## # ... with 24 more rows, and 4 more variables: ## # WorldTotalUse <int>, PriceRecievedFarmers <dbl>, ## # USStockUse <dbl>, WorldStockUse <dbl>
The first model for the data contained the US stockuse data with an inverse transformation. Let's sidestep the polynomial model tuning for now and just fit a quadratic model. To make things easier, I'll define a function that can be used to fit the model when given an rsplit
object and return the holdout mean squared error (MSE):
glm_results < function(split, ...) { # Get the data used ot fit the model aka the "analysis" set # and fit the model with a formula given in the ... mod < glm(data = analysis(split), ...) # Get predictions on the other data (aka the "assessment" set # and compute some metrics holdout < assessment(split) # Compute performance using the yardstick package rmse < holdout %>% mutate(pred = predict(mod, holdout)) %>% rmse(truth = PriceRecievedFarmers, estimate = pred) rmse^2 }
We can use this with any formula since it is just passed to glm
using the ellipses. For example to get the holdout MSE estimate for the first fold:
glm_results(resamp_info$splits[[1]], formula = PriceRecievedFarmers ~ poly(1 / USStockUse, 2))
## [1] 1.67
To get these statistics for all folds, purrr::map_dbl
is used to add another column:
resamp_info < resamp_info %>% mutate( model_1_deg_2 = map_dbl( splits, glm_results, PriceRecievedFarmers ~ poly(1 / USStockUse, 2) ) ) resamp_info
## # 5fold crossvalidation ## # A tibble: 5 x 3 ## splits id model_1_deg_2 ## * <list> <chr> <dbl> ## 1 <S3: rsplit> Fold1 1.67 ## 2 <S3: rsplit> Fold2 0.558 ## 3 <S3: rsplit> Fold3 0.596 ## 4 <S3: rsplit> Fold4 9.58 ## 5 <S3: rsplit> Fold5 0.248
That's a lot of variation in the outcome! The mean value is fairly consistent with the blog post though:
resamp_info %>% select(model_1_deg_2) %>% colMeans()
## model_1_deg_2 ## 2.53
# MSE = 2.675 in the blog post
Kfold crossvalidation is one of the noisiest resampling methods so this difference isn't too surprising.
This same process could be repeated for each polynomial degree to get new columns for this model (we'll discuss this below). The good things about doing things this way:
 It is a lot cleaner (so far) than doing
for
loops.  Other tidyverse infrastructure can be used. For example,
tidyposterior
is a great way to do model comparisons with resampling and Bayesian analysis.  It is simple to change resampling methods. Suppose you wanted to change to a larger number of bootstrap resamples (given the variance shown above). The same infrastructure can be easily exchanged;
resample::bootstraps
is used in place ofrsample::vfold_cv
.
Tidy Model Specification (maybe)
The model specification part is, for me, a lot more difficult to tidy. It would be good to be able to state what predictors that we want, specify the polynomial degree, and have a function to generate the appropriate formula. The original post sensibly just types the terms out.
I spent some time thinking about how we could use expressions and tidy evaluation (video from Hadley) to make it a little less scriptlike. The problem is that the solution took me a while to write and, arguably, it doesn't really buy you much more than the original code (apart from the potential copy/paste duplication errors).
In any case, the function uses rlang
to manipulate expressions to make the formula. The inputs are expressions of how the predictors are used (i.e. log, inverse, etc.) and the degree. It captures the expression (without evaluating it), substitutes it into the polynomial function, then creates the formula:
make_formula < function(..., degree = 1) { # Capture the expression so it is not evaluated var_expr < exprs(...) # Create a template expression inv_poly < quote(poly(x = x, degree = degree)) # Use a wrapper around rlang::call_modify to reverse the # order of the arguments so that we can map over the # predictor expressions add_args < function(arg, call, ...) call_modify(call, x = arg, ...) # Add the variables and the degree into the template poly_expr < map(var_expr, add_args, call = inv_poly, degree = degree) # Convert to character poly_char < map_chr(poly_expr, deparse) # Convert to a formula poly_char < paste(poly_char, collapse = " + ") as.formula(paste("PriceRecievedFarmers ~", poly_char)) } # For example: make_formula(1/USStockUse, log(WorldStockUse), degree = 3)
## PriceRecievedFarmers ~ poly(x = 1/USStockUse, degree = 3) + poly(x = log(WorldStockUse), ## degree = 3) ## <environment: 0x7fdf40b64840>
From here, we could use a bunch of mutate
commands like the one shown above or write a slightly smaller for
loop to work across the polynomial degrees. While the function above works well, the overall approach to working across models isn't particularly satisfying.
Any suggestions? I'm on twitter @topepos
!
(edit  fix fixed the formula call!)
Rbloggers.com offers daily email updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.