CDSS Hackathon 2020

DK & SE — Time Series on Covid-19 Data Using Stan

Background

How has COVID-19 affected mobility of people and goods across the globe?

We are interested in the pandemic’s effects on society. With so many different and disparate approaches, it is difficult to see what measures correlated with what behavior. We decided to compare two regions, similar in geographic location, population size, culture, and history: Sweden and Denmark. Their mobility behavior and COVID-19 statistics are reflected in Google Community Mobility Reports and ECDC case/death series (see data/README.md).

We chose these regions because the approach to mitigating the effects of COVID-19 in terms of mobility and the resultant economic and health implications were considerably different. While both effectively closed their borders, mobility was severely curtailed in Denmark (lockdown) whereas Sweden only issued recommendations (non-mandatory). These are contrasting responses for geographies that are physically connected and whose populations are virtually homogeneous to each other.

An important consideration for these datasets is to treat time itself as a variable. Here we consider time in terms of before-and-after border closures, and we employ a before-and-after modeling technique called difference-in-differences (pretest-posttest with comparison group). For this to be a valid method, we assume parallel trends — reflected in our choice of Sweden and Denmark.

To incorporate uncertainty we use Bayesian modeling via Stan in R (rstanarm::stan_glm), which gives a posterior distribution for our estimates rather than single point estimates.

mobility_raw <- read_csv("data/google_mobility_swe_den.csv", show_col_types = FALSE)
covid_raw <- read_csv("data/covid19_ecdc.csv", show_col_types = FALSE)

covid19 <- covid_raw |>
  filter(countries_and_territories %in% c("Sweden", "Denmark"), year == 2020) |>
  arrange(date) |>
  transmute(
    country = countries_and_territories,
    date,
    daily_deaths,
    daily_confirmed_cases,
    deaths,
    pop = pop_data_2019,
    daily_deaths_per_mil = 1e6 * daily_deaths / pop_data_2019,
    daily_confirmed_cases_per_mil = 1e6 * daily_confirmed_cases / pop_data_2019,
    deaths_per_mil = 1e6 * deaths / pop_data_2019
  )

mobility_region <- mobility_raw |>
  filter(country_region %in% c("Sweden", "Denmark")) |>
  transmute(
    country = country_region,
    region = sub_region_1,
    date,
    retail_rec = retail_and_recreation_percent_change_from_baseline,
    grocery_pharm = grocery_and_pharmacy_percent_change_from_baseline,
    parks = parks_percent_change_from_baseline,
    transit = transit_stations_percent_change_from_baseline,
    workplace = workplaces_percent_change_from_baseline,
    residential = residential_percent_change_from_baseline
  ) |>
  mutate(
    yr_month = format(as.Date(date), "%Y-%m"),
    week = week(ymd(date))
  )

# Country-day mobility means (avoid broadcasting country COVID onto every region).
mobility_cd <- mobility_region |>
  group_by(country, date, yr_month, week) |>
  summarise(
    across(
      c(retail_rec, grocery_pharm, parks, transit, workplace, residential),
      ~ mean(.x, na.rm = TRUE)
    ),
    .groups = "drop"
  )

panel <- mobility_cd |>
  left_join(covid19, by = c("country", "date")) |>
  mutate(
    # Shared post period: after Sweden's non-EU border cutoff (DK closed Mar 14; SE Mar 19).
    closed = as.integer(date > as.Date("2020-03-19")),
    country = factor(country, levels = c("Denmark", "Sweden"))
  )

tibble(
  n_region_day = nrow(mobility_region),
  n_country_day = nrow(panel),
  date_min = min(panel$date),
  date_max = max(panel$date)
)

# A tibble: 1 × 4
  n_region_day n_country_day date_min   date_max  
         <int>         <int> <date>     <date>    
1        68976           420 2020-02-15 2020-09-11

Time series models are sensitive to missingness, so check the region-day join (before the country-day collapse used in the DiD).

mobility_vis <- mobility_region |>
  left_join(covid19, by = c("country", "date"))

# Your vis_miss on the region-day join; column order kept; labels under the plot with margin so long names are not clipped.
p_miss <- vis_miss(mobility_vis, warn_large_data = FALSE)
x_labs <- miss_var_summary(mobility_vis) |>
  mutate(lab = paste0(variable, " (", round(as.numeric(pct_miss), 1), "%)")) |>
  arrange(match(variable, names(mobility_vis))) |>
  pull(lab)
