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-04-14

2020-04-02

  • Model the \(\beta_i\) as lognormal since they are strictly positive.
  • There is enough data to loosen the informative priors, and we have done so.
  • Added posterior predictions for effective reproduction number, \(R_e\), along with some text about it.
  • Reparametrised to use a noncentered parametrisation to aid Stan in performing NUTS.

2020-03-31

  • Reparametrised the prior for S in terms of \(\mu\) the mean and \(\kappa\) the sample size.

2020-03-28

  • Weak priors on \(a_S\) and \(b_S\) to help the NUTS sampler.
  • Rewrote priors for scale parameters.

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

Parametrisation

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 parametrisation 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 model \(\beta_i\) on the \(\log\) scale to impose the constraints \(\beta_i > 0\). We put hierarchical priors on \(\alpha_i\) and \(\beta_i\) so that for each country, \(i\),

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

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

\[ \begin{aligned} \mu_\alpha &\sim \mathrm{Normal}(-2.5 , 3^2) \\ \mu_\beta &\sim \mathrm{Normal}(-3, 1^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 we parametrise in terms of the mean, \(\mu_S\), and prior sample size, \(\kappa_S\)

\[ \begin{aligned} \mu_S &\sim \mathrm{Beta}(1, 99) \\ \kappa_S &\sim \mathrm{Exponential}(0.001). \end{aligned} \]

We can then back-transform into \(a_S\) and \(b_S\) as

\[ \begin{aligned} a_S &= \mu_S \cdot \kappa_S \\ b_S &= (1 - \mu_S) \cdot \kappa_S \end{aligned} \]

We parametrise 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{Exponential}(\sigma_\phi) \\ \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{LogNormal}(\mu_\beta, \sigma^2_\beta) \\ \alpha_i &\sim \mathrm{Normal}(\mu_\alpha, \sigma^2_\alpha) \\ \mu_\beta &\sim \mathrm{Normal}(-3, 1^2) \\ \mu_\alpha &\sim \mathrm{Normal}(-2.5, 3^2) \\ p(\sigma_\alpha) &\sim \mathrm{Exponential}(1) \\ p(\sigma_\beta) &\sim \mathrm{Exponential}(1) \\ \sqrt\phi_i &\sim \mathrm{Exponential}(\sigma_\phi^2)\\ \sigma_\phi &\sim \mathrm{Exponential}(1) \\ S_i &\sim \mathrm{Beta}(a_S, b_S) \\ a_S &= \mu_S \cdot \kappa_S \\ b_S &= (1 - \mu_S) \cdot \kappa_S \\ \mu_S &\sim \mathrm{Beta}(1, 99) \\ \kappa_S &\sim \mathrm{Exponential}(0.001) \end{aligned} \]

Data

We use Icelandic data from the Directorate of Health and supplement with data on other countries from the European CDC, obtained daily.

Data used for modeling were filtered according to

  • Only use data for countries with a population above 200.000.
  • Only keep data for which cases per 1000 inhabitants \(\geq\) 0.1.
  • Only keep data for countries for which the first observation after this filtering is \(\leq\) 0.2
  • Only keep data for countries with more than 6 days of follow-up after the above steps.
  • Countries that are not included for other reasons
    • China
    • South Korea
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-04-04 11 0.106 0.162
Armenia 2020-03-27 19 0.111 0.351
Australia 2020-03-26 20 0.111 0.253
Austria 2020-03-17 29 0.113 1.568
Bahrain 2020-03-14 32 0.128 0.829
Barbados 2020-03-30 16 0.115 0.251
Belgium 2020-03-18 28 0.108 2.651
Bosnia And Herzegovina 2020-03-31 15 0.107 0.311
Canada 2020-03-27 19 0.107 0.686
Chile 2020-03-29 17 0.101 0.397
Croatia 2020-03-26 20 0.101 0.399
Cyprus 2020-03-25 21 0.105 0.561
Denmark 2020-03-13 33 0.117 1.095
Djibouti 2020-04-08 7 0.124 0.306
Dominican Republic 2020-04-01 14 0.103 0.295
Ecuador 2020-03-29 17 0.106 0.433
Estonia 2020-03-16 30 0.129 1.005
Finland 2020-03-23 23 0.113 0.554
France 2020-03-17 29 0.102 1.506
French Polynesia 2020-03-27 19 0.107 0.197
Germany 2020-03-20 26 0.169 1.498
Greece 2020-03-29 17 0.101 0.205
Iceland 2020-03-05 39 0.109 5.047
Iran 2020-03-12 34 0.109 0.884
Ireland 2020-03-20 26 0.114 2.181
Israel 2020-03-22 24 0.104 1.360
Italy 2020-03-09 37 0.122 2.634
Kuwait 2020-04-05 10 0.114 0.309
Latvia 2020-03-25 21 0.103 0.344
Lithuania 2020-03-27 19 0.108 0.388
Luxembourg 2020-03-16 30 0.125 5.347
Malaysia 2020-04-04 11 0.104 0.151
Malta 2020-03-19 27 0.109 0.872
Mauritius 2020-03-31 15 0.107 0.270
Moldova 2020-04-02 13 0.105 0.423
Montenegro 2020-03-27 19 0.107 0.436
Netherlands 2020-03-19 27 0.120 1.553
New Zealand 2020-03-30 16 0.115 0.224
North Macedonia 2020-03-28 18 0.105 0.410
Norway 2020-03-13 33 0.115 1.206
Panama 2020-03-25 21 0.104 0.818
Poland 2020-04-06 9 0.108 0.183
Portugal 2020-03-22 24 0.125 1.656
Puerto Rico 2020-04-03 12 0.108 0.308
Qatar 2020-03-14 32 0.113 1.141
Romania 2020-03-31 15 0.101 0.343
Serbia 2020-04-01 14 0.103 0.414
Singapore 2020-03-27 19 0.102 0.503
Slovakia 2020-04-08 7 0.106 0.141
Slovenia 2020-03-16 30 0.105 0.583
Spain 2020-03-15 31 0.123 3.627
Sweden 2020-03-16 30 0.103 1.091
Switzerland 2020-03-14 32 0.130 2.968
Turkey 2020-03-30 16 0.110 0.732
United Arab Emirates 2020-04-03 12 0.105 0.463
United Kingdom 2020-03-25 21 0.120 1.312
United States 2020-03-23 23 0.107 1.770
Uruguay 2020-04-03 12 0.107 0.151 
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, a total of 4000 iterations. The slowest chain converged in 388 seconds.
  • Adapt_delta was set to 0.99

Warnings

  • There were 7 divergent transitions after warmup.
  • There were 2 transitions after warmup that exceeded maximum treedepth.

We have examined the model carefully to find what causes these divergent transitions, and our current belief is country-specific correlations between \(\beta_i\) and \(S_i\) when the model predicts a high probability of the country being close to its maximum. We have now removed South Korea from the data as the model has a hard time fitting the daily number of new cases in a country after it has reached its asymptote. We are working on steps to overcome this flaw.

Country Level Effects

results %>% 
  ungroup %>% 
  filter(par %in% c("beta", "alpha", "S", "phi_inv", "nu")) %>%
  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
S Albania 0.0003 0.0002 0.0016 1.0005 3278
Armenia 0.0006 0.0004 0.0015 1.0008 2393
Australia 0.0004 0.0002 0.0017 1.0009 3349
Austria 0.0018 0.0016 0.0021 0.9997 6207
Bahrain 0.0018 0.0010 0.0042 1.0005 5004
Barbados 0.0005 0.0002 0.0021 1.0003 4048
Belgium 0.0035 0.0030 0.0044 1.0001 6182
Bosnia And Herzegovina 0.0006 0.0004 0.0014 1.0019 2102
Canada 0.0010 0.0008 0.0014 1.0010 3050
Chile 0.0008 0.0005 0.0016 1.0015 2413
Croatia 0.0007 0.0005 0.0013 1.0014 1958
Cyprus 0.0007 0.0005 0.0015 1.0001 2344
Denmark 0.0018 0.0013 0.0034 1.0002 3144
Djibouti 0.0017 0.0007 0.0046 0.9998 7187
Dominican Republic 0.0008 0.0005 0.0028 1.0003 2980
Ecuador 0.0011 0.0005 0.0034 0.9999 4580
Estonia 0.0014 0.0010 0.0022 1.0002 4190
Finland 0.0010 0.0006 0.0025 1.0014 2564
France 0.0019 0.0016 0.0023 1.0005 5189
French Polynesia 0.0003 0.0002 0.0019 0.9999 3123
Germany 0.0022 0.0019 0.0029 0.9998 5262
Greece 0.0003 0.0002 0.0011 1.0021 2123
Iceland 0.0052 0.0046 0.0059 0.9999 12848
Iran 0.0012 0.0009 0.0021 1.0005 3308
Ireland 0.0054 0.0036 0.0091 1.0001 6899
Israel 0.0020 0.0015 0.0029 0.9999 4363
Italy 0.0033 0.0030 0.0037 0.9997 8613
Kuwait 0.0010 0.0005 0.0026 1.0008 3080
Latvia 0.0005 0.0003 0.0014 1.0009 2369
Lithuania 0.0007 0.0004 0.0021 1.0023 2994
Luxembourg 0.0059 0.0048 0.0075 1.0004 7079
Malaysia 0.0002 0.0001 0.0008 1.0021 1759
Malta 0.0015 0.0010 0.0032 1.0006 3322
Mauritius 0.0006 0.0003 0.0019 1.0006 3511
Moldova 0.0009 0.0005 0.0022 1.0025 2419
Montenegro 0.0007 0.0004 0.0016 1.0008 3204
Netherlands 0.0025 0.0021 0.0033 1.0019 1321
New Zealand 0.0005 0.0002 0.0018 1.0026 3512
North Macedonia 0.0010 0.0005 0.0030 1.0001 4356
Norway 0.0014 0.0012 0.0019 1.0000 4107
Panama 0.0017 0.0011 0.0034 0.9998 3451
Poland 0.0004 0.0002 0.0013 1.0021 1769
Portugal 0.0024 0.0019 0.0034 1.0022 4309
Puerto Rico 0.0009 0.0004 0.0029 0.9998 4308
Qatar 0.0023 0.0013 0.0051 1.0001 5724
Romania 0.0007 0.0005 0.0021 1.0002 2857
Serbia 0.0011 0.0005 0.0032 1.0004 3366
Singapore 0.0014 0.0006 0.0038 0.9999 5636
Slovakia 0.0005 0.0002 0.0026 1.0003 5212
Slovenia 0.0007 0.0006 0.0009 0.9997 2651
Spain 0.0042 0.0039 0.0045 1.0000 11303
Sweden 0.0016 0.0013 0.0025 1.0002 4052
Switzerland 0.0036 0.0028 0.0047 1.0007 6171
Turkey 0.0020 0.0013 0.0040 1.0004 4118
United Arab Emirates 0.0016 0.0007 0.0042 1.0004 5241
United Kingdom 0.0025 0.0019 0.0045 1.0010 3489
United States 0.0025 0.0023 0.0029 0.9998 6755
Uruguay 0.0002 0.0001 0.0015 1.0013 2452
Beta Albania 0.1080 0.0550 0.1717 1.0003 4323
Armenia 0.1139 0.0595 0.1788 1.0009 3730
Australia 0.1590 0.1142 0.2176 1.0003 4536
Austria 0.1657 0.1311 0.1958 0.9999 5990
Bahrain 0.0925 0.0609 0.1299 1.0002 4567
Barbados 0.1155 0.0601 0.1927 0.9998 4393
Belgium 0.1562 0.1230 0.1846 0.9999 7126
Bosnia And Herzegovina 0.1130 0.0589 0.1823 1.0004 3669
Canada 0.1538 0.1072 0.1879 1.0006 4475
Chile 0.1085 0.0612 0.1550 1.0009 3275
Croatia 0.0997 0.0537 0.1458 1.0007 2932
Cyprus 0.1491 0.0823 0.2183 1.0001 4396
Denmark 0.0981 0.0687 0.1289 1.0000 4359
Djibouti 0.1246 0.0724 0.2076 1.0002 8513
Dominican Republic 0.0975 0.0550 0.1525 1.0007 3282
Ecuador 0.1172 0.0653 0.1871 0.9997 8009
Estonia 0.1141 0.0697 0.1593 0.9999 5039
Finland 0.1043 0.0592 0.1602 1.0004 4092
France 0.1440 0.1098 0.1731 1.0000 6135
French Polynesia 0.1053 0.0526 0.1694 1.0005 3739
Germany 0.1172 0.0843 0.1483 0.9999 5952
Greece 0.1144 0.0645 0.1791 1.0010 3540
Iceland 0.1761 0.1524 0.1976 0.9998 8462
Iran 0.1049 0.0646 0.1451 1.0003 4055
Ireland 0.1085 0.0858 0.1362 1.0000 8975
Israel 0.1329 0.0856 0.1830 0.9999 6538
Italy 0.1132 0.0972 0.1280 0.9996 8721
Kuwait 0.1113 0.0594 0.1810 1.0004 4367
Latvia 0.1160 0.0643 0.1820 0.9999 4227
Lithuania 0.1174 0.0615 0.1929 1.0014 4641
Luxembourg 0.1700 0.1210 0.2171 1.0000 5135
Malaysia 0.1001 0.0499 0.1573 1.0009 2931
Malta 0.0990 0.0580 0.1455 1.0000 3824
Mauritius 0.1224 0.0646 0.2197 0.9998 5789
Moldova 0.1239 0.0650 0.2182 1.0011 4882
Montenegro 0.1398 0.0743 0.2500 1.0001 4767
Netherlands 0.1105 0.0836 0.1346 1.0000 4548
New Zealand 0.1378 0.0805 0.2222 1.0003 6892
North Macedonia 0.1116 0.0626 0.1783 1.0001 5144
Norway 0.1305 0.0887 0.1670 0.9997 4970
Panama 0.1055 0.0651 0.1522 0.9999 4258
Poland 0.1122 0.0589 0.1844 1.0002 3313
Portugal 0.1362 0.0905 0.1814 1.0016 5561
Puerto Rico 0.1111 0.0616 0.1784 1.0002 4823
Qatar 0.1062 0.0756 0.1409 1.0000 7407
Romania 0.1072 0.0584 0.1700 0.9997 3890
Serbia 0.1171 0.0627 0.1939 0.9999 5518
Singapore 0.1152 0.0697 0.1726 0.9998 7010
Slovakia 0.1106 0.0560 0.1857 0.9997 5483
Slovenia 0.1194 0.0814 0.1509 0.9999 3800
Spain 0.1649 0.1476 0.1804 1.0001 8830
Sweden 0.1270 0.0988 0.1531 0.9997 7216
Switzerland 0.1363 0.0943 0.1780 0.9998 6763
Turkey 0.1228 0.0900 0.1609 1.0006 6738
United Arab Emirates 0.1186 0.0680 0.1931 0.9998 8642
United Kingdom 0.1294 0.0930 0.1674 0.9998 6391
United States 0.1568 0.1361 0.1751 0.9998 8161
Uruguay 0.1071 0.0526 0.1794 0.9999 3577
Alpha Albania 0.0212 -2.9751 2.2339 1.0024 2695
Armenia -0.7226 -2.0032 0.7906 1.0018 5052
Australia 0.5378 -0.6015 2.2163 1.0006 3448
Austria -1.8257 -2.3121 -1.2184 0.9998 6306
Bahrain -3.1192 -4.1161 -2.1459 0.9998 9602
Barbados -0.5299 -2.9991 1.6497 1.0002 3479
Belgium -3.0561 -3.4386 -2.6216 0.9997 8549
Bosnia And Herzegovina -0.6707 -1.8380 0.5950 0.9998 5124
Canada -1.8513 -2.1459 -1.4875 1.0001 7122
Chile -1.3507 -1.8037 -0.8728 0.9995 7306
Croatia -0.7633 -1.3949 0.0947 1.0002 4665
Cyprus -2.0114 -2.7202 -1.2155 1.0001 8065
Denmark -2.7255 -3.2518 -2.1975 0.9998 9419
Djibouti -1.8387 -3.3482 0.1490 0.9999 6895
Dominican Republic -1.5882 -2.9868 -0.5000 1.0002 5474
Ecuador -2.2731 -3.7243 -0.6551 1.0000 6938
Estonia -1.8065 -2.5224 -0.9301 1.0001 7007
Finland -1.7791 -2.8631 -0.7011 1.0000 6588
France -2.5255 -2.9380 -2.0351 0.9997 7520
French Polynesia -1.1476 -3.6423 1.4159 1.0008 3892
Germany -1.2714 -1.7274 -0.7471 0.9998 6778
Greece 0.0229 -0.9941 1.7671 1.0017 3030
Iceland -3.8031 -4.2669 -3.3143 1.0000 8480
Iran -2.3076 -3.0241 -1.5106 0.9998 7489
Ireland -3.1006 -3.5803 -2.6177 1.0001 9512
Israel -1.7796 -2.3970 -1.0765 1.0000 8134
Italy -2.3708 -2.6599 -2.0685 0.9996 9062
Kuwait -0.5374 -2.3279 1.0360 1.0002 4736
Latvia -0.5511 -1.7179 1.0485 1.0008 3703
Lithuania -0.7776 -2.4927 1.0708 1.0007 4662
Luxembourg -2.1820 -2.9349 -1.3408 0.9999 5850
Malaysia -0.1058 -1.5546 1.4126 1.0017 2636
Malta -2.0028 -2.8426 -1.1010 0.9998 7449
Mauritius -0.4458 -2.6722 1.5530 1.0014 3956
Moldova -1.0193 -1.9813 -0.0174 0.9998 7338
Montenegro -0.8595 -2.0550 0.6337 0.9998 4830
Netherlands -2.1186 -2.4086 -1.8033 0.9998 8218
New Zealand 0.6978 -0.8583 2.4153 1.0019 3485
North Macedonia -2.0322 -3.5614 -0.4178 1.0001 6735
Norway -1.9365 -2.5894 -1.1495 0.9998 5340
Panama -1.8687 -2.4948 -1.2829 0.9998 9039
Poland 0.1181 -1.2729 1.6168 1.0031 2089
Portugal -1.9801 -2.5129 -1.3904 1.0006 7032
Puerto Rico -1.4369 -3.1323 0.4097 0.9999 7036
Qatar -3.6285 -4.6062 -2.6295 1.0003 7947
Romania -1.4619 -2.5006 -0.6761 0.9997 6759
Serbia -1.2217 -2.8989 0.5733 1.0002 6173
Singapore -2.7081 -4.0789 -1.2047 1.0005 7015
Slovakia -1.3974 -3.7577 2.2537 1.0001 2823
Slovenia -1.5846 -2.1160 -0.9247 0.9996 4216
Spain -2.9230 -3.1736 -2.6499 0.9998 8898
Sweden -3.1729 -3.5025 -2.8207 1.0000 9404
Switzerland -2.1894 -2.9339 -1.3934 0.9999 7413
Turkey -2.2422 -2.7397 -1.9548 0.9999 5481
United Arab Emirates -1.5082 -3.0952 0.3842 1.0000 5989
United Kingdom -2.4576 -2.8390 -2.0819 1.0004 7902
United States -2.5087 -2.6876 -2.3042 0.9998 9397
Uruguay -0.3874 -3.3562 2.1710 1.0028 2314
Phi Albania 0.1965 0.0420 0.7034 1.0004 5175
Armenia 0.4448 0.2182 0.9592 1.0005 8396
Australia 0.1274 0.0667 0.2768 0.9998 8946
Austria 0.0929 0.0551 0.1672 1.0001 7963
Bahrain 0.8760 0.5357 1.4871 0.9998 9357
Barbados 0.8805 0.2181 2.5657 0.9999 9388
Belgium 0.0844 0.0503 0.1564 0.9998 7755
Bosnia And Herzegovina 0.1333 0.0561 0.3463 1.0003 7437
Canada 0.0203 0.0102 0.0478 0.9999 5871
Chile 0.0298 0.0140 0.0721 0.9999 7408
Croatia 0.0676 0.0300 0.1635 0.9997 8421
Cyprus 0.1058 0.0399 0.2597 0.9998 7126
Denmark 0.1965 0.1215 0.3354 1.0001 8991
Djibouti 0.5523 0.2009 1.7102 1.0001 6643
Dominican Republic 0.1381 0.0628 0.3442 0.9998 6674
Ecuador 0.8048 0.4614 1.5333 1.0000 9045
Estonia 0.3654 0.2193 0.6429 0.9998 8855
Finland 0.4744 0.2611 0.9104 1.0000 10008
France 0.0967 0.0579 0.1758 0.9999 7741
French Polynesia 0.8518 0.0967 2.6644 1.0000 8837
Germany 0.0682 0.0404 0.1270 1.0001 8225
Greece 0.1341 0.0629 0.3198 0.9996 7589
Iceland 0.1171 0.0676 0.2103 1.0007 8291
Iran 0.4439 0.2743 0.7419 0.9998 8068
Ireland 0.1051 0.0586 0.2077 1.0002 7558
Israel 0.1591 0.0913 0.3023 1.0003 7924
Italy 0.0468 0.0297 0.0797 1.0001 8140
Kuwait 0.1259 0.0476 0.3924 0.9999 5336
Latvia 0.4672 0.2355 0.9802 0.9999 8493
Lithuania 0.9137 0.4656 1.8776 1.0007 8920
Luxembourg 0.3263 0.1919 0.5744 1.0000 8675
Malaysia 0.0414 0.0139 0.1500 1.0010 3684
Malta 0.3840 0.2066 0.7410 0.9998 10421
Mauritius 0.7002 0.2995 1.6635 1.0000 7183
Moldova 0.0767 0.0294 0.2157 0.9998 6572
Montenegro 0.3352 0.1284 0.8105 1.0000 7407
Netherlands 0.0352 0.0205 0.0672 0.9997 7922
New Zealand 0.4056 0.1702 1.0355 1.0000 6386
North Macedonia 0.7750 0.3887 1.6516 0.9999 8833
Norway 0.2791 0.1656 0.5051 1.0002 8130
Panama 0.1009 0.0541 0.2057 1.0005 8939
Poland 0.0253 0.0075 0.1219 1.0013 2479
Portugal 0.1080 0.0618 0.2094 1.0005 5922
Puerto Rico 0.3627 0.1645 0.8971 1.0002 7861
Qatar 0.9327 0.5890 1.5304 0.9996 8720
Romania 0.0830 0.0396 0.2090 1.0000 7123
Serbia 0.6694 0.3192 1.5242 1.0001 9484
Singapore 1.1054 0.5899 2.1577 1.0004 8652
Slovakia 0.6610 0.2317 1.9882 1.0000 7210
Slovenia 0.0807 0.0360 0.1687 1.0001 9020
Spain 0.0340 0.0207 0.0604 0.9998 8411
Sweden 0.0723 0.0418 0.1328 0.9999 8924
Switzerland 0.3621 0.2174 0.6461 1.0000 8717
Turkey 0.0133 0.0065 0.0323 1.0007 6647
United Arab Emirates 0.8113 0.3765 1.9294 1.0008 7742
United Kingdom 0.0467 0.0261 0.0953 1.0000 7154
United States 0.0130 0.0073 0.0266 1.0000 6943
Uruguay 0.2582 0.0832 0.7877 1.0006 6415 
p <- results %>% 
  ungroup %>% 
  filter(par %in% c("beta", "alpha", "S", "phi_inv", "nu")) %>% 
  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_beta", 
                    "mu_alpha", "sigma_alpha", 
                    "mu_s", "kappa_s",
                    "sigma_phi_inv_sqrt", "mu_nu", "sigma_nu")) %>% 
  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 -2.1191 -2.3060 -2.0002 1.0012 1218
