Fit the model for more countries

I’d like to get predictions for more countries so we can compare how the USA is doing with, say, Italy or South Korea. I’m going back to reported case data from JHU for this, and I’m going to use the Negative Binomial model because then I don’t have to worry about the odd zero close to the start of the time series. I’ll process the other countries separately because of the funky county/state business in the USA, then stick them together afterwards.

countries <- c("Canada", "Australia", "Korea, South", "Italy", "Iran")
other_confirmed_total <- jhu_wide %>% 
  rename(province = "Province/State", 
         country_region = "Country/Region") %>% 
  pivot_longer(-c(province, country_region, Lat, Long), 
               names_to = "Date", values_to = "cumulative_cases") %>% 
  filter(country_region %in% countries,
         # have to trap the rows with missing province (most other countries)
         # otherwise str_detect(province ...) is missing and dropped by filter()
         is.na(province) | str_detect(province, "Princess", negate = TRUE)) %>% 
  mutate(Date = mdy(Date)) %>% 
  # filter out state rows prior to march 9, and county rows after that. 
  group_by(country_region, Date) %>% # then group by country and Date and sum
  summarize(cumulative_cases = sum(cumulative_cases)) %>% 
  group_by(country_region) %>% 
  mutate(incident_cases = c(0, diff(cumulative_cases))) %>% 
  ungroup()
us_confirmed_total <- mutate(us_confirmed_total, country_region = "USA")
all_confirmed_total <- bind_rows(other_confirmed_total, us_confirmed_total)

pall <- ggplot(data = all_confirmed_total,
       mapping = aes(x = Date)) + # don't add y = here, so we can change variables later for ribbons etc
  geom_point(mapping = aes(y = incident_cases)) + 
  facet_wrap(~country_region, scales = "free_y")
pall

Phew. That took far too long to sort out, because I tried to filter() all the countries including the USA, and that was just not working. Treating the USA separately was easy because I already had that dataframe built above. The y axis scales are wildly different on those, but just eyeballing the shape it’s clear that South Korea is doing something radically different. However, the exponential growth in South Korea started earlier too, so that trajectory might not be impossible. There are also a couple weird outliers in Australia and Italy where incident cases dropped to zero suddenly for one day. Canada also has a strange 2 day oscillation going on. I think these are likely reporting artifacts.

Now fit the models and see what that looks like. The early period of the time series for all these countries is mostly driven by travel related cases, so as before I’ll trim off the data prior to Feb 15. I picked that date by eyeballing the data and looking for a date prior to where case loads started accelerating for all countries. It would be better to have that set uniquely for each country, and I’ve got some thoughts on how to do that, but I’ll leave that for tomorrow.

all_results <- all_confirmed_total %>%
  mutate(day = as.numeric(Date - ymd("2020-02-27"))) %>% # get day relative to Feb 27
  filter(Date > "2020-02-15", Date <= "2020-03-11") %>% # trim data to the dates we want
  group_by(country_region) %>% 
  # nest() puts the data for each country into a "list column" called "data" 
  # so that I can do the next step
  nest() %>% 
  # fit the models to each dataframe in the list column
  mutate(models = map(data, ~gam(incident_cases ~ day, data = ., family = nb)),
         data = map2(models, data, ~bind_cols(.y, tibble(.fitted = fitted.values(.x))))) ## adds regression output to data 
all_fitted <- all_results %>% 
  select(-models) %>%  # remove the models column before the next step
  unnest(cols = data)  # pulls the "data" column back out to make a single data frame
pall + 
  geom_line(data = all_fitted,
            mapping = aes(y = .fitted))

That’s a first cut, because broom functions don’t work on mgcv::gam objects. So I need to write a custom function to get what I want. But it certainly looks like all countries are showing less growth in cases than expected from a simple exponential growth model. That’s good news.

augment_gam <- function(x, newdata = NULL, se.fit = FALSE){
  if (is.null(newdata)){
    result <- predict(x, se.fit = se.fit)
    original_data <- x$model
  } else {
    result <- predict(x, newdata = newdata, se.fit = se.fit)
    original_data <- newdata
  }
  if (se.fit){
    result <- bind_cols(original_data,
                     tibble(.fitted = result$fit,
                     .se.fit = result$se.fit))
  } else {
    result <- bind_cols(original_data,
                     tibble(.fitted = result))
  }
  return(result)
}

