karimn / covid-19-transmission

0 stars 0 forks source link

Denmark (a singleton) with no hierarchical parameters fails catastrophically #21

Open karimn opened 4 years ago

karimn commented 4 years ago

This is a good case to investigate to identify non-identification or weak identification problems. It's relatively small and quick to run, and has no sub-regions so we can run out the hierarchical parameters as the cause of divergence.

karimn commented 4 years ago

@wwiecek I'm now investigating this since it is relatively easy to play around with. Take a look at the report for this country using the typical configuration of priors and trend. mobility_report.pdf

Two things to notice: the pairs() plot at the end shows some bimodal distributions that I think are good signs of degeneracy in the model, and the "new cases" plot showing the new cases data (which isn't used in the model vs the predicted level). The median predicted level seems outrageously high.

I'm going to try to run with even tighter imputed cases. Let me know if you have any other ideas.

karimn commented 4 years ago

No improvement with tighter prior on imputed cases.

wwiecek commented 4 years ago

This new report template is amazing and very informative -- thanks!

At a glance, everything looks good. If you didn't tell me that the model is struggling, I don't know if it would be obvious to me! Only when looking at the posterior I can see that some parameters are multimodal, implying chains getting stuck rather than true multimodality. Is it 4 chains that we are looking at?

There is no "deeper" problem with the median cases to me: some of the crazier upper ranges are due to IFR getting stuck at zero, but this could be because of problems elsewhere. You can expect 5-10k new infections every day for this level of mortality. Median is random because you can have 2 chains stuck with high IFR, 2 chains with sensible IFR. For New Cases visual, maybe you can multiply the observed by max_predicted/max_observed? Just to align the two shapes. (Observed < predicted by definition due to testing.)

What's up with this perk in late February:

image

Is it because you have t* in late February? So that a lot of slopes crowd around that point?

I'd start by fixing log(imputed_cases[1]). Seems like the model wants either exp(4.5) or exp(6.5). Maybe we could see if there is still a bi-modal behaviour if we fix it at 100 cases? (=exp(4.6)). I'm not saying that imputed cases are at fault here, but (literally) taking some paramters out of the equation may be a way to investigate. Or maybe even fixing imputed cases AND ifr_noise == 0? What do you think?

karimn commented 4 years ago

I found an insanely simple bug in the code that seems to have been causing all the problem. Makes my head explode! I start getting suspicious by the very high number of divergences. This has happened to me before. When almost all iterations are divergent there is something wrong with how one of the parameters is coded.

wwiecek commented 4 years ago

Can you ref the commit here? Just curious about the diff. What I understood was that you had a "free" parameter defined in the model just floating around? Happened to me with my meta-analysis package when enabling/disabling various hierarchies. But then I'd get DTs but no convergence problems. Maybe it's an issue here because our model is so complex.

Glad this works now!

karimn commented 4 years ago

The merge is referenced above, #22.

Yes, a floating parameter that is not used anywhere. It doesn't make a lot of sense to me. I would imagine such a parameter to just float freely and not affect the model at all. But I've been bit by this before so I should know better!

Anyway, this doesn't solve the problem 100%. I occasionally still see some high Rhat/low ESS (divergence can just been taken care by raising adapt-delta high enough, so it looks like false positives). What does see so far to get rid of this problem is removing the IFR noise while leaving imputed cases as they were (not fixed).

karimn commented 4 years ago

Here's the report for a successful run mobility_report.pdf. Notice that the cases look like they're a much more reasonable level.

wwiecek commented 4 years ago

The merge is referenced above, #22.

Thanks & sorry -- confused myself with too many tabs at once. This looks exactly like the problem I had with hierarchical models & conditional variable definitions/dimensions in the past. I wonder if there is a theoretical explanation why having a parameter just "floating free" like that can destroy the model.

Two thoughts:

1) Fixing IFR is not a bad price to pay for convergence. We can always try investigating IFR a bit differently.

2) The report looks really good now. Note how this model is typically bad at having non-smooth predicted cases (e.g. the perk in March followed by a drop):

image

This is a feature of all simple epi models, nothing we can do about it. Of course there should be some mismatch over time because early on you'd have little testing. So the ratio of observed to predicted will keep getting higher over time as countries get better at testing.

BTW I thought the peak of observed would match the peak of median? But any rescaling is fine, this is mostly to be able to compare the shapes.

karimn commented 4 years ago

@wwiecek, about the cases plot, the model isn't really using any data on cases it can't respond to peaks in observed data. Maybe I'm missing something here. Also, the new cases out of the model are assumed to be the "true" ones generating observed deaths, while the observed cases can suffer from under-reporting.

wwiecek commented 4 years ago

@wwiecek, about the cases plot, the model isn't really using any data on cases it can't respond to peaks in observed data. Maybe I'm missing something here. Also, the new cases out of the model are assumed to be the "true" ones generating observed deaths, while the observed cases can suffer from under-reporting.

No, it can't respond to the peaks in cases because it doesn't "see" them (and what you say about under-reporting), but I was just pointing out that with these models you really have to have very extreme changes in R to see the shape of cases time series such as in observed data. Our model will always give you an "overly smoothed" version of reality.

karimn commented 4 years ago

Oh, I see what you mean. Yes, I agree and that will be the problem with any model we use. If I saw otherwise I would conclude that we're overfitting.