names(x_labs) <- names(mobility_vis)
p_miss +
  scale_x_discrete(limits = names(mobility_vis), labels = x_labs, position = "bottom") +
  theme(
    axis.text.x = element_text(angle = 55, hjust = 1, vjust = 1, size = 8.5),
    legend.position = "top",
    plot.margin = unit(c(0.5, 1.2, 4.2, 0.5), "cm")
  )

The data seems to have been collected sparsely for some variables on some days, so we will aggregate into weeks for the charts below, and the DiD models use the country-day panel. The jitter from day to day is likely too granular for this model anyway.

Data exploration

After border closures

Denmark closed its border on March 14, 2020. Sweden closed its borders to non-EU member states on March 19, 2020, and virtually all EU countries imposed travel restrictions at about the same time, effectively prohibiting travel into Sweden. Both are closed within a week of each other; we treat 2020-03-19 as the shared post-period cutoff.

Change in mobility

mob_week <- panel |>
  group_by(country, week) |>
  summarise(
    across(
      c(retail_rec, grocery_pharm, transit, workplace, residential),
      ~ mean(.x, na.rm = TRUE)
    ),
    .groups = "drop"
  ) |>
  pivot_longer(
    cols = -c(country, week),
    names_to = "series",
    values_to = "mean"
  ) |>
  mutate(
    series = recode(
      series,
      retail_rec = "Retail",
      grocery_pharm = "Grocery & Pharmacy",
      transit = "Transit",
      workplace = "Workplace",
      residential = "Residential"
    )
  )

