Project: Grain yields

We’ll be working with grain yield data from the United States Department of Agriculture, National Agricultural Statistics Service. Unfortunately, they report all areas in acres. So, the first thing you need to do is write some utility functions to convert areas in acres to areas in hectares.

To solve this exercise, you need to know the following:

• There are 4840 square yards in an acre.
• There are 36 inches in a yard and one inch is 0.0254 meters.
• There are 10000 square meters in a hectare.

Converting areas to metric

Acres to squared yards

# Write a function to convert acres to sq. yards
acres_to_sq_yards <- function(acres) {
  acres * 4840
}

yard to meters

# Write a function to convert yards to meters
yards_to_meters <- function(yards) {
  (yards * 36) * 0.0254
}

meters to hectares

sq_meters_to_hectares <- function(sq_meters) {
  sq_meters / 10000
}

We can bring everything together to write the overall acres-to-hectares conversion function. magrittr’s raise_to_power() will be useful here.

library(magrittr)

Squared yards to squared meters

# Write a function to convert sq. yards to sq. meters
sq_yards_to_sq_meters <- function(sq_yards) {
  sq_yards %>%
    # Take the square root
    sqrt() %>%
    # Convert yards to meters
    yards_to_meters() %>%
    # Square it
    raise_to_power(2)
}

Acres to hectares

# Write a function to convert acres to hectares
acres_to_hectares <- function(acres) {
  acres %>%
    # Convert acres to sq yards
    acres_to_sq_yards() %>%
    # Convert sq yards to sq meters
    sq_yards_to_sq_meters() %>%
    # Convert sq meters to hectares
    sq_meters_to_hectares()
}

harmonic acres to hectares

Function to harmonically convert areas in acres to hectares. The function should get the reciprocal of the input, then convert from acres to hectares, then get the reciprocal again.

# reciprocal function
get_reciprocal <- function(x) {
 1/x
}

# Define a harmonic acres to hectares function
harmonic_acres_to_hectares <- function(acres) {
  acres %>% 
    # Get the reciprocal
    get_reciprocal() %>%
    # Convert acres to hectares
    acres_to_hectares() %>% 
    # Get the reciprocal again
    get_reciprocal()
}

Converting yields to metric

The yields in the NASS corn data are also given in US units, namely bushels per acre. You’ll need to write some more utility functions to convert this unit to the metric unit of kg per hectare.

Bushels historically meant a volume of 8 gallons, but in the context of grain, they are now defined as masses. This mass differs for each grain!

• One pound (lb) is 0.45359237 kilograms (kg).
• One bushel is 48 lbs of barley, 56 lbs of corn, or 60 lbs of wheat.

# Write a function to convert lb to kg
lbs_to_kgs <- function(lbs){
  lbs * 0.45359237
}

# Write a function to convert bushels to lbs
bushels_to_lbs <- function(bushels, crop) {
  # Define a lookup table of scale factors
  c(barley = 48, corn = 56, wheat = 60) %>%
    # Extract the value for the crop
    extract(crop) %>%
    # Multiply by the no. of bushels
    multiply_by(bushels)
}

# Write a function to convert bushels to kg
bushels_to_kgs <- function(bushels, crop) {
  bushels %>%
    # Convert bushels to lbs for this crop
    bushels_to_lbs(crop) %>%
    # Convert lbs to kgs
    lbs_to_kgs()
}

# Write a function to convert bushels/acre to kg/ha
bushels_per_acre_to_kgs_per_hectare <- function(bushels_per_acre, crop = c("barley", "corn", "wheat")) {
  # Match the crop argument
  crop <- match.arg(crop)
  bushels_per_acre %>%
    # Convert bushels to kgs for this crop
    bushels_to_kgs(crop) %>%
    # Convert harmonic acres to ha
    harmonic_acres_to_hectares()
}

Applying the unit conversion

Now that we’ve written some functions, it’s time to apply them! The NASS corn dataset is available, and you can fortify it (jargon for “adding new columns”) with metrics areas and yields.

This fortification process can also be turned in to a function, so you’ll define a function for this, and test it on the NASS wheat dataset.

library(dplyr)
library(agridat)
corn <- nass.corn
wheat <- nass.wheat
barley <- nass.barley
glimpse(corn)

