library(tidyverse); library(knitr); library(kableExtra); library(broom); library(cowplot); 
library(rstan); library(tidybayes); library(scales); library(DT); library(plotly)



theme_set(theme_classic(base_size = 12) + 
            background_grid(color.major = "grey90", 
                            color.minor = "grey95", 
                            minor = "xy", major = "xy") +
            theme(legend.position = "none"))
m <- read_rds("Hierarchical_Model_NegBin.rds")
d <- read_csv("Interactive Model Checking/stan_data.csv")
countries <- d %>% distinct(country, country_id)
which_iceland <- d %>% filter(country == "Iceland") %>% .$country_id %>% unique
n_countries <- max(d$country_id)
results <- tidyMCMC(m, conf.int = T, rhat = T, ess = T, estimate.method = "median", conf.method = "quantile") %>% 
  mutate(par = str_match(term, "[a-zA-Z_2]+")) %>% 
  group_by(par) %>% 
  mutate(num = row_number() %>% as.numeric)

Updates

2020-03-27

  • We now place informative priors on \(\mu_\alpha\), \(\mu_\beta\) and \(\lambda_\phi\) to ensure that the model can be identified.

2020-03-26

  • Implemented a negative binomial likelihood with country-specific overdispersion.
  • Changed filtering criteria for inclusion in model data in order to include more countries.
  • Did some visual inspection of data in order to make sure we don’t lose countries to inner_join strategies. Saw that we were missing USA among others due to different names in different datasets, but is now fixed.

Methods

Parametrization

Let \(E_i\) and \(I_{i, t}\) be the population and number of infected in country \(i\) at time \(t\). Then the percent of infected can be calculated as

\[ P_{i, t} = \frac{I_{i, t}}{E_i}. \]

In the usual Logistic Regression GLM we could model the percent of infected, as a function of time (in days), with

\[ \log\left(\frac{P_{i, t}}{1 - P_{i, t}}\right) = \alpha_i + \beta_i \cdot t, \]

where \(\alpha_i\) is a measure of how many have been infected in country \(i\) at time \(t = 0\) and \(\beta_i\) is a measure of growth. In the case of COVID-19 infections we don’t know the maximum percent of populations that will be infected, so we have another unknown parameter, the saturation percent at which a country will reach its maximum number of infected, \(S_i\). Thus our model looks like

\[ \log\left(\frac{P_{i, t}}{S_i - P_{i, t}}\right) = \alpha_i + \beta_i \cdot t. \]

These parameters are hard to estimate when data from only one country are used. However, if we were to pool information about them between countries, as in a hierarchical Bayesian model, estimation might be possible. Let

\[ z_{i, t} = \alpha_i + \beta_i \cdot t, \]

where \(\alpha_i\) is a measure of how many have been infected in country \(i\) at time \(t = 0\) and \(\beta_i\) is a measure of growth. Then

\[ P_{i, t} = \frac{S_i}{1 + \exp(-z_{i, t})}, \]

and conditional on some sampling distribution, \(f\), we could write

\[ I_{i, t} \sim \mathrm{f_\theta}(P_{i, t}, E_i), \]

where \(\theta\) contains all relevant parameters to be estimated.

Bayesian Inference

Bayesian inference is a great tool when small amounts of data are to be shared from various sources. In this case the sources are different countries, and the data are cumulative numbers of cases. If we utilize a Negative Binomial likelihood for the observed cumulative cases, then

\[ I_{i, t} \sim \mathrm{NegBin}(P_{i, t} \cdot E_i, \phi_i), \]

where \(\phi_i\) is a country-level effect specifying the amount of overdispersion in that country.

Modeling daily counts

However, there is a lot of implicit correlation in the values of \(I_{i, t}\) for different values of \(t\). Thus, a better parametrization would be to model the daily number of cases.

Let

\[ z_{i, t} = \alpha_i + \beta_i \cdot t, \]

so that the percent of infected, \(P_i\), is

\[ P_{i, t} = \frac{S_i}{1 + \exp(-z_{i,t})}. \]

If we furthermore write

\[ z^*_{i, t - 1} = \alpha_i + \beta_i \cdot (t - 1), \]

and

\[ P^*_{i, t - 1} = \frac{S_i}{1 + \exp(-z^*_{i, t-1})}, \]

the change in rates between days is

\[ C_{i, t} = P_{i, t} - P^*_{i, t - 1}. \]

Since \(C_{i, t}\) is simply the first derivative of \(P_{i, t}\) with respect to \(t\), we can skip the differencing step and directly model the derivative

\[ C_{i, t} = \frac{d}{dt}P_{i, t} = \beta_i S_i \frac{\exp{(-z_{i, t})}}{(\exp(-z_{i, t}) + 1))^2} \]

Then, conditional on a Negative Binomial likelihood and population size \(E_i\), the daily number of observed cases, \(D_{i, t}\), can be written as

\[ D_{i, t} \sim \mathrm{NegBin}(C_{i, t} \cdot E_i, \phi_i) \]

Parameters

The parameters, \(\alpha\) and \(\beta\), are treated as in a generalized linear model, except that we use the constraints \(\beta_i > 0\) and \(\mu_\beta > 0\). We put hierarchical priors on \(\alpha_i\) and \(\beta_i\) so that for each country, \(i\),