all_results <- all_confirmed_total %>%
  mutate(day = as.numeric(Date - ymd("2020-02-27"))) %>% # get day relative to Feb 27
  filter(Date > "2020-02-15", Date <= "2020-03-11") %>% # trim data to the dates we want
  group_by(country_region) %>% 
  # nest() puts the data for each country into a "list column" called "data" 
  # so that I can do the next step
  nest() %>% 
  # fit the models to each dataframe in the list column
  mutate(models = map(data, ~gam(incident_cases ~ day, data = ., family = nb)),
         data = map(models, augment_gam, se.fit = TRUE), ## adds regression output to data
         theta = map(models, ~.x$family$getTheta(TRUE)),
         df = map_dbl(models, "df.residual"))  

all_fitted <- all_results %>% 
  select(-models) %>%  # remove the models column before the next step
  unnest(cols = data) %>%   # pulls the "data" column back out to make a single data frame
  mutate(Date = ymd("2020-02-27")+days(day), # recreate the Date column
         fitted_cases = exp(.fitted),        # .fitted on log scale
         lcl_cases = exp(.fitted - qt(0.975, df)*.se.fit),
         ucl_cases = exp(.fitted + qt(0.975, df)*.se.fit))

predict_data <- tibble(day = 14:20)
predicted <- all_results %>% 
  mutate(predictions = map(models, augment_gam, se.fit = TRUE, newdata = predict_data),
         theta = map(models, ~.x$family$getTheta(TRUE)),
         df = map_dbl(models, "df.residual")) %>% 
  select(-models, -data) %>% 
  unnest(cols = c("predictions", "theta")) %>% 
  mutate(Date = ymd("2020-02-27")+days(day), # recreate the Date column
         expected_cases = exp(.fitted),
         lpl = qnbinom(0.025, mu = expected_cases, size = theta),
         upl = qnbinom(0.975, mu = expected_cases, size = theta))
pall2 <- ggplot(data = filter(all_confirmed_total, Date >= ymd("2020-02-15")),
       mapping = aes(x = Date)) + # don't add y = here, so we can change variables later for ribbons etc
  geom_point(mapping = aes(y = incident_cases)) + 
  facet_wrap(~country_region, scales = "free_y")

pall2 + 
  geom_line(data = all_fitted,
            mapping = aes(y = fitted_cases), 
            col = "blue", size = 1.5) + 
  geom_ribbon(data = all_fitted,
              mapping = aes(ymin = lcl_cases, ymax = ucl_cases),
              alpha = 0.25) + 
  geom_line(data = predicted,
            mapping = aes(y = expected_cases),
            linetype = 2) + 
  geom_ribbon(data = predicted,
              mapping = aes(ymin = lpl, ymax = upl),
              alpha = 0.25) +
  scale_y_log10()

Clearly I need to figure out country specific ranges to fit the model. South Korea is clearly not growing anymore, so the terrible model fit is not unexpected. Iran is either doing great, or lying about the number of new cases, opinion out there seems divided on this point.

One excellent point raised by Brian McGill in the open discussion on Dyanmic Ecology is that the extreme response of South Korea and other countries which works in the short term, might not be the best strategy in the long term. The “squash it” strategy works for a disease like SARS or MERS which are less contagious, and can therefore really be put to bed for good. But seems like COVID-19 has exceeded global capacity to squash it everywhere. Which means we need to figure out how to live with it long term. With an estimated R0 = 2.5 COVID-19 will continue to be an issue until ~ 60% of the population are immune. I think the real goal of the “stretch it out” strategy by social distancing etc. is to buy us enough time to get a vaccine and build herd immunity artificially.

Next steps

I want to add some more models to the mix, I’m thinking maybe a logistic curve for the cumulative number of cases as suggested by my friend Scott Bradshaw (an MD), and then maybe a proper epidemiological model.