Train and analyze many models for #TidyTuesday crop yields

Lately I’ve been publishing screencasts demonstrating how to use the tidymodels framework, from just getting started to tuning more complex models. Today’s screencast explores how to fluently apply tidy data principles to the task of building many models using with this week’s #TidyTuesday dataset on crop yields. 🌽

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

Explore the data

Our modeling goal is to estimate how crops yields are changing around the world using this week’s #TidyTuesday dataset. We can build many models for the country-crop combinations we are interested in.

First, let’s read in two of the datasets for this week.

library(tidyverse)

key_crop_yields <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-09-01/key_crop_yields.csv")
land_use <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-09-01/land_use_vs_yield_change_in_cereal_production.csv")

I’m going to use the land_use dataset only to find the top population countries. Let’s create a vector of their names.

top_countries <- land_use %>%
  janitor::clean_names() %>%
  filter(!is.na(code), entity != "World") %>%
  group_by(entity) %>%
  filter(year == max(year)) %>%
  ungroup() %>%
  slice_max(total_population_gapminder, n = 30) %>%
  pull(entity)

top_countries

##  [1] "China"                        "India"                       
##  [3] "United States"                "Indonesia"                   
##  [5] "Pakistan"                     "Brazil"                      
##  [7] "Nigeria"                      "Bangladesh"                  
##  [9] "Russia"                       "Mexico"                      
## [11] "Japan"                        "Ethiopia"                    
## [13] "Philippines"                  "Egypt"                       
## [15] "Vietnam"                      "Democratic Republic of Congo"
## [17] "Germany"                      "Turkey"                      
## [19] "Iran"                         "Thailand"                    
## [21] "United Kingdom"               "France"                      
## [23] "Italy"                        "South Africa"                
## [25] "Tanzania"                     "Myanmar"                     
## [27] "Kenya"                        "South Korea"                 
## [29] "Colombia"                     "Spain"

Now let’s create a tidy version of the crop yields data, for the countries and crops I am interested in.

tidy_yields <- key_crop_yields %>%
  janitor::clean_names() %>%
  pivot_longer(wheat_tonnes_per_hectare:bananas_tonnes_per_hectare,
    names_to = "crop", values_to = "yield"
  ) %>%
  mutate(crop = str_remove(crop, "_tonnes_per_hectare")) %>%
  filter(
    crop %in% c("wheat", "rice", "maize", "barley"),
    entity %in% top_countries,
    !is.na(yield)
  )

tidy_yields

## # A tibble: 6,032 x 5
##    entity     code   year crop   yield
##    <chr>      <chr> <dbl> <chr>  <dbl>
##  1 Bangladesh BGD    1961 wheat  0.574
##  2 Bangladesh BGD    1961 rice   1.70 
##  3 Bangladesh BGD    1961 maize  0.799
##  4 Bangladesh BGD    1961 barley 0.577
##  5 Bangladesh BGD    1962 wheat  0.675
##  6 Bangladesh BGD    1962 rice   1.53 
##  7 Bangladesh BGD    1962 maize  0.738
##  8 Bangladesh BGD    1962 barley 0.544
##  9 Bangladesh BGD    1963 wheat  0.607
## 10 Bangladesh BGD    1963 rice   1.77 
## # … with 6,022 more rows

This data structure is just right for plotting crop yield over time!

tidy_yields %>%
  ggplot(aes(year, yield, color = crop)) +
  geom_line(alpha = 0.7, size = 1.5) +
  geom_point() +
  facet_wrap(~entity, ncol = 5) +
  scale_x_continuous(guide = guide_axis(angle = 90)) +
  labs(x = NULL, y = "yield (tons per hectare)")

Notice that not all countries produce all crops, but that overall crop yields are increasing.

Many models

Now let’s fit a linear model to each country-crop combination.

library(tidymodels)

tidy_lm <- tidy_yields %>%
  nest(yields = c(year, yield)) %>%
  mutate(model = map(yields, ~ lm(yield ~ year, data = .x)))

tidy_lm