\[ \begin{aligned} \beta_i &\sim \mathrm{Normal}(\mu_\beta, \sigma^2_\beta) \\ \alpha_i &\sim \mathrm{Normal}(\mu_\alpha, \sigma^2_\alpha) \end{aligned} \]

The \(\mu\) parameters are given informed prior distribution based on our previous Poisson model

\[ \begin{aligned} \mu_\alpha &\sim \mathrm{Normal}(-3.6 , 0.5^2) \\ \mu_\beta &\sim \mathrm{Normal}(0.13, 0.05^2) \end{aligned} \]

We chose to put priors on the standard deviations, so that \(\sigma_\beta\) and \(\sigma_\alpha\) are given \(\mathrm{Exponential}\) prior distributions

\[ \begin{aligned} \sigma_\alpha &\sim \mathrm{Exponential}(1) \\ \sigma_\beta &\sim \mathrm{Exponential}(2). \end{aligned} \]

The \(S_i\) parameters take on values in \((0, 1)\), so we thought it proper to model them with Beta distributions. By putting hierarchical priors on them we could also share information between countries on the estimated saturation points. Thus the saturation parameter for country \(i\) is sampled as

\[ S_i \sim \mathrm{Beta}(a_S, b_S), \]

where \(a_S\) and \(b_S\) are given vague Half-Cauchy prior distributions to help with the posterior geometry

\[ \begin{aligned} a_S &\sim \mathrm{Cauchy}(0, 5) \\ b_s &\sim \mathrm{Cauchy}(0, 5). \end{aligned} \]

We parametrize the negative binomial likelihood in the form of mean and overdispersion. If we write

\[ D_{i, t} \sim \mathrm{NegBin}(\mu_{i, t}, \phi_i), \]

where

\[ \mu_{i, t} = C_{i, t} \cdot E_i, \]

we can write

\[ E[D_{i, t}] = \mu_{i, t} \qquad Var[D_{i, t}] = \mu_{i, t} + \phi_i \cdot \mu_{i, t}^2. \]

Based on Dan Simpson’s excellent post, which is also linked in Stan-Dev’s post on prior choices we put a hierarchical exponential prior on the \(\phi_i\) parameters as follows:

\[ \begin{aligned} \sqrt{\phi_i} &\sim \mathrm{Normal}(0, \sigma_\phi^2) \\ \sigma_\phi &\sim \mathrm{Exponential}(1) \end{aligned} \]

Putting it all together we get

\[ \begin{aligned} D_{i, t} &\sim \mathrm{NegBin}(C_{i, t} \cdot E_i, \phi_i) \\ C_{i, t} = \frac{d}{dt}P_{i, t} &= \beta_i S_i \frac{\exp{(-z_{i, t})}}{(\exp(-z_{i, t}) + 1))^2} \\ z_{i, t} &= \alpha_i + \beta_i \cdot t \\ \beta_i &\sim \mathrm{Normal}(\mu_\beta, \sigma^2_\beta) \\ \alpha_i &\sim \mathrm{Normal}(\mu_\alpha, \sigma^2_\alpha) \\ \mu_\beta &\sim \mathrm{Normal}(0.13, 0.05^2) \\ \mu_\alpha &\sim \mathrm{Normal}(-3.6, 0.5^2) \\ p(\sigma_\alpha) &\sim \mathrm{Exponential}(1) \\ p(\sigma_\beta) &\sim \mathrm{Exponential}(2) \\ \sqrt\phi_i &\sim \mathrm{Normal}(0, \lambda_\phi^2)\\ \lambda_\phi &\sim \mathrm{Exponential}(1) \\ S_i &\sim \mathrm{Beta}(a_S, b_S) \\ p(a_S) &\sim \mathrm{Cauchy}(0, 5) \\ p(b_S) &\sim \mathrm{Cauchy}(0, 5) \end{aligned} \]

Data

Data used for modeling were filtered according to

  • Only keep data for which cases per 1000 inhabitants \(\geq\) 0.02.
  • Only keep data for countries for which the first observation after this filtering is \(\leq\) 0.06
  • Only keep data for countries with more than 4 days of follow-up after the above steps.
