As I’ve discussed previously, we sometimes don’t have enough data where doing a train/test split makes sense. As such, we are better off building our model using cross-validation. In previous blog articles, I’ve talked about how to build models using cross-validation within the {**tidymodels**} framework (see **HERE** and **HERE**). In my prior examples, we fit the model over the cross-validation folds and then constructed the final model that we could then use to make predictions with, later on.

Recently, I ran into a situation where I wanted to see what the model coefficients look like across all of the cross-validation folds. So, I decided to make a quick blog post on how to do this, in case it is useful to others.

**Load Packages & Data**

We will use the {**mtcars**} package from R and build a regression model, using several independent variables, to predict miles per gallon (mpg).

### Packages ------------------------------------------------------- library(tidyverse) library(tidymodels) ### Data ------------------------------------------------------- dat <- mtcars dat %>% head()

**Create Cross-Validation Folds of the Data**

I’ll use 10-fold cross validation.

### Modelling ------------------------------------------------------- ## Create 10 Cross Validation Folds set.seed(1) cv_folds <- vfold_cv(dat, v = 10) cv_folds

**Specify a linear model and set up the model formula**

## Specify the linear regression engine ## model specs lm_spec <- linear_reg() %>% set_engine("lm") ## Model formula mpg_formula <- mpg ~ cyl + disp + wt + drat

**Set up the model workflow and fit the model to the cross-validated folds**

## Set up workflow lm_wf <- workflow() %>% add_formula(mpg_formula) %>% add_model(lm_spec) ## Fit the model to the cross validation folds lm_fit <- lm_wf %>% fit_resamples( resamples = cv_folds, control = control_resamples(extract = extract_model, save_pred = TRUE) )

**Extract the model coefficients for each of the 10 folds (this is the fun part!)**

Looking at the **lm_fit** output above, we see that it is a tibble consisting of various nested lists. The **id** column indicates which cross-validation fold the lists in each row pertain to. The model coefficients for each fold are stored in the **.extracts** column of lists. Instead of printing out all 10, let’s just have a look at the first 3 folds to see what they look like.

lm_fit$.extracts %>% .[1:3]

There we see in the **.extracts** column, **<lm****> **indicating the linear model for each fold. With a series of unnesting we can snag the model coefficients and then put them into a tidy format using the {**broom**} package. I’ve commented out each line of code below so that you know exactly what is happening.

# Let's unnest this and get the coefficients out model_coefs <- lm_fit %>% select(id, .extracts) %>% # get the id and .extracts columns unnest(cols = .extracts) %>% # unnest .extracts, which produces the model in a list mutate(coefs = map(.extracts, tidy)) %>% # use map() to apply the tidy function and get the coefficients in their own column unnest(coefs) # unnest the coefs column you just made to get the coefficients for each fold model_coefs

Now that we have a table of estimates, we can plot the coefficient estimates and their 95% confidence intervals. The **term** column indicates each variable. We will remove the **(Intercept)** for plotting purposes.

**Plot the Coefficients**

## Plot the model coefficients and 2*SE across all folds model_coefs %>% filter(term != "(Intercept)") %>% select(id, term, estimate, std.error) %>% group_by(term) %>% mutate(avg_estimate = mean(estimate)) %>% ggplot(aes(x = id, y = estimate)) + geom_hline(aes(yintercept = avg_estimate), size = 1.2, linetype = "dashed") + geom_point(size = 4) + geom_errorbar(aes(ymin = estimate - 2*std.error, ymax = estimate + 2*std.error), width = 0.1, size = 1.2) + facet_wrap(~term, scales = "free_y") + labs(x = "CV Folds", y = "Estimate ± 95% CI", title = "Regression Coefficients ± 95% CI for 10-fold CV", subtitle = "Dashed Line = Average Coefficient Estimate over 10 CV Folds per Independent Variable") + theme_classic() + theme(strip.background = element_rect(fill = "black"), strip.text = element_text(face = "bold", size = 12, color = "white"), axis.title = element_text(size = 14, face = "bold"), axis.text.x = element_text(angle = 60, hjust = 1, face = "bold", size = 12), axis.text.y = element_text(face = "bold", size = 12), plot.title = element_text(size = 18), plot.subtitle = element_text(size = 16))

Now we can clearly see the model coefficients and confidence intervals for each of the 10 cross validated folds.

**Wrapping Up**

This was just a quick and easy way of fitting a model using cross-validation to extract out the model coefficients for each fold. Often, this is probably not necessary as you will fit your model, evaluate your model, and be off and running. However, there may be times where more specific interrogation of the model is required or, you might want to dig a little deeper into the various outputs of the cross-validated folds.

All of the code is available on my **GitHub page**.

If you notice any errors in code, please reach out!