## # A tibble: 111 x 5
##    entity     code  crop   yields            model 
##    <chr>      <chr> <chr>  <list>            <list>
##  1 Bangladesh BGD   wheat  <tibble [58 × 2]> <lm>  
##  2 Bangladesh BGD   rice   <tibble [58 × 2]> <lm>  
##  3 Bangladesh BGD   maize  <tibble [58 × 2]> <lm>  
##  4 Bangladesh BGD   barley <tibble [58 × 2]> <lm>  
##  5 Brazil     BRA   wheat  <tibble [58 × 2]> <lm>  
##  6 Brazil     BRA   rice   <tibble [58 × 2]> <lm>  
##  7 Brazil     BRA   maize  <tibble [58 × 2]> <lm>  
##  8 Brazil     BRA   barley <tibble [58 × 2]> <lm>  
##  9 China      CHN   wheat  <tibble [58 × 2]> <lm>  
## 10 China      CHN   rice   <tibble [58 × 2]> <lm>  
## # … with 101 more rows

Next, let’s tidy() those models to get out the coefficients, and adjust the p-values for multiple comparisons while we’re at it.

slopes <- tidy_lm %>%
  mutate(coefs = map(model, tidy)) %>%
  unnest(coefs) %>%
  filter(term == "year") %>%
  mutate(p.value = p.adjust(p.value))

slopes

## # A tibble: 111 x 10
##    entity  code  crop  yields  model term  estimate std.error statistic  p.value
##    <chr>   <chr> <chr> <list>  <lis> <chr>    <dbl>     <dbl>     <dbl>    <dbl>
##  1 Bangla… BGD   wheat <tibbl… <lm>  year   0.0389   0.00253      15.4  5.11e-20
##  2 Bangla… BGD   rice  <tibbl… <lm>  year   0.0600   0.00231      26.0  6.05e-31
##  3 Bangla… BGD   maize <tibbl… <lm>  year   0.122    0.0107       11.3  1.82e-14
##  4 Bangla… BGD   barl… <tibbl… <lm>  year   0.00505  0.000596      8.47 4.34e-10
##  5 Brazil  BRA   wheat <tibbl… <lm>  year   0.0366   0.00222      16.5  2.55e-21
##  6 Brazil  BRA   rice  <tibbl… <lm>  year   0.0755   0.00490      15.4  4.96e-20
##  7 Brazil  BRA   maize <tibbl… <lm>  year   0.0709   0.00395      18.0  4.37e-23
##  8 Brazil  BRA   barl… <tibbl… <lm>  year   0.0466   0.00319      14.6  5.05e-19
##  9 China   CHN   wheat <tibbl… <lm>  year   0.0880   0.00141      62.6  1.72e-51
## 10 China   CHN   rice  <tibbl… <lm>  year   0.0843   0.00289      29.2  1.47e-33
## # … with 101 more rows

Explore results

Now we can visualize the results of this modeling, which is estimating how crop yields are changing around the world.

library(ggrepel)
slopes %>%
  ggplot(aes(estimate, p.value, label = entity)) +
  geom_vline(
    xintercept = 0, lty = 2,
    size = 1.5, alpha = 0.7, color = "gray50"
  ) +
  geom_point(aes(color = crop), alpha = 0.8, size = 2.5, show.legend = FALSE) +
  scale_y_log10() +
  facet_wrap(~crop) +
  geom_text_repel(size = 3, family = "IBMPlexSans") +
  theme_light(base_family = "IBMPlexSans") +
  theme(strip.text = element_text(family = "IBMPlexSans-Bold", size = 12)) +
  labs(x = "increase in tons per hectare per year")

• On the x-axis is the slope of these models. Notice that most countries are on the positive side, with increasing crop yields. The further to the right a country is, the larger the increase in crop yield over this time period. Corn yields have increased the most.

• On the y-axis is the p-value, a measure of how surprising the effect we see is under the assumption of no relationship (no change with time). Countries lower in the plots have smaller p-values; we are more certain those are real relationships.

We can extend this to check out how well these models fit the data with glance(). This approach for using statistical models to estimate changes in many subgroups at once has been so helpful to me in many situations!