d %>% 
  group_by(country) %>% 
  summarise(First = min(date),
            Days_In_Data = n(),
            Start_Rate = min(case_rate),
            End_Rate = max(case_rate)) %>% 
  set_names(c("Country", "Entry", "Days in data", "Start", "End")) %>% 
  kable(caption = "Table 1. Summary information about countries used in modeling",
        align = c("l", rep("c", ncol(.) - 1)),
        digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(which_iceland, bold = T) %>% 
  add_header_above(c("", "", "", "Rate per 1000" = 2)) %>%
  scroll_box(height = "500px")
Table 1. Summary information about countries used in modeling
Rate per 1000
Country Entry Days in data Start End
Albania 2020-03-19 10 0.020 0.065
Andorra 2020-03-14 14 0.026 3.461
Armenia 2020-03-18 11 0.026 0.126
Aruba 2020-03-20 6 0.038 0.263
Australia 2020-03-19 10 0.022 0.134
Austria 2020-03-11 18 0.020 0.860
Bahrain 2020-02-27 30 0.020 0.284
Barbados 2020-03-21 8 0.021 0.084
Belgium 2020-03-10 19 0.021 0.631
Bermuda 2020-03-20 9 0.032 0.272
Bosnia And Herzegovina 2020-03-22 7 0.028 0.069
Bulgaria 2020-03-22 7 0.023 0.042
Canada 2020-03-20 9 0.023 0.125
Cayman Islands 2020-03-21 8 0.046 0.123
Chile 2020-03-21 8 0.023 0.085
Costa Rica 2020-03-21 8 0.022 0.052
Croatia 2020-03-20 9 0.025 0.142
Cyprus 2020-03-16 13 0.025 0.137
Czech Republic 2020-03-15 14 0.020 0.213
Denmark 2020-03-11 18 0.046 0.354
Dominican Republic 2020-03-24 5 0.023 0.054
Ecuador 2020-03-21 8 0.025 0.094
Estonia 2020-03-13 16 0.020 0.434
Finland 2020-03-13 15 0.028 0.185
France 2020-03-10 19 0.022 0.506
French Polynesia 2020-03-21 8 0.054 0.107
Germany 2020-03-13 16 0.028 0.582
Greece 2020-03-15 14 0.022 0.092
Greenland 2020-03-20 9 0.035 0.159
Iceland 2020-03-02 26 0.029 2.840
Iran 2020-03-04 25 0.028 0.390
Ireland 2020-03-15 14 0.026 0.434
Isle Of Man 2020-03-22 7 0.024 0.343
Israel 2020-03-15 14 0.021 0.356
Italy 2020-03-02 27 0.028 1.429
Kuwait 2020-03-14 15 0.024 0.053
Latvia 2020-03-18 11 0.032 0.147
Lebanon 2020-03-20 9 0.022 0.057
Liechtenstein 2020-03-05 18 0.026 1.578
Lithuania 2020-03-21 8 0.025 0.130
Luxembourg 2020-03-13 15 0.042 2.607
Malaysia 2020-03-18 11 0.021 0.068
Maldives 2020-03-16 13 0.024 0.026
Malta 2020-03-13 16 0.020 0.316
Mauritius 2020-03-23 6 0.023 0.078
Moldova 2020-03-23 6 0.023 0.049
Monaco 2020-02-29 17 0.026 1.078
Montenegro 2020-03-20 9 0.021 0.119
Netherlands 2020-03-11 18 0.022 0.503
New Caledonia 2020-03-24 5 0.025 0.053
New Zealand 2020-03-23 6 0.021 0.087
North Macedonia 2020-03-19 10 0.020 0.105
Norway 2020-03-07 22 0.021 0.666
Panama 2020-03-18 11 0.020 0.185
Portugal 2020-03-16 13 0.024 0.417
Romania 2020-03-23 6 0.022 0.067
San Marino 2020-02-28 29 0.030 6.586
Serbia 2020-03-23 6 0.021 0.052
Seychelles 2020-03-15 14 0.020 0.072
Singapore 2020-03-06 23 0.020 0.126
Slovakia 2020-03-20 9 0.023 0.054
Slovenia 2020-03-12 17 0.027 0.304
South Korea 2020-02-26 32 0.022 0.185
Spain 2020-03-10 19 0.026 1.371
Sweden 2020-03-09 20 0.020 0.303
Switzerland 2020-03-07 22 0.024 1.409
Trinidad And Tobago 2020-03-22 7 0.035 0.047
United Arab Emirates 2020-03-24 5 0.020 0.041
United Kingdom 2020-03-16 13 0.021 0.215
United States 2020-03-19 10 0.029 0.318
Uruguay 2020-03-19 10 0.023 0.069 
p <- d %>% 
  ggplot(aes(days, case_rate, group = country, col = country == "Iceland")) +
  geom_line() +
  scale_y_log10() +
  scale_colour_manual(values = c("grey", "blue")) +
  labs(x = "Days since rate reached 0.02 per 1000",
       y = "Cases per 1000",
       title = "Observed trends for countries used in data",
       subtitle = "Shown as days since a country entered the modeling data")
ggplotly(p)

Software

The model is fit using Stan’s R interface, and the model code can be found in the appendix.

All code is available at https://github.com/bgautijonsson/covid19. We’re sorry for the mess and we’re working on it. We are also going to translate the README files into English.

Results

Convergence

Sampler configuration

  • Four chains were run for 2000 iterations after a warm-up period of 2000 iterations. The slowest chain converged in 430 seconds.
  • Adapt_delta was set to 0.99

Warnings

  • There were 535 divergent transitions after warmup.

Country Level Effects

results %>% 
  ungroup %>% 
  filter(par %in% c("beta", "alpha", "S", "phi_inv")) %>%
  mutate(par = str_replace(par, "phi_inv", "phi")) %>% 
  inner_join(countries, by = c("num" = "country_id")) %>% 
  mutate(par = str_to_title(par)) %>% 
  select(par, country, num, everything(), -num, -term, -std.error) %>% 
  set_names(c("Parameter", "Country", "Median", "Lower", "Upper", "Rhat", "ESS")) %>% 
  kable(digits = 4, align = c("l", "l", rep("c", ncol(.) - 2)),
        caption = "Table 2. Summary of posterior samples of country level effects") %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  add_header_above(c("", "", "", "95% PI" = 2, "Convergence" = 2)) %>% 
  column_spec(1, bold = T) %>% 
  row_spec(which_iceland + c(0, 1, 2) * n_countries, bold = T) %>% 
  collapse_rows(1, valign = "top") %>% 
  scroll_box(height = "600px")
Table 2. Summary of posterior samples of country level effects
95% PI
Convergence
Parameter Country Median Lower Upper Rhat ESS
Beta Albania 0.1320 0.0469 0.2220 0.9999 5999
Andorra 0.2149 0.1364 0.3096 0.9998 5391
Armenia 0.1064 0.0429 0.1895 1.0005 7021
Aruba 0.1536 0.0692 0.2554 1.0002 6589
Australia 0.1314 0.0544 0.2322 1.0001 7040
Austria 0.2108 0.1635 0.2620 1.0004 3821
Bahrain 0.0844 0.0408 0.1540 1.0002 5070
Barbados 0.1129 0.0379 0.2188 1.0006 6502
Belgium 0.2013 0.1648 0.2457 1.0001 5180
Bermuda 0.1546 0.0711 0.2529 1.0002 5864
Bosnia And Herzegovina 0.1181 0.0396 0.2224 1.0007 6712
Bulgaria 0.0634 0.0152 0.1548 1.0000 6262
Canada 0.1732 0.0826 0.2636 0.9996 6308
Cayman Islands 0.1376 0.0529 0.2422 0.9998 7177
Chile 0.1581 0.0919 0.2174 1.0004 5604
Costa Rica 0.1057 0.0367 0.1900 0.9999 7038
Croatia 0.1362 0.0629 0.2189 1.0003 6849
Cyprus 0.0999 0.0376 0.1806 1.0000 6753
Czech Republic 0.1270 0.0762 0.1853 1.0005 6949
Denmark 0.0543 0.0226 0.1159 1.0000 4430
Dominican Republic 0.1533 0.0685 0.2340 1.0001 5296
Ecuador 0.0916 0.0320 0.1789 0.9998 6265
Estonia 0.0927 0.0441 0.1703 0.9997 6485
Finland 0.0737 0.0312 0.1348 1.0001 6944
France 0.1582 0.1327 0.1945 1.0000 4317
French Polynesia 0.1139 0.0376 0.2154 1.0003 7926
Germany 0.1745 0.1192 0.2463 1.0003 5149
Greece 0.0639 0.0241 0.1200 1.0006 6589
Greenland 0.1404 0.0581 0.2381 1.0001 6524
Iceland 0.1673 0.1296 0.2072 1.0003 4884
Iran 0.0547 0.0391 0.0767 1.0004 5435
Ireland 0.1741 0.1294 0.2315 1.0005 4480
Isle Of Man 0.1609 0.0759 0.2638 0.9998 7610
Israel 0.1607 0.0881 0.2420 1.0005 6333
Italy 0.1599 0.1134 0.2082 0.9998 5359
Kuwait 0.0404 0.0099 0.1040 1.0000 5646
Latvia 0.0938 0.0405 0.1610 1.0001 6513
Lebanon 0.0857 0.0236 0.1731 1.0011 6501
Liechtenstein 0.1745 0.0972 0.2745 1.0002 5170
Lithuania 0.1214 0.0543 0.2021 1.0004 7467
Luxembourg 0.2336 0.1550 0.3237 1.0003 4200
Malaysia 0.0639 0.0262 0.1153 0.9999 7808
Maldives 0.0419 0.0047 0.1476 1.0002 6511
Malta 0.1110 0.0523 0.2788 1.0030 1543
Mauritius 0.1351 0.0527 0.2330 1.0000 6955
Moldova 0.1201 0.0428 0.2092 0.9998 6349
Monaco 0.1614 0.0963 0.2434 1.0001 6414
Montenegro 0.1452 0.0618 0.2424 1.0000 7359
Netherlands 0.1773 0.1526 0.2103 0.9998 3856
New Caledonia 0.1177 0.0376 0.2207 1.0001 8389
New Zealand 0.1461 0.0651 0.2344 1.0003 6363
North Macedonia 0.1049 0.0424 0.1825 0.9998 6060
Norway 0.1292 0.0775 0.2027 1.0002 6075
Panama 0.1711 0.1040 0.2405 1.0001 6001
Portugal 0.2003 0.1597 0.2466 1.0006 4654
Romania 0.1387 0.0549 0.2345 0.9997 6182
San Marino 0.1659 0.1097 0.2260 1.0003 5888
Serbia 0.1433 0.0602 0.2312 1.0001 6088
Seychelles 0.0624 0.0162 0.1479 1.0003 6938
Singapore 0.1080 0.0733 0.1454 1.0005 6339
Slovakia 0.1019 0.0288 0.1953 0.9999 6556
Slovenia 0.0517 0.0251 0.1026 1.0002 4980
South Korea 0.1136 0.0806 0.1590 1.0007 2273
Spain 0.1939 0.1667 0.2275 1.0002 6132
Sweden 0.0810 0.0503 0.1299 1.0012 5008
Switzerland 0.1767 0.1163 0.2447 0.9996 6475
Trinidad And Tobago 0.1052 0.0289 0.2099 1.0000 6979
United Arab Emirates 0.1224 0.0355 0.2342 1.0002 7439
United Kingdom 0.1901 0.1396 0.2417 1.0002 6187
United States 0.2015 0.1564 0.2515 0.9998 3866
Uruguay 0.0867 0.0245 0.1784 1.0007 7232
Alpha Albania -4.4836 -5.8582 -2.9738 1.0002 3868
Andorra -2.9658 -4.0032 -1.8188 1.0005 6076
Armenia -3.7871 -5.0257 -2.3545 1.0001 3995
Aruba -3.0118 -4.2729 -1.5625 1.0007 3877
Australia -3.8526 -5.1933 -2.3428 1.0004 3770
Austria -4.0512 -4.7587 -3.5929 1.0002 3667
Bahrain -4.0314 -5.3528 -2.7273 0.9999 4485
Barbados -3.6510 -4.9880 -2.1452 1.0004 4234
Belgium -4.7645 -5.6409 -4.0867 1.0002 3756
Bermuda -3.4840 -4.7681 -2.0101 1.0010 4739
Bosnia And Herzegovina -3.7892 -5.1189 -2.3456 1.0000 5065
Bulgaria -3.9117 -5.2925 -2.3329 1.0007 3404
Canada -4.2407 -5.5571 -2.8797 0.9999 4580
Cayman Islands -3.7856 -5.1664 -2.3399 0.9998 4384
Chile -4.2943 -5.4266 -3.0228 1.0016 3375
Costa Rica -4.2327 -5.5031 -2.7944 0.9997 3866
Croatia -3.7864 -4.9926 -2.3808 0.9998 3931
Cyprus -3.8049 -5.1494 -2.2848 1.0010 4264
Czech Republic -3.9381 -5.0557 -2.7361 1.0005 3544
Denmark -2.5308 -3.8284 -0.8942 1.0004 3094
Dominican Republic -4.1703 -5.3095 -2.8616 1.0006 3244
Ecuador -3.4933 -4.7350 -1.9448 1.0002 4642
Estonia -2.9726 -4.1737 -1.6196 0.9998 3863
Finland -3.3446 -4.5456 -1.9384 1.0013 3502
France -4.1202 -5.0192 -3.5730 1.0003 2646
French Polynesia -3.5075 -4.7668 -1.9951 1.0001 4275
Germany -3.9045 -4.9129 -2.9811 1.0004 4271
Greece -3.7663 -5.0113 -2.3242 0.9999 4207
Greenland -3.6843 -4.9995 -2.1948 1.0003 3858
Iceland -4.1103 -4.6779 -3.5534 0.9999 6781
Iran -3.1773 -4.1087 -2.2986 1.0000 3191
Ireland -3.8881 -4.8526 -3.0653 1.0001 3709
Isle Of Man -3.0681 -4.3137 -1.5748 0.9998 4195
Israel -3.9744 -5.1786 -2.6725 1.0005 4803
Italy -4.0714 -4.7532 -3.4012 1.0004 7355
Kuwait -3.9089 -5.3668 -2.3502 1.0000 4364
Latvia -3.5399 -4.7166 -2.1777 1.0002 3871
Lebanon -4.0484 -5.3899 -2.5241 0.9998 4133
Liechtenstein -3.1635 -4.3621 -1.8617 1.0001 5134
Lithuania -3.5867 -4.7609 -2.2117 1.0005 3291
Luxembourg -3.3920 -4.2410 -2.4724 1.0007 5786
Malaysia -3.8214 -5.0704 -2.3990 0.9999 3459
Maldives -4.3795 -6.0359 -2.6475 0.9996 4930
Malta -3.0808 -4.4616 -1.8780 1.0001 2392
Mauritius -3.7862 -5.0379 -2.3527 1.0005 3971
Moldova -4.2103 -5.4743 -2.7709 1.0006 3989
Monaco -3.4887 -4.6457 -2.2638 1.0010 5712
Montenegro -3.9643 -5.2409 -2.5129 1.0001 4739
Netherlands -4.2784 -5.1487 -3.7934 1.0012 2774
New Caledonia -3.8651 -5.1946 -2.3586 1.0000 4576
New Zealand -3.8290 -4.9670 -2.4762 1.0005 4663
North Macedonia -3.7943 -5.0547 -2.3699 1.0000 3501
Norway -3.6947 -4.7689 -2.6347 1.0007 4641
Panama -4.1748 -5.3377 -2.8955 1.0002 4392
Portugal -4.1642 -5.1052 -3.4025 1.0010 3170
Romania -4.0441 -5.2488 -2.6464 0.9999 4112
San Marino -3.3475 -4.3019 -2.3553 1.0004 6354
Serbia -4.2790 -5.5261 -2.8605 1.0008 3702
Seychelles -3.7331 -5.1936 -2.0356 1.0000 4632
Singapore -4.8260 -6.1067 -3.5614 1.0003 3534
Slovakia -4.2371 -5.6551 -2.7066 0.9998 4329
Slovenia -2.6456 -3.8217 -1.3404 1.0012 3611
South Korea -0.1533 -1.3549 1.4050 1.0006 1865
Spain -4.3218 -4.9059 -3.8946 1.0007 4555
Sweden -3.4677 -4.5274 -2.3699 1.0003 3308
Switzerland -4.1032 -5.0423 -3.0924 1.0001 5607
Trinidad And Tobago -3.7788 -5.1145 -2.2470 0.9997 5029
United Arab Emirates -3.9780 -5.3797 -2.4570 1.0004 5105
United Kingdom -4.5651 -5.6428 -3.4712 1.0004 3993
United States -3.9122 -4.8696 -3.0270 1.0002 3369
Uruguay -3.8470 -5.2258 -2.2216 1.0008 3828
S Albania 0.0019 0.0004 0.0059 1.0002 3544
Andorra 0.0064 0.0035 0.0133 1.0006 3100
Armenia 0.0028 0.0007 0.0080 1.0005 3607
Aruba 0.0041 0.0014 0.0108 1.0015 2994
Australia 0.0028 0.0007 0.0082 1.0005 4132
Austria 0.0021 0.0013 0.0057 1.0009 2436
Bahrain 0.0020 0.0004 0.0068 1.0009 3795
Barbados 0.0030 0.0008 0.0086 1.0002 3760
Belgium 0.0024 0.0012 0.0063 1.0003 2781
Bermuda 0.0035 0.0011 0.0092 1.0009 3741
Bosnia And Herzegovina 0.0029 0.0008 0.0082 1.0004 4014
Bulgaria 0.0026 0.0006 0.0077 1.0011 4282
Canada 0.0024 0.0007 0.0069 1.0001 3674
Cayman Islands 0.0029 0.0008 0.0083 1.0000 3725
Chile 0.0023 0.0007 0.0067 1.0013 3444
Costa Rica 0.0023 0.0006 0.0066 1.0000 3278
Croatia 0.0028 0.0008 0.0079 1.0001 3395
Cyprus 0.0027 0.0007 0.0078 1.0002 4031
Czech Republic 0.0027 0.0008 0.0073 1.0004 3502
Denmark 0.0039 0.0010 0.0107 1.0007 3087
Dominican Republic 0.0024 0.0007 0.0067 1.0001 3462
Ecuador 0.0032 0.0009 0.0090 1.0003 3279
Estonia 0.0036 0.0011 0.0100 1.0001 3112
Finland 0.0032 0.0008 0.0090 1.0009 3153
France 0.0023 0.0011 0.0063 1.0005 2292
French Polynesia 0.0033 0.0010 0.0087 1.0003 3333
Germany 0.0027 0.0010 0.0076 1.0004 2957
Greece 0.0027 0.0007 0.0078 1.0001 3673
Greenland 0.0030 0.0009 0.0088 1.0006 3241
Iceland 0.0052 0.0036 0.0094 1.0007 3193
Iran 0.0039 0.0014 0.0098 1.0004 2984
Ireland 0.0028 0.0010 0.0076 1.0003 3121
Isle Of Man 0.0041 0.0015 0.0102 1.0003 3536
Israel 0.0028 0.0009 0.0075 0.9996 3560
Italy 0.0027 0.0017 0.0063 1.0012 3288
Kuwait 0.0026 0.0006 0.0080 1.0001 5176
Latvia 0.0031 0.0009 0.0086 1.0006 2780
Lebanon 0.0024 0.0005 0.0073 1.0016 3501
Liechtenstein 0.0038 0.0017 0.0095 1.0005 3409
Lithuania 0.0031 0.0009 0.0088 1.0006 3723
Luxembourg 0.0051 0.0030 0.0108 1.0011 3128
Malaysia 0.0027 0.0007 0.0079 1.0004 3650
Maldives 0.0020 0.0003 0.0068 1.0009 4076
Malta 0.0023 0.0004 0.0082 1.0011 2069
Mauritius 0.0030 0.0009 0.0082 1.0005 3355
Moldova 0.0023 0.0006 0.0067 1.0002 4051
Monaco 0.0038 0.0017 0.0094 1.0026 3326
Montenegro 0.0026 0.0007 0.0076 1.0000 3459
Netherlands 0.0022 0.0012 0.0059 1.0015 2481
New Caledonia 0.0027 0.0007 0.0079 1.0002 3948
New Zealand 0.0029 0.0009 0.0080 1.0002 3942
North Macedonia 0.0028 0.0007 0.0080 1.0008 3490
Norway 0.0027 0.0010 0.0080 1.0006 3078
Panama 0.0025 0.0007 0.0071 1.0006 3470
Portugal 0.0028 0.0012 0.0073 1.0004 2827
Romania 0.0026 0.0007 0.0071 1.0001 3842
San Marino 0.0078 0.0052 0.0128 1.0005 3541
Serbia 0.0023 0.0006 0.0067 1.0006 3585
Seychelles 0.0028 0.0006 0.0082 1.0002 4087
Singapore 0.0015 0.0004 0.0049 1.0004 3568
Slovakia 0.0022 0.0005 0.0066 1.0000 3870
Slovenia 0.0040 0.0010 0.0108 1.0013 3187
South Korea 0.0003 0.0002 0.0009 1.0007 1748
Spain 0.0043 0.0026 0.0085 1.0014 3485
Sweden 0.0030 0.0008 0.0084 1.0004 3268
Switzerland 0.0033 0.0017 0.0083 1.0006 3225
Trinidad And Tobago 0.0029 0.0008 0.0084 0.9999 3668
United Arab Emirates 0.0026 0.0006 0.0077 1.0000 3618
United Kingdom 0.0022 0.0007 0.0063 1.0001 3404
United States 0.0030 0.0011 0.0079 1.0001 2887
Uruguay 0.0026 0.0006 0.0082 1.0004 3364
Phi Albania 0.3413 0.0697 1.1530 0.9999 7328
Andorra 0.8871 0.3815 2.0713 1.0000 8489
Armenia 0.4365 0.1818 1.1083 0.9999 7474
Aruba 0.6294 0.0656 2.2710 1.0001 9047
Australia 1.2303 0.6056 2.7074 0.9998 9093
Austria 0.0601 0.0291 0.1402 0.9998 5770
Bahrain 1.0134 0.5914 1.7470 0.9999 10124
Barbados 1.2773 0.3169 3.5522 1.0001 9116
Belgium 0.0926 0.0463 0.2131 0.9996 7316
Bermuda 1.3496 0.2673 3.9355 0.9997 10276
Bosnia And Herzegovina 1.1993 0.4869 2.9722 1.0003 8292
Bulgaria 0.0817 0.0017 0.5140 1.0010 6337
Canada 0.3205 0.1264 0.9162 0.9999 6419
Cayman Islands 0.3854 0.0011 2.7661 1.0000 9145
Chile 0.0417 0.0106 0.2064 1.0000 4516
Costa Rica 0.1513 0.0125 0.6731 0.9999 7043
Croatia 0.1963 0.0732 0.6169 1.0002 6275
Cyprus 0.5453 0.1649 1.5338 0.9999 8128
Czech Republic 0.1699 0.0810 0.4149 1.0001 6636
Denmark 0.4014 0.2064 0.8658 0.9997 6713
Dominican Republic 0.0140 0.0001 0.2294 1.0010 4043
Ecuador 0.1995 0.0671 0.6801 0.9996 5943
Estonia 0.4657 0.2307 0.9834 0.9997 8802
Finland 0.2952 0.1330 0.7174 1.0001 6778
France 0.0349 0.0182 0.0784 0.9999 7386
French Polynesia 1.1982 0.3510 3.2280 0.9998 9127
Germany 0.1876 0.0961 0.4212 1.0000 6923
Greece 0.2071 0.0899 0.5093 1.0001 7257
Greenland 0.1820 0.0006 1.6772 0.9999 9634
Iceland 0.1536 0.0769 0.3133 1.0004 8208
Iran 0.0702 0.0401 0.1355 0.9998 7301
Ireland 0.0675 0.0278 0.1868 1.0000 6384
Isle Of Man 1.1309 0.3232 3.2632 0.9998 9780
Israel 0.6537 0.3379 1.3635 1.0005 8697
Italy 0.2925 0.1782 0.5307 1.0000 8131
Kuwait 0.3938 0.1395 1.0058 1.0001 7892
Latvia 0.1275 0.0286 0.4243 1.0004 7173
Lebanon 0.2836 0.1013 0.8465 0.9998 7433
Liechtenstein 1.3227 0.5624 2.8958 0.9996 10258
Lithuania 0.1051 0.0243 0.4466 0.9999 6200
Luxembourg 0.3271 0.1530 0.7798 1.0000 7077
Malaysia 0.0625 0.0234 0.2101 1.0004 5696
Maldives 1.6402 0.0260 6.2196 0.9999 7253
Malta 0.1130 0.0014 0.4202 1.0002 4329
Mauritius 0.3361 0.0788 1.2637 1.0003 7241
Moldova 0.0144 0.0001 0.2130 0.9999 7194
Monaco 0.3132 0.0050 1.3514 0.9999 8934
Montenegro 0.4516 0.0837 1.5848 1.0001 8380
Netherlands 0.0208 0.0086 0.0549 1.0002 6639
New Caledonia 0.3172 0.0012 1.9866 1.0000 8523
New Zealand 0.0744 0.0118 0.4269 0.9997 5985
North Macedonia 0.1185 0.0145 0.4564 0.9999 7723
Norway 0.5452 0.3020 1.0772 0.9996 8383
Panama 0.1314 0.0484 0.3904 1.0001 6517
Portugal 0.0386 0.0155 0.1182 1.0005 4831
Romania 0.1297 0.0368 0.5895 1.0001 4915
San Marino 1.0518 0.5829 1.9479 0.9997 10271
Serbia 0.0744 0.0050 0.4390 1.0002 5259
Seychelles 0.4172 0.0011 3.0311 0.9998 9839
Singapore 0.2270 0.1154 0.4782 0.9998 7644
Slovakia 0.5711 0.2393 1.4812 1.0001 7482
Slovenia 0.1062 0.0343 0.3028 1.0003 6876
South Korea 0.2553 0.1477 0.4820 1.0005 4304
Spain 0.0451 0.0237 0.0999 1.0012 6129
Sweden 0.1528 0.0823 0.3172 1.0000 7436
Switzerland 0.6067 0.3429 1.1487 0.9998 8357
Trinidad And Tobago 1.5410 0.7055 3.4637 1.0004 8031
United Arab Emirates 1.2866 0.4722 3.4591 1.0003 8561
United Kingdom 0.0919 0.0416 0.2502 1.0004 5852
United States 0.0236 0.0091 0.0845 1.0012 3733
Uruguay 0.7388 0.2765 1.9223 0.9999 7562 
p <- results %>% 
  ungroup %>% 
  filter(par %in% c("beta", "alpha", "S", "phi_inv")) %>% 
  mutate(par = str_replace(par, "phi_inv", "phi")) %>% 
  inner_join(countries, by = c("num" = "country_id")) %>% 
  mutate(par = str_to_title(par)) %>% 
  mutate(plot_var = str_c(par, "_", country)) %>% 
  mutate(country = fct_reorder(country, estimate)) %>% 
  ggplot(aes(country, estimate, ymin = conf.low, ymax = conf.high, col = country == "Iceland")) +
  geom_linerange() +
  geom_point() +
  coord_flip() +
  facet_wrap("par", scales = "free_x") +
  scale_colour_manual(values = c("grey", "blue")) +
  labs(title = "Posterior medians and 95% PIs") +
  theme(axis.title = element_blank())

ggplotly(p)

Hyperparameters

results %>% 
  ungroup %>% 
  filter(par %in% c("mu_beta", "sigma_sq_beta", 
                    "mu_alpha", "sigma_sq_alpha", 
                    "a_s", "b_s",
                    "sigma_phi_inv_sqrt")) %>% 
  mutate(par = str_replace(par, "sigma_phi_inv_sqrt", "sigma_phi_sqrt")) %>% 
  select(par, everything(), -num, -term, -std.error) %>% 
  set_names(c("Parameter",  "Median", "Lower", "Upper", "Rhat", "ESS")) %>% 
  kable(digits = 4, align = c("l", "l", rep("c", ncol(.) - 2)),
        caption = "Table 3. Summary of posterior samples of hyperparameters") %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  add_header_above(c("", "", "95% PI" = 2, "Convergence" = 2)) %>% 
  column_spec(1, bold = T)
Table 3. Summary of posterior samples of hyperparameters
95% PI
Convergence
Parameter Median Lower Upper Rhat ESS
mu_beta 0.1318 0.1132 0.1519 1.0028 2013
mu_alpha -3.7447 -4.1372 -3.3593 1.0014 1077
a_s 2.7974 1.6225 5.1249 1.0051 1324
b_s 847.0816 413.8872 1680.5013 1.0061 1315
sigma_phi_sqrt 0.7027 0.5790 0.8633 0.9999 7967

Posterior Predictions

We obtain predictions by conditioning on the first observed number of cumulative cases in a country (as seen in the modeling data), perform simulations of new daily cases and sum those up to get predicted cumulative cases.

knitr::include_app("https://bgautijonsson.shinyapps.io/Hierarchical_Report_NegBin/", height = "600px")

Some things to keep in mind

Comparisons between countries

This model does not predict the total amount of COVID-19 cases, it predicts the total amount of DIAGNOSED COVID-19 cases. It then depends on each country’s specific diagnostic criteria what that number means. Iceland does many tests per capita so there it is predicting a broader spectrum of cases, not all serious, whereas in a country like Italy or Spain the model’s predictions have another meaning as the diagnostic criteria are different from those of Iceland.

Changes in policies

As of yet, the model only implements one slope for each country, which is assumed to be fixed throughout the epidemic. Each country’s policies should of course affect this growth rate and the date at which a country implements a policy should affect when the slope changes. Finding solutions to this is on our next actions list.

Next actions

Definite

  • Time-varying growth rates.
  • Include country-level information in the model

Maybe

  • Include information on number of recovered to estimate country-specific recovery rates conditional on being diagnosed.
  • Include information on deaths to estimate country-specific death rates conditional on being diagnosed.

Appendix

Stan Code

data {
  int<lower = 0> N_obs;
  int country[N_obs];
  vector[N_obs] days;
  int new_cases[N_obs];
  
  int<lower = 0> N_countries;
  int<lower = 0> total_cases[N_countries];
  int<lower = 0> total_deaths[N_countries];
  vector[N_countries] pop;
}

parameters {
  vector<lower = 0>[N_countries] beta;
  vector[N_countries] alpha;
  vector<lower = 0, upper = 1>[N_countries] S;
  
  real<lower = 0> mu_beta;
  real<lower = 0> sigma_beta;
  
  real mu_alpha;
  real<lower = 0> sigma_alpha;
  
  real<lower = 0> a_s;
  real<lower = 0> b_s;
  
  vector<lower = 0>[N_countries] phi_inv_sqrt;
  real<lower = 0> sigma_phi_inv_sqrt;
  
}

transformed parameters {
  vector[N_obs] linear = alpha[country] + beta[country] .* days;
  vector<lower = 0, upper = 1>[N_obs] difference;
  vector<lower = 0>[N_countries] phi_inv = square(phi_inv_sqrt);
  vector<lower = 0>[N_countries] phi = inv(phi_inv);
  for (i in 1:N_obs) {
    difference[i] = beta[country[i]] * S[country[i]] * exp(-linear[i]) / square(exp(-linear[i]) + 1);
  }
}

model {
  
  mu_alpha ~ normal(-3.6, 0.5);
  sigma_alpha ~ exponential(1);
  alpha ~ normal(mu_alpha, sigma_alpha);
  
  mu_beta ~ normal(0.13, 0.05);
  sigma_beta ~ exponential(2);
  beta ~ normal(mu_beta, sigma_beta);
  
  S ~ beta(a_s, b_s);
  
  a_s ~ cauchy(0, 5);
  b_s ~ cauchy(0, 5);
  
  phi_inv_sqrt ~ normal(0, sigma_phi_inv_sqrt);
  sigma_phi_inv_sqrt ~ exponential(1);
  
  new_cases ~ neg_binomial_2(difference .* pop[country], phi[country]);
}