Analyzing US Dairy Consumption with Time Series Prediction (sweep package used)

Milk products facts

milk_pivot <- milk_products_facts %>%
  pivot_longer(-year, names_to = "product", values_to = "consumption")

milk_pivot %>%
  mutate(product = str_replace_all(product, "_", " ")) %>%
  ggplot(aes(year, consumption, color = product)) +
  geom_line(aes(group = 1), size = 2, show.legend = FALSE) +
  facet_wrap(~product, scales = "free_y") +
  expand_limits(y = 0) +
  theme(
    strip.text = element_text(size = 13, face = "bold")
  )  +
  labs(y = "lbs per person",
       title = "Milk Product Consumption (lbs per person)") 

Time series prediction

This section is inspired by David Robinson’s code, and this is also the first time of me using sweep and timetk packages.

milk_productes_tidied <- milk_pivot %>%
  separate(product, c("category", "product"), sep = "_",
           extra = "merge", fill = "right")%>%
  mutate(
    product = coalesce(product, category),
    product = str_to_title(str_replace_all(product, "_", " ")),
    category = str_to_title(category)
  )

milk_ts <- milk_productes_tidied %>%
  mutate(year = make_date(year)) %>% 
  nest(-c(category, product)) %>%
  mutate(ts = map(data, tk_ts, start = 1975))


milk_ets <- milk_ts %>% 
  mutate(model = map(ts, ets))

# the following commented code does not work due to tidyr update
# milk_ets %>%
#   unnest(map(model, sw_glance))

# This is the correct way to do it
milk_ets %>%
  mutate(glanced = map(model, sw_glance)) %>%
  unnest(glanced)

## # A tibble: 17 x 17
##    category product  data   ts    model model.desc  sigma  logLik    AIC     BIC
##    <chr>    <chr>    <list> <lis> <lis> <chr>       <dbl>   <dbl>  <dbl>   <dbl>
##  1 Fluid    Milk     <tibb~ <ts ~ <ets> ETS(A,A,N) 2.05   -110.    229.   238.  
##  2 Fluid    Yogurt   <tibb~ <ts ~ <ets> ETS(M,A,N) 0.0694  -40.3    90.7   99.5 
##  3 Butter   Butter   <tibb~ <ts ~ <ets> ETS(A,N,N) 0.192    -8.92   23.8   29.1 
##  4 Cheese   American <tibb~ <ts ~ <ets> ETS(A,Ad,~ 0.352   -33.3    78.7   89.3 
##  5 Cheese   Other    <tibb~ <ts ~ <ets> ETS(M,A,N) 0.0214  -26.4    62.7   71.5 
##  6 Cheese   Cottage  <tibb~ <ts ~ <ets> ETS(A,Ad,~ 0.0799   30.5   -48.9  -38.4 
##  7 Evap     Cnd Can~ <tibb~ <ts ~ <ets> ETS(A,Ad,~ 0.186    -5.91   23.8   34.4 
##  8 Evap     Cnd Bul~ <tibb~ <ts ~ <ets> ETS(A,N,N) 0.128     8.42  -10.8   -5.56
##  9 Evap     Cnd Bul~ <tibb~ <ts ~ <ets> ETS(A,N,N) 0.378   -38.0    81.9   87.2 
## 10 Frozen   Ice Cre~ <tibb~ <ts ~ <ets> ETS(M,N,N) 0.0286  -45.4    96.7  102.  
## 11 Frozen   Ice Cre~ <tibb~ <ts ~ <ets> ETS(M,N,N) 0.0586  -37.6    81.1   86.4 
## 12 Frozen   Sherbet  <tibb~ <ts ~ <ets> ETS(A,N,N) 0.0510   48.1   -90.2  -84.9 
## 13 Frozen   Other    <tibb~ <ts ~ <ets> ETS(M,N,N) 0.172   -49.0   104.   109.  
## 14 Dry      Whole M~ <tibb~ <ts ~ <ets> ETS(A,N,N) 0.0782   29.8   -53.5  -48.2 
## 15 Dry      Nonfat ~ <tibb~ <ts ~ <ets> ETS(M,N,N) 0.148   -44.9    95.7  101.  
## 16 Dry      Butterm~ <tibb~ <ts ~ <ets> ETS(M,N,N) 0.141    69.7  -133.  -128.  
## 17 Dry      Whey     <tibb~ <ts ~ <ets> ETS(A,N,N) 0.260   -22.0    50.0   55.3 
## # ... with 7 more variables: ME <dbl>, RMSE <dbl>, MAE <dbl>, MPE <dbl>,
## #   MAPE <dbl>, MASE <dbl>, ACF1 <dbl>

