Investigate why model fails when adding Poland to a number of smaller countries

[karimn]

A run including AR AU CA PT runs fine and a separate run for PL also runs fine, but combining results in high Rhat and very low ESS measures. @wwiecek suspects that the root of the problem is the beta coefficients of the mobility model.

[karimn]

Just ran with no hierarchical structure to the mobility model (no beta subnational and national effects) and the problem still fails. The beta coefficients seem fine in this run.

[karimn]

Running with no hierarchical structure at all, R0 and the beta parameters, works.

[karimn]

Allowing hierarchical effects for the beta parameters but fully pooling the R0 still works. It is looking more likely that the hierarchical nature of the R0 parameters is the problem.

[karimn]

Averaging all the mobility variables as

(- g_residential + g_transit_stations + g_grocery_and_pharmacy + g_parks + g_retail_and_recreation + g_workplaces) / 6

Model still has problems:

Sampling: 9801.823 sec elapsed Warning messages: 1: The largest R-hat is NA, indicating chains have not mixed. Running the chains for more iterations may help. See http://mc-stan.org/misc/warnings.html#r-hat 2: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. Running the chains for more iterations may help. See http://mc-stan.org/misc/warnings.html#bulk-ess 3: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. Running the chains for more iterations may help. See http://mc-stan.org/misc/warnings.html#tail-ess

Extracting results...

Maximum Rhat = 1.531146 Minimum ESS Bulk = 7.146713 Minimum ESS Tail = 26.46728

[karimn]

others_beta.txt poland_beta.txt

I'm attaching posterior distributions of the beta parameters for Poland and, separately, the other countries. One thing that jumps out is the second coefficient (for g_transit_stations) is consistently positive and with greater magnitude.

[karimn]

Below is a look at how R0 differs.

Poland:

> pl_results %>% unnest(param_results) %>% filter(parameter == "R0") %>% select(sub_region, starts_with("per_"), mean)
# A tibble: 12 x 11
   sub_region        per_0.025 per_0.05 per_0.1 per_0.25 per_0.5 per_0.75 per_0.9 per_0.95 per_0.975  mean
   <chr>                 <dbl>    <dbl>   <dbl>    <dbl>   <dbl>    <dbl>   <dbl>    <dbl>     <dbl> <dbl>
 1 Kuiavia-Pomerania      1.99     2.08    2.17     2.33    2.49     2.64    2.80     2.90      2.99  2.49
 2 Lublin                 1.93     2.05    2.16     2.32    2.50     2.65    2.82     2.94      3.04  2.49
 3 Łódź                   2.05     2.15    2.21     2.34    2.51     2.67    2.81     2.91      3.01  2.51
 4 Lesser Poland          2.08     2.18    2.28     2.42    2.59     2.77    2.97     3.10      3.27  2.61
 5 Mazovia                2.24     2.31    2.37     2.51    2.68     2.89    3.13     3.29      3.42  2.72
 6 Opole                  2.02     2.13    2.20     2.35    2.51     2.67    2.83     2.94      3.06  2.52
 7 Subcarpathia           2.06     2.15    2.25     2.39    2.55     2.73    2.89     3.02      3.11  2.56
 8 Podlaskie              1.75     1.86    1.97     2.17    2.36     2.53    2.67     2.76      2.84  2.34
 9 Silesia                2.04     2.11    2.20     2.32    2.47     2.60    2.75     2.83      2.90  2.47
10 Swietokrzyskie         1.86     1.97    2.09     2.28    2.45     2.62    2.77     2.86      2.95  2.44
11 Greater Poland         2.13     2.22    2.30     2.42    2.58     2.74    2.91     3.01      3.10  2.59
12 West Pomeranian        1.75     1.86    1.98     2.19    2.39     2.57    2.71     2.81      2.90  2.37

Other countries:

> incr_results %>% unnest(param_results) %>% filter(parameter == "R0") %>% select(sub_region, starts_with("per_"), mean)
# A tibble: 14 x 11
   sub_region           per_0.025 per_0.05 per_0.1 per_0.25 per_0.5 per_0.75 per_0.9 per_0.95 per_0.975  mean
   <chr>                    <dbl>    <dbl>   <dbl>    <dbl>   <dbl>    <dbl>   <dbl>    <dbl>     <dbl> <dbl>
 1 City of Buenos Aires     2.38      2.57    2.82     3.28    3.86     4.58    5.38     5.89      6.31  4.00
 2 Chaco                    2.36      2.53    2.79     3.25    3.84     4.53    5.28     5.88      6.36  3.96
 3 Río Negro                2.21      2.40    2.67     3.16    3.75     4.45    5.20     5.75      6.32  3.87
 4 Córdoba                  2.40      2.59    2.85     3.30    3.89     4.60    5.38     5.95      6.50  4.03
 5 New South Wales          1.17      1.23    1.30     1.43    1.61     1.80    1.98     2.10      2.19  1.63
 6 Tasmania                 1.08      1.15    1.24     1.40    1.60     1.84    2.07     2.23      2.39  1.64
 7 Victoria                 0.938     1.02    1.10     1.26    1.45     1.67    1.88     1.99      2.10  1.48
 8 Alberta                  1.19      1.21    1.24     1.29    1.35     1.41    1.46     1.49      1.52  1.35
 9 British Columbia         1.24      1.26    1.28     1.33    1.37     1.42    1.47     1.50      1.52  1.38