sigma_beta 0.2582 0.1579 0.4232 1.0029 1117
mu_alpha -1.6156 -1.9800 -1.2044 1.0038 1454
sigma_alpha 1.1992 0.9136 1.6315 1.0071 1268
mu_s 0.0017 0.0013 0.0022 1.0007 2644
kappa_s 985.3293 616.0072 1505.5048 1.0000 5989
sigma_phi_sqrt 0.5115 0.3977 0.6733 1.0072 668

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.

Posterior predictions for \(R_e\) are obtained by looking at part of the SIR differential equation on daily new cases

\[ \frac{d}{dt} I = \beta \cdot I \cdot \frac{S}{N} - \text{Some rate of recovery} \]

We do not model the rate of recovery, but based on literature review we assume that diagnosed cases are healthy after 15 days. As we are not interested in the negative part of \(\frac{d}{dt}I\) we write \(\frac{d}{dt} I_{i, t} = D_{i, t}\), the number of new cases. We can then simplify to estimate \(\beta\) in country \(i\) at time \(t\) as

\[ \beta_{i, t} = D_{i, t} \cdot \frac{N}{I_{i, t} \cdot S_{i, t}} \]

where \(I_{i, t}\) is the cumulative number of diagnosed cases at time \(t\) minus the cumulative number of cases at time \(t - 15\). \(S_{i, t}\) is simply the population minus the cumulative number of cases at time \(t\). Then we utilize