ggplot(mob_week, aes(week, mean, colour = country, group = country)) +
  geom_point(size = 1.2) +
  geom_line(linewidth = 0.6) +
  geom_vline(xintercept = 12, linetype = "dotted", colour = col_close, linewidth = 0.85) +
  facet_wrap(~series, ncol = 1, scales = "free_y") +
  scale_color_manual(values = c(Denmark = col_dnk, Sweden = col_swe)) +
  labs(
    title = "Changes in mobility over time by week",
    x = "Week",
    y = "% change from baseline",
    colour = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

Change in COVID metrics

covid_week <- panel |>
  group_by(country, week) |>
  summarise(
    daily_confirmed_cases_per_mil = mean(daily_confirmed_cases_per_mil, na.rm = TRUE),
    daily_deaths_per_mil = mean(daily_deaths_per_mil, na.rm = TRUE),
    deaths_per_mil = mean(deaths_per_mil, na.rm = TRUE),
    .groups = "drop"
  ) |>
  pivot_longer(
    cols = -c(country, week),
    names_to = "series",
    values_to = "mean"
  ) |>
  mutate(
    series = recode(
      series,
      daily_confirmed_cases_per_mil = "Daily confirmed cases per million",
      daily_deaths_per_mil = "Daily deaths per million",
      deaths_per_mil = "Cumulative deaths per million"
    )
  )

ggplot(covid_week, aes(week, mean, colour = country, group = country)) +
  geom_point(size = 1.2) +
  geom_line(linewidth = 0.6) +
  geom_vline(xintercept = 12, linetype = "dotted", colour = col_close, linewidth = 0.85) +
  facet_wrap(~series, ncol = 1, scales = "free_y") +
  scale_color_manual(values = c(Denmark = col_dnk, Sweden = col_swe)) +
  labs(
    title = "Weekly averages of COVID-19 statistics",
    x = "Week",
    y = NULL,
    colour = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

Difference in differences

What is the effect of a mobility-related policy on COVID-related metrics?

A difference-in-differences model compares COVID metrics before-and-after border closure for Sweden to the before-and-after for Denmark. For two countries that can be considered quite easily the same country, we would expect the closure of borders to have comparable effects on COVID metrics under parallel trends. The DiD interaction closed:countrySweden is the quantity of interest.

In conditional mean terms

did_means <- function(y) {
  panel |>
    filter(!is.na(.data[[y]])) |>
    group_by(country, closed) |>
    summarise(mean = mean(.data[[y]], na.rm = TRUE), .groups = "drop") |>
    mutate(outcome = y)
}

means_cases <- did_means("daily_confirmed_cases_per_mil")
means_deaths <- did_means("daily_deaths_per_mil")

cell <- function(tbl, ctry, cl) {
  tbl$mean[tbl$country == ctry & tbl$closed == cl]
}

delta <- function(tbl, ctry) cell(tbl, ctry, 1L) - cell(tbl, ctry, 0L)

tibble(
  outcome = c("daily_confirmed_cases_per_mil", "daily_deaths_per_mil"),
  sweden_before_after = c(delta(means_cases, "Sweden"), delta(means_deaths, "Sweden")),
  denmark_before_after = c(delta(means_cases, "Denmark"), delta(means_deaths, "Denmark")),
  naive_did = c(
    delta(means_cases, "Sweden") - delta(means_cases, "Denmark"),
    delta(means_deaths, "Sweden") - delta(means_deaths, "Denmark")
  )
)

# A tibble: 2 × 4
  outcome                     sweden_before_after denmark_before_after naive_did
  <chr>                                     <dbl>                <dbl>     <dbl>
1 daily_confirmed_cases_per_…               43.2                11.8       31.5 
2 daily_deaths_per_mil                       3.21                0.591      2.62

Observing only differences in means of daily confirmed cases per million, Sweden’s before→after change exceeds Denmark’s by the naive DiD above (same construction for daily deaths per million).

In regression terms

Primary specification on the country-day panel:

\[ y_{ct} = \alpha + \beta_1\,\mathrm{closed}_t + \beta_2\,\mathrm{Sweden}_c + \beta_3\,(\mathrm{closed}_t\times\mathrm{Sweden}_c) + \varepsilon_{ct}. \]

\(\beta_3\) is the DiD estimate. A second model adds retail and grocery/pharmacy mobility as covariates.

reg_cases <- stan_glm(
  daily_confirmed_cases_per_mil ~ closed * country,
  data = panel,
  family = gaussian(),
  seed = 42,
  refresh = 0
)

reg_cases_mob <- stan_glm(
  daily_confirmed_cases_per_mil ~ closed * country + retail_rec * grocery_pharm,
  data = panel,
  family = gaussian(),
  seed = 42,
  refresh = 0
)

print(reg_cases, digits = 3)

stan_glm
 family:       gaussian [identity]
 formula:      daily_confirmed_cases_per_mil ~ closed * country
 observations: 420
 predictors:   4
------
                     Median MAD_SD
(Intercept)           5.704  4.081
closed               11.787  4.266
countrySweden        -1.649  5.738
closed:countrySweden 31.507  6.201

Auxiliary parameter(s):
      Median MAD_SD
sigma 23.525  0.802

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

print(reg_cases_mob, digits = 3)

stan_glm
 family:       gaussian [identity]
 formula:      daily_confirmed_cases_per_mil ~ closed * country + retail_rec * 
       grocery_pharm
 observations: 420
 predictors:   7
------
                         Median MAD_SD
(Intercept)               7.230  3.733
closed                    7.864  4.076
countrySweden            -7.199  5.332
retail_rec               -0.674  0.098
grocery_pharm             1.280  0.199
closed:countrySweden     35.808  5.611
retail_rec:grocery_pharm  0.016  0.005

Auxiliary parameter(s):
      Median MAD_SD
sigma 21.998  0.775

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

diag_cases <- tibble(
  model = c("DiD", "DiD + mobility"),
  max_rhat = c(max(rhat(reg_cases), na.rm = TRUE), max(rhat(reg_cases_mob), na.rm = TRUE)),
  min_neff_ratio = c(
    min(neff_ratio(reg_cases), na.rm = TRUE),
    min(neff_ratio(reg_cases_mob), na.rm = TRUE)
  )
)
diag_cases

# A tibble: 2 × 3
  model          max_rhat min_neff_ratio
  <chr>             <dbl>          <dbl>
1 DiD                1.00          0.334
2 DiD + mobility     1.00          0.472

posterior_interval(reg_cases, prob = 0.9)

                             5%      95%
(Intercept)           -1.002335 12.22430
closed                 4.738444 19.03943
countrySweden        -11.018812  7.72573
closed:countrySweden  21.321477 41.93866
sigma                 22.262350 24.91557

posterior_interval(reg_cases_mob, prob = 0.9)

                                    5%         95%
(Intercept)                0.911912142 13.62949550
closed                     0.858495245 14.86935987
countrySweden            -16.213380644  1.94812856
retail_rec                -0.837353127 -0.50703969
grocery_pharm              0.941374972  1.61682293
closed:countrySweden      26.125772471 45.76337813
retail_rec:grocery_pharm   0.008065712  0.02397155
sigma                     20.792979235 23.38723987

pp_check(reg_cases) + ggtitle("PPC: daily cases per million (DiD)")

reg_deaths <- stan_glm(
  daily_deaths_per_mil ~ closed * country,
  data = panel,
  family = gaussian(),
  seed = 42,
  refresh = 0
)

reg_deaths_mob <- stan_glm(
  daily_deaths_per_mil ~ closed * country + retail_rec * grocery_pharm,
  data = panel,
  family = gaussian(),
  seed = 42,
  refresh = 0
)

print(reg_deaths, digits = 3)

stan_glm
 family:       gaussian [identity]
 formula:      daily_deaths_per_mil ~ closed * country
 observations: 420
 predictors:   4
------
                     Median MAD_SD
(Intercept)          0.033  0.345 
closed               0.583  0.391 
countrySweden        0.020  0.498 
closed:countrySweden 2.612  0.524 

Auxiliary parameter(s):
      Median MAD_SD
sigma 2.075  0.070 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

print(reg_deaths_mob, digits = 3)

stan_glm
 family:       gaussian [identity]
 formula:      daily_deaths_per_mil ~ closed * country + retail_rec * grocery_pharm
 observations: 420
 predictors:   7
------
                         Median MAD_SD
(Intercept)              -0.015  0.326
closed                    0.727  0.361
countrySweden            -0.192  0.472
retail_rec               -0.067  0.009
grocery_pharm             0.034  0.018
closed:countrySweden      2.659  0.518
retail_rec:grocery_pharm  0.000  0.000

Auxiliary parameter(s):
      Median MAD_SD
sigma 1.860  0.066 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

diag_deaths <- tibble(
  model = c("DiD", "DiD + mobility"),
  max_rhat = c(max(rhat(reg_deaths), na.rm = TRUE), max(rhat(reg_deaths_mob), na.rm = TRUE)),
  min_neff_ratio = c(
    min(neff_ratio(reg_deaths), na.rm = TRUE),
    min(neff_ratio(reg_deaths_mob), na.rm = TRUE)
  )
)
diag_deaths

# A tibble: 2 × 3
  model          max_rhat min_neff_ratio
  <chr>             <dbl>          <dbl>
1 DiD                1.00          0.345
2 DiD + mobility     1.00          0.498

posterior_interval(reg_deaths, prob = 0.9)

                              5%       95%
(Intercept)          -0.56408869 0.5931228
closed               -0.02948765 1.2326987
countrySweden        -0.80928110 0.8609287
closed:countrySweden  1.69659041 3.5184665
sigma                 1.96725738 2.1948489

posterior_interval(reg_deaths_mob, prob = 0.9)

                                    5%           95%
(Intercept)              -0.5403777211  0.5234285830
closed                    0.1181225439  1.3005178006
countrySweden            -0.9952220459  0.5766464594
retail_rec               -0.0807647635 -0.0525108928
grocery_pharm             0.0059167308  0.0637226929
closed:countrySweden      1.8469136099  3.5260204477
retail_rec:grocery_pharm -0.0009485852  0.0004203318
sigma                     1.7579943764  1.9731282330

pp_check(reg_deaths) + ggtitle("PPC: daily deaths per million (DiD)")

mcmc_areas(as.matrix(reg_cases), pars = "closed:countrySweden") +
  labs(title = "Posterior: DiD interaction (cases per million)")

mcmc_areas(as.matrix(reg_deaths), pars = "closed:countrySweden") +
  labs(title = "Posterior: DiD interaction (deaths per million)")

Read the interaction posteriors as Sweden’s excess before→after change relative to Denmark, after the shared border-closure cutoff. Mobility covariates in the secondary models absorb some of the time-varying behavior associated with retail/grocery activity; the interaction remains the DiD contrast of interest.

Discussion and next steps

Our results speak to whether, irrespective of mobility curtailment and measures undertaken — from border closures to restaurant shutdowns — both countries’ health metrics eventually moved in ways that leave only marginal differences through time, while much of the observed contrast sits in diminished mobility and increased economic disruption. Next up, it would be prudent to examine data pertaining to the characteristics of the deceased and infected. Should additional data confirm what the reportage suggests, i.e. that the majority of the deaths occurred concomitant with underlying health issues, especially in higher age brackets, the rationale for mobility restrictions would face serious credibility problems.