10 Nova Scotia              1.24      1.30    1.38     1.51    1.66     1.84    2.03     2.17      2.29  1.69
11 Ontario                  1.73      1.76    1.80     1.86    1.93     2.01    2.08     2.12      2.16  1.94
12 Quebec                   2.29      2.36    2.44     2.62    2.84     3.12    3.45     3.68      3.88  2.90
13 Saskatchewan             1.16      1.27    1.38     1.59    1.86     2.15    2.48     2.70      2.91  1.90
14 Azores                   1.19      1.34    1.55     1.90    2.36     2.85    3.42     3.79      4.15  2.44

[wwiecek]

That's not the way we would want to average over mobility. Sorry, totally my fault as I didn't say what I've meant. Residential always goes up when other metrics go down. We have to average over the other 5 dimensions or just the 3 dimensions that were in the Vollmer model.

[karimn]

No problem, re-running with newer specification. However, the wrong averaging one have resulted in a single variable that even if measuring the wrong thing would not suffer from correlation. So eliminating correlation does not seem to do it. What do you think?

Also, holding R0 fixed and allowing betas to vary over national/subnational units sometimes works. Above I reported it works, but after several runs, I noticed that it does not consistently work.

[wwiecek]

I think you're right. Residential is only one of 6 variables, so it shouldn't have so much of an impact. In any case I think it's better to proceed with just 1 average mobility variable for now, to make that part of the model simpler. What do you think?

So we know that the no pooling model (without hierarchical structure) runs and the model with full pooling (completely fixed R0) does too. (Do I remember this correctly?) Do we know if a model with very large SD on both R0 and betas runs? When SDs are sufficiently high, it should revert to no pooling case.

[karimn]

Attempted three scenarios:

1. More diffuse priors for both R0 and betas, did not result in convergence. I will try with greater diffusion and perhaps leaving the betas at their default level.

   
   hyperparam_tau_beta_toplevel: 0.5
   hyperparam_tau_beta_national_sd: 1 
   hyperparam_tau_beta_subnational_sd: 0.25

hyperparam_toplevel_R0_sd: 0.25 hyperparam_tau_national_effect_log_R0_sd: 0.5 hyperparam_tau_subnational_effect_log_R0_sd: 0.08

2. A single averaged mobility variable of the five non-residential variables, did not result in convergence.

3. More concentrated priors for R0 and beta resulted in convergence. I will try loosening a bit and see what happens.

hyperparam_tau_beta_toplevel: 0.5 hyperparam_tau_beta_national_sd: 0.0025 hyperparam_tau_beta_subnational_sd: 0.25

hyperparam_toplevel_R0_sd: 0.25 hyperparam_tau_national_effect_log_R0_sd: 0.0008 hyperparam_tau_subnational_effect_log_R0_sd: 0.08

FYI, 
a) To run with a specific hyperparameters for R0 and the betas use a YAML file (as above) and pass in using the `--hyperparam=<file> option.`
b) To specify a particular formula for the design matrix use the `--mobility-model=<formula>` option. For example, for (2) above I used `--mobility-model=~0+average_all_mob`

[karimn]

Running with a more relaxed prior of betas (but not R0) converges nicely.

Hyperparameters:

hyperparam_tau_beta_toplevel: 0.5
hyperparam_tau_beta_national_sd: 0.25
hyperparam_tau_beta_subnational_sd: 0.125

hyperparam_toplevel_R0_sd: 0.25
hyperparam_tau_national_effect_log_R0_sd: 0.0008
hyperparam_tau_subnational_effect_log_R0_sd: 0.08

[wwiecek]

Running with a more relaxed prior of betas (but not R0) converges nicely.

Hyperparameters:

hyperparam_tau_beta_toplevel: 0.5
hyperparam_tau_beta_national_sd: 0.25
hyperparam_tau_beta_subnational_sd: 0.125

hyperparam_toplevel_R0_sd: 0.25
hyperparam_tau_national_effect_log_R0_sd: 0.0008
hyperparam_tau_subnational_effect_log_R0_sd: 0.08

That's a relief, because the other day I tried doing it and it didn't work, which was very suprising. Probably I did something wrong!

[karimn]

@wwiecek I ran the top 10 in terms of death singleton countries (with a single sub-region) and I'm not getting any decent convergence. High Rhat and low ESS. I ran with the same settings as the separate countries run.

I'm seeing similar results with some of the other non-singleton like Canada.

[wwiecek]

Could it be that it fails because of some recent updates? Or hypervariances parameters are to blame, because when I sequentially, one-by-one, ran the "top" 15 countries (wrt mortality), only Ireland and Turkey struggled. I did it about ~10 days ago.

[karimn]

That's probably it: the hyper variance parameters. If you ran one-by-one you a non-pooling model.

On Thu, Jul 9, 2020 at 6:03 AM Witold Wiecek notifications@github.com wrote:

Could it be that it fails because of some recent updates? Or hypervariances parameters are to blame, because when I sequentially, one-by-one, ran the "top" 15 countries (wrt mortality), only Ireland and Turkey struggled. I did it about ~10 days ago.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/karimn/covid-19-transmission/issues/9#issuecomment-656035351, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABB5525SX4XIP66SVMICUTR2WIWJANCNFSM4N7SVJ5Q .

[karimn]

Running Poland by itself with "--no-pooling" turned on does not fix things; three chains finished while one of them is still at 10%.

[wwiecek]

I pasted outputs from PL in imputed cases thread, I think. You will see that with better imputed cases prior it was only the one region that was failing, number 4 I think. Can you look at Rhat for imputed cases and/or R0's and check whether it's all failing or just a single region? If it's the latter, we should just troubleshoot by re-running and tweaking for that 1 region repeatedly

[karimn]

Last run of countries separated PL ran without a problem. Need to check what's going on when run with others.