\[ R_e = \beta \gamma^{-1} \]

where \(\gamma^{-1}\) is the length of the serial interval. After a literature review we put a \(\mathrm{Gamma}(30, 2.5)\) distibution on the length of serial interval.

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

Validation

GIFs with historical predictions for some countries made by the model can be seen here.

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

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

Estimating the burden on the Icelandic health care system

We use the current age distribution of cases in Iceland to split our posterior predictive distritubion of daily cases into age categories. We then use numbers from Ferguson et al. at Imperial College, (Table 1), to sample numbers of patients in hospital and from those numbers, how many will end up in ICU. Along with the numbers above, we use the following parameters in our simulations

  • lag_infected_to_hospital <- 7
  • lag_hospital_to_icu <- 3
  • days_from_infection_to_healthy <- 21
  • days_in_hospital <- 14
  • days_in_icu <- 10

The results from our predictions are then shared at covid.hi.is

Appendix

Stan Code

data {
  int<lower = 0> N_obs;
  // country[N_obs] takes on value i for country number i etc.
  int country[N_obs];
  // days is how many days since case rate per 1000 exceeded a limit we have chosen
  vector[N_obs] days;
  // We model the daily number of newly diagnosed cases instead of cumulative numbers
  int new_cases[N_obs];
  
  int<lower = 0> N_countries;
  vector[N_countries] pop;
}