Rows: 6,381
Columns: 4
$ year  <int> 1866, 1866, 1866, 1866, 1866, 1866, 1866, 1866, 1866, 18…
$ state <fct> Alabama, Arkansas, California, Connecticut, Delaware, Fl…
$ acres <dbl> 1050000, 280000, 42000, 57000, 200000, 125000, 1770000, …
$ yield <dbl> 9.0, 18.0, 28.0, 34.0, 23.0, 9.0, 6.0, 29.0, 36.5, 32.0,…

corn <- corn %>%
  # Add some columns
  mutate(
    # Convert farmed area from acres to ha
    farmed_area_ha = acres_to_hectares(acres),
    # Convert yield from bushels/acre to kg/ha
    yield_kg_per_ha = bushels_per_acre_to_kgs_per_hectare(
      yield,
      crop = "corn"
    )
  )

head(corn)

We can wrap the mutation code into a function, fortify_with_metric_units.

# Wrap this code into a function
fortify_with_metric_units <- function(data, crop) {
  data %>%
    mutate(
      farmed_area_ha = acres_to_hectares(acres),
      yield_kg_per_ha = bushels_per_acre_to_kgs_per_hectare(
        yield, 
        crop = crop
      )
    )
}

# Try it on the wheat dataset
wheat <- fortify_with_metric_units(wheat, crop = "wheat")
head(wheat)

Plotting yields over time

Now that the units have been dealt with, it’s time to explore the datasets. An obvious question to ask about each crop is, “how do the yields change over time in each US state?” Let’s draw a line plot to find out.

library(ggplot2)

# Using corn, plot yield (kg/ha) vs. year
ggplot(corn, aes(x = year, y = yield_kg_per_ha)) +
  # Add a line layer, grouped by state
  geom_line(aes(group = state)) +
  # Add a smooth trend layer
  geom_smooth()

Create plotting function

# Wrap this plotting code into a function
plot_yield_vs_year <- function(data){
  ggplot(data, aes(year, yield_kg_per_ha)) +
    geom_line(aes(group = state)) +
    geom_smooth()
}

We can now plot wheat

# Test it on the wheat dataset
plot_yield_vs_year(wheat)

Look at the huge increase in yields from the time of the Green Revolution in the 1950s.

A nation divided

The USA has a varied climate, so we might expect yields to differ between states. Rather than trying to reason about 50 states separately, we can use the USA Census Regions to get 9 groups.

The “Corn Belt”, where most US corn is grown is in the “West North Central” and “East North Central” regions. The “Wheat Belt” is in the “West South Central” region.

Let´s load the USA census region

usa_census_regions <- read_csv("us_data_census.csv")

Parsed with column specification:
cols(
  State = col_character(),
  `State Code` = col_character(),
  Region = col_character(),
  Division = col_character()
)

usa_census_regions <- usa_census_regions[ , -c(2, 3)]
colnames(usa_census_regions) <- c("state", "census_region")
usa_census_regions <- usa_census_regions[ , c(2, 1)]
head(usa_census_regions)

#creating census regions variable
census_regions <- unique(usa_census_regions[,1])
#("census_region" = region)
census_regions

# Inner join the corn dataset to usa_census_regions by state
corn <- corn %>%
  inner_join(usa_census_regions, by = "state")

Column `state` joining factor and character vector, coercing into character vector

head(corn)

We can turn the code into a function

fortify_with_census_region <- function(data){
  data %>%
    inner_join(usa_census_regions, by = "state")
}

wheat <- fortify_with_census_region(wheat)

Column `state` joining factor and character vector, coercing into character vector

head(wheat)

With the census data incorporated into the crop datasets, we can now look at yield differences between the regions.

Plotting yields over time by region

So far, you have a function to plot yields over time for each crop, and you’ve added a census_region column to the crop datasets. Now we are ready to look at how the yields change over time in each region of the USA.

# Plot yield vs. year for the corn dataset
plot_yield_vs_year(corn) +
  # Facet, wrapped by census region
  facet_wrap(vars(census_region))

# Wrap this code into a function
plot_yield_vs_year_by_region <- function(data){
  plot_yield_vs_year(data) +
    facet_wrap(vars(census_region))
}

# Try it on the wheat dataset
plot_yield_vs_year_by_region(wheat)

Radical regional yield analysis! Reassuringly, the corn yields are highest in the West North Central region, the heart of the Corn Belt. For wheat, it looks like the yields are highest in the Wheat Belt (West South Central region) have been overtaken by some other regions.

Running a model

