# Estimate change in #TidyTuesday CEO departures with bootstrap resampling

**rstats | Julia Silge**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This is the latest in my series of
screencasts demonstrating how to use the
tidymodels packages, from just starting out to tuning more complex models with many hyperparameters. Today’s screencast walks through how to use bootstrap resampling, with this week’s
`#TidyTuesday`

dataset on CEO departures. 👋

Here is the code I used in the video, for those who prefer reading instead of or in addition to video.

## Explore data

Our modeling goal is to estimate how involuntary CEO departures are changing with time. Let’s start by reading in the data.

library(tidyverse) departures_raw <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-04-27/departures.csv")

How are involuntary departures changing with time? What about the rest of the CEO departures?

departures_raw %>% filter(departure_code < 9) %>% mutate(involuntary = if_else(departure_code %in% 3:4, "involuntary", "other")) %>% filter(fyear > 1995, fyear < 2019) %>% count(fyear, involuntary) %>% ggplot(aes(fyear, n, color = involuntary)) + geom_line(size = 1.2, alpha = 0.5) + geom_point(size = 2) + geom_smooth(method = "lm", lty = 2) + scale_y_continuous(limits = c(0, NA)) + labs(x = NULL, y = "Number of CEO departures", color = NULL)

Looks like proportionally more departures are involuntary over time, but that is what we’ll work on estimating. Let’s create a data set to use for modeling.

departures <- departures_raw %>% filter(departure_code < 9) %>% mutate(involuntary = if_else(departure_code %in% 3:4, "involuntary", "other")) %>% filter(fyear > 1995, fyear < 2019) departures ## # A tibble: 6,942 x 20 ## dismissal_datase… coname gvkey fyear co_per_rol exec_fullname departure_code ## <dbl> <chr> <dbl> <dbl> <dbl> <chr> <dbl> ## 1 559043 SONICB… 27903 2002 -1 L. Gregory B… 7 ## 2 12 AMERIC… 1045 1997 1 Robert L. Cr… 5 ## 3 13 AMERIC… 1045 2002 3 Donald J. Ca… 3 ## 4 31 ABBOTT… 1078 1998 6 Duane L. Bur… 5 ## 5 43 ADVANC… 1161 2001 11 Walter Jerem… 5 ## 6 51 AETNA … 1177 1997 16 Ronald Edwar… 5 ## 7 63 AHMANS… 1194 1997 22 Charles R. R… 7 ## 8 65 AIR PR… 1209 2000 28 Harold A. Wa… 5 ## 9 76 ALBERT… 1239 2007 34 Howard B. Be… 5 ## 10 78 ALBERT… 1240 2000 38 Gary Glenn M… 3 ## # … with 6,932 more rows, and 13 more variables: ceo_dismissal <dbl>, ## # interim_coceo <chr>, tenure_no_ceodb <dbl>, max_tenure_ceodb <dbl>, ## # fyear_gone <dbl>, leftofc <dttm>, still_there <chr>, notes <chr>, ## # sources <chr>, eight_ks <chr>, cik <dbl>, _merge <chr>, involuntary <chr>

## Bootstrapping a model

We can count up the two kinds of departures per financial year and fit the model once, for the whole data set.

library(broom) df <- departures %>% count(fyear, involuntary) %>% pivot_wider(names_from = involuntary, values_from = n) mod <- glm(cbind(involuntary, other) ~ fyear, data = df, family = "binomial") summary(mod) ## ## Call: ## glm(formula = cbind(involuntary, other) ~ fyear, family = "binomial", ## data = df) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.9858 -1.2075 -0.1947 0.7302 3.6816 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -33.236731 8.949722 -3.714 0.000204 *** ## fyear 0.015875 0.004459 3.560 0.000370 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 78.421 on 22 degrees of freedom ## Residual deviance: 65.722 on 21 degrees of freedom ## AIC: 200.86 ## ## Number of Fisher Scoring iterations: 4 tidy(mod, exponentiate = TRUE) ## # A tibble: 2 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) 3.68e-15 8.95 -3.71 0.000204 ## 2 fyear 1.02e+ 0 0.00446 3.56 0.000370