parameters {
  // Since we use a non-centered parametrisation we first create normal(0, 1) variables for alpha and beta
  
  // Beta parameters
  vector[N_countries] z_beta;
  real mu_beta;
  real<lower = 0> sigma_beta;
  
  // Alpha parameters
  vector[N_countries] z_alpha;
  real mu_alpha;
  real<lower = 0> sigma_alpha;
  
  //  Asymptote/saturation parameters
  // We model the beta distribution in terms of mean and sample size instead of a and b.
  vector<lower = 0, upper = 1>[N_countries] S;
  real<lower = 0, upper = 1> mu_s;
  real<lower = 0> kappa_s;
  
  //  Overdispersion parameters. One for each country which is dependent on hyperparameter
  vector<lower = 0>[N_countries] z_phi_inv_sqrt;
  real<lower = 0> sigma_phi_inv_sqrt;
  
}

transformed parameters {
   // Non-Centerd parametrisations
   // If B ~ normal(mu_b, sigma_b) then B = mu_b + sigma_b * normal(0, 1)
  vector<lower = 0>[N_countries] beta = exp(mu_beta + sigma_beta * z_beta);
  vector[N_countries] alpha = mu_alpha + sigma_alpha * z_alpha;
   // If X ~ exponential(lambda) then X ~ lambda * exponential(1)
  vector<lower = 0>[N_countries] phi_inv_sqrt = sigma_phi_inv_sqrt * z_phi_inv_sqrt;
   // Overdispersion parameters
  vector<lower = 0>[N_countries] phi_inv = square(phi_inv_sqrt);
  vector<lower = 0>[N_countries] phi = inv(phi_inv);
   // Asymptote hyperparameters
  real<lower = 0> a_s = mu_s * kappa_s;
  real<lower = 0> b_s = (1 - mu_s) * kappa_s;
   // Logistic equation calculations
  vector[N_obs] linear = alpha[country] + beta[country] .* days;
  vector<lower = 0>[N_obs] difference;
  for (i in 1:N_obs) {
    difference[i] = beta[country[i]] * S[country[i]] * exp(-linear[i]) / square(exp(-linear[i]) + 1);
  }
}

model {
   // Alpha parameters
  mu_alpha ~ normal(-2.5, 3);
  sigma_alpha ~ exponential(1);
  z_alpha ~ std_normal();
  
   // Beta parameters
  mu_beta ~ normal(-3, 1);
  sigma_beta ~ exponential(1);
  z_beta ~ std_normal();
  
   // Asymptote parameters
  mu_s ~ beta(1, 99);
  kappa_s ~ exponential(0.001);
  S ~ beta(a_s, b_s);
  
   // Overdispersion parameters
  z_phi_inv_sqrt ~ exponential(1);
  sigma_phi_inv_sqrt ~ exponential(1);
  
  //  Likelihood
  new_cases ~ neg_binomial_2(difference .* pop[country], phi[country]);
}