milk_ts %>%
  crossing(model_name = c("auto.arima", "ets")) %>%
  mutate(model = map2(model_name, ts, ~ invoke(.x, list(.y))),
         forecast = map(model, forecast, h = 10)) %>%
  mutate(sweeped = map(forecast, sw_sweep)) %>%
  unnest(sweeped)  %>%
  ggplot(aes(index, consumption, color = model_name, lty = key)) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = lo.95, ymax = hi.95, fill = model_name), alpha = .5) +
  facet_wrap(~ product, scales = "free_y") +
  theme(strip.text = element_text(size = 13, face = "bold")) +
  expand_limits(y = 0) +
  scale_x_continuous(breaks = c(1980, 2000, 2020)) +
  scale_linetype_discrete(guide = FALSE) +
  labs(x = "Year",
       y = "Average US consumption (lbs per person)",
       title = "Forecasted consumption of dairy products",
       subtitle = "Based on USDA data 1975-2017. Showing 95% prediction intervals.",
       color = "Model")

Cheese

cheese %>%
  rename(year = Year) %>%
  pivot_longer(-year, names_to = "product", values_to = "consumption") %>% 
  ggplot(aes(year, consumption, color = product)) +
  geom_line(size = 2, show.legend = FALSE) +
  facet_wrap(~product, scales = "free_y") +
  expand_limits(y = 0) +
  theme(
    strip.text = element_text(size = 15, face = "bold")
  )  +
  labs(y = "lbs per person",
       title = "Cheese Product Consumption (lbs per person)") 

Milk sales

milk_sales %>%
  ggplot(aes(year, pounds, color = milk_type)) +
  geom_line(show.legend = FALSE, size = 2) +
  facet_wrap(~milk_type) +
  expand_limits(y = 0) +
  theme(
    strip.text = element_text(size = 15, face = "bold")
  )  +
  labs(title = "Fluid Milk Sales (Pounds)") 

milk_sales %>%
  ggplot(aes(year, pounds, color = milk_type)) +
  geom_line(show.legend = FALSE, size = 2) +
  facet_wrap(~milk_type, scales = "free") +
  expand_limits(y = 0) +
  theme(
    strip.text = element_text(size = 15, face = "bold")
  )  +
  labs(title = "Free-scale Fluid Milk Sales (Pounds)") 

Cows facts

cows %>%
  ggplot(aes(year, milk_per_cow)) +
  geom_line() +
  expand_limits(y = 0)

cows %>%
  ggplot(aes(year, avg_milk_cow_number)) +
  geom_line() +
  expand_limits(y = 0)

Each cow produced more milk annually from 1980 to 2010, although the # of cows decreased in the same time interval.

Milk production in each state

state_milk %>%
  ggplot(aes(year, state, size = milk_produced, color = state)) +
  geom_point() +
  scale_color_discrete(guide = FALSE) +
  labs(y = NULL)

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

If you think the above plot is too fuzzy, the following heatmap is created to convey the same information.

state_milk %>%
  ggplot(aes(factor(year), state, fill = milk_produced)) +
  geom_tile() +
  scale_fill_gradient2(high = "red",
                       low = "white",
                       mid = "blue",
                       midpoint = mean(state_milk$milk_produced)) +
  labs(
    x = NULL,
    y = NULL,
    fill = "pounds of milk produced",
    title = "State-wise Milk Produced (pounds)"
  ) +
  scale_x_discrete(breaks = seq(1970, 2017, by = 3))