The smooth trend line we see in the plots of yield over time use a generalized additive model (GAM) to determine where the line should lie. This sort of model is ideal for fitting nonlinear curves. So we can make predictions about future yields, let’s explicitly run the model. The syntax for running this GAM takes the following form.

gam(response ~ s(explanatory_var1) + explanatory_var2, data = dataset)

Here, s() means “make the variable smooth”, where smooth very roughly means nonlinear.

library(mgcv)
# Run a generalized additive model of 
# yield vs. smoothed year and census region
gam(yield_kg_per_ha ~ s(year) + census_region, data = corn)

Family: gaussian 
Link function: identity 

Formula:
yield_kg_per_ha ~ s(year) + census_region

Estimated degrees of freedom:
7.64  total = 16.64 

GCV score: 836963     

We can wrap the code into a function and apply to wheat as well.

# Wrap the model code into a function
run_gam_yield_vs_year_by_region <- function(data){
  gam(yield_kg_per_ha ~ s(year) + census_region, data = data)
}

# Try it on the wheat dataset
run_gam_yield_vs_year_by_region(wheat)

Family: gaussian 
Link function: identity 

Formula:
yield_kg_per_ha ~ s(year) + census_region

Estimated degrees of freedom:
7.03  total = 16.03 

GCV score: 316685.5     

Making yield predictions

We can make predictions using a call to predict(), in the following form.

predict(model, cases_to_predict, type = "response")

# we can assign our model to a variable
corn_model = run_gam_yield_vs_year_by_region(corn)
wheat_model = run_gam_yield_vs_year_by_region(wheat)

# Make predictions in 2050  
predict_this <- data.frame(
  year = 2050,
  census_region = census_regions
)
predict_this

# Predict the yield
pred_yield_kg_per_ha <- predict(corn_model,
                                predict_this,
                                type = "response")

corn_prediction <- predict_this %>%
  # Add the prediction as a column of predict_this 
  mutate(pred_yield_kg_per_ha = c(pred_yield_kg_per_ha))
corn_prediction

We can wrap the script into a function

# Wrap this prediction code into a function
predict_yields <- function(model, year){
  predict_this <- data.frame(
    year = year,
    census_region = census_regions
  ) 
  pred_yield_kg_per_ha <- predict(model, predict_this, type = "response")
  predict_this %>%
    mutate(pred_yield_kg_per_ha = c(pred_yield_kg_per_ha))
}

# Try it on the wheat dataset
wheat_prediction <- predict_yields(wheat_model,
               year = 2050)
wheat_prediction

The models predict that in 2050, the highest yields will be in the Pacific region for both corn and wheat.

Do it all over again

Now you are going to rerun the whole analysis from this chapter on a new crop, barley. Since all the infrastructure is in place, that’s less effort than it sounds!

Barley prefers a cooler climate compared to corn and wheat and is commonly grown in the US mountain states of Idaho and Montana.

fortified_barley <- barley %>% 
  # Fortify with metric units
  fortify_with_metric_units("barley") %>%
  # Fortify with census regions
  fortify_with_census_region()

Column `state` joining factor and character vector, coercing into character vector

# See the result
glimpse(fortified_barley)

Rows: 4,839
Columns: 7
$ year            <int> 1866, 1866, 1866, 1866, 1866, 1866, 1866, 1866…
$ state           <chr> "Connecticut", "Illinois", "Indiana", "Iowa", …
$ acres           <dbl> 1000, 96000, 11000, 66000, 2000, 10000, 34000,…
$ yield           <dbl> 22.5, 23.4, 23.0, 22.0, 23.0, 23.5, 21.5, 25.5…
$ farmed_area_ha  <dbl> 404.6856, 38849.8217, 4451.5421, 26709.2524, 8…
$ yield_kg_per_ha <dbl> 1210.5192, 1258.9400, 1237.4197, 1183.6188, 12…
$ census_region   <chr> "New England", "East North Central", "East Nor…

# Plot yield vs. year by region
plot_yield_vs_year_by_region(fortified_barley)

fortified_barley %>% 
  # Run a GAM of yield vs. year by region
  run_gam_yield_vs_year_by_region  %>% 
  # Make predictions of yields in 2050
  predict_yields(2050)

Since all your analysis code was contained in functions, it was really simple to apply it to another dataset. Here you can see that yields are highest in the Mountain region, and the model predicts that this will still be the case in 2050.