When we use `exponentiate = TRUE`

, we get the model coefficients on the linear scale instead of the logistic scale.

What we want to do is fit a model like this a whole bunch of times, instead of just once. Let’s create bootstrap resamples.

library(rsample) set.seed(123) ceo_folds <- bootstraps(departures, times = 1e3) ceo_folds ## # Bootstrap sampling ## # A tibble: 1,000 x 2 ## splits id ## <list> <chr> ## 1 <split [6942/2543]> Bootstrap0001 ## 2 <split [6942/2557]> Bootstrap0002 ## 3 <split [6942/2509]> Bootstrap0003 ## 4 <split [6942/2554]> Bootstrap0004 ## 5 <split [6942/2542]> Bootstrap0005 ## 6 <split [6942/2530]> Bootstrap0006 ## 7 <split [6942/2509]> Bootstrap0007 ## 8 <split [6942/2553]> Bootstrap0008 ## 9 <split [6942/2586]> Bootstrap0009 ## 10 <split [6942/2625]> Bootstrap0010 ## # … with 990 more rows

Now we need to make a function to count up the departures by year and type, fit our model, and return the coefficients we want.

fit_binom <- function(split) { df <- analysis(split) %>% count(fyear, involuntary) %>% pivot_wider(names_from = involuntary, values_from = n) mod <- glm(cbind(involuntary, other) ~ fyear, data = df, family = "binomial") tidy(mod, exponentiate = TRUE) }

We can apply that function to all our bootstrap resamples with `purrr::map()`

.

boot_models <- ceo_folds %>% mutate(coef_info = map(splits, fit_binom)) boot_models ## # Bootstrap sampling ## # A tibble: 1,000 x 3 ## splits id coef_info ## <list> <chr> <list> ## 1 <split [6942/2543]> Bootstrap0001 <tibble [2 × 5]> ## 2 <split [6942/2557]> Bootstrap0002 <tibble [2 × 5]> ## 3 <split [6942/2509]> Bootstrap0003 <tibble [2 × 5]> ## 4 <split [6942/2554]> Bootstrap0004 <tibble [2 × 5]> ## 5 <split [6942/2542]> Bootstrap0005 <tibble [2 × 5]> ## 6 <split [6942/2530]> Bootstrap0006 <tibble [2 × 5]> ## 7 <split [6942/2509]> Bootstrap0007 <tibble [2 × 5]> ## 8 <split [6942/2553]> Bootstrap0008 <tibble [2 × 5]> ## 9 <split [6942/2586]> Bootstrap0009 <tibble [2 × 5]> ## 10 <split [6942/2625]> Bootstrap0010 <tibble [2 × 5]> ## # … with 990 more rows

## Explore results

What did we find? We can compute bootstrap confidence intervals with `int_pctl()`

.

percentile_intervals <- int_pctl(boot_models, coef_info) percentile_intervals ## # A tibble: 2 x 6 ## term .lower .estimate .upper .alpha .method ## <chr> <dbl> <dbl> <dbl> <dbl> <chr> ## 1 (Intercept) 6.03e-23 0.0000273 0.000000246 0.05 percentile ## 2 fyear 1.01e+ 0 1.02 1.03 0.05 percentile

We can also visualize the results as well.

boot_models %>% unnest(coef_info) %>% filter(term == "fyear") %>% ggplot(aes(estimate)) + geom_vline(xintercept = 1, lty = 2, color = "gray50", size = 2) + geom_histogram() + labs( x = "Annual increase in involuntary CEO departures", title = "Over this time period, CEO departures are increasingly involuntary", subtitle = "Each passing year corresponds to a departure being 1-2% more likely to be involuntary" )

**leave a comment**for the author, please follow the link and comment on their blog:

**rstats | Julia Silge**.

R-bloggers.com offers

**daily e-mail 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/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.