Ok gotcha.
Good, although I would be worried if there are any genuine tail ESS warnings given that your report has all these 95% credible intervals whose endpoints are not estimated very precisely if the tail ESS is low.
1 Like
Me too! Is there an equivalent to the code that @saudiwin posted above to check for tail ESS warnings after the fact?
Just look at monitor(fit)
@bgoodri Any thoughts on the banana pairs plot above? Does that geometry simply mean sampling will be less efficient, or is it a concern for inference more generally?
If it gets more extreme, then there could be divergent transitions. But if there are no divergent transitions for now, then Stan is presumably managing it with some cost to the runtime.
1 Like
Thanks both. monitor(fit) for what was in our report last week looks pretty good to me–4000 samples is giving > 3000 ESS (bulk and tail). Rhats all very very close to 1.0. I’ve seen guidance somewhere about numbers of effective samples to give stable estimates of 95% CIs but I don’t have the reference handy.
Note that monitor(fit) didn’t actually work because there are almost 6K parameters including all the generated quantities stuff so when I print that out it’s overwhelming. I did this, and assume it’s the same thing, wonder if there’s a better way:
`sims = as.array(fit)
rhats=vapply(1:dim(sims)[3], function(x) Rhat(sims[,x]), 0)
range(rhats[(rhats)])
`
result=monitor(fit,print=F) looks good, but doesn’t report Rhat?
You can do
monitor(extract(fit, permuted = FALSE, inc_warmup = TRUE,
pars = c("mu", "E_deaths"))
to only compute it for a subset of pars and then Rhat is in the third column for the right. The best reference is
Rstan is relatively slow with so many outputs. You should consider to calculate the gq in a second go (maybe with the gq facility). That should speed up things. Another option is to use cmdstan, which does not suffer from this bottleneck as I recall.
1 Like
On a separate note, I’ve been working with a team collecting a massive amount of this kind of intervention data (every country in the world, currently at 5,000 records). As part of that I’ve designed a simpler model that does something similar to your’s except without modeling the infectious disease process. Basically I use prior information about the ratio between COVID-19 tests and the unobserved number of infected to identify the latent attack rate. Paper here:
https://osf.io/preprints/socarxiv/jp4wk/
It seems like the model I designed has a relatively tight CI around the attack rate while yours models deaths with more certainty and the count of infected with more uncertainty. Anyways appreciate any comments on the approach! It looks like the challenge with scaling your’s up further is obtaining estimates of the age-stratified IFR (though this model is not really to our purpose as we aren’t trying to model the epidemic).
1 Like
As quantile estimators are essentially based on the difference of the probabilities on either side of the quantile value you want about 100 effective samples on each side of the quantile. For a 95% interval you would probably look at 2.5% and 97.5% quantiles, which by definition mean that 2.5% of the effective samples have to be on one side of each boundary quantile so you’d need at least 4000 effective samples total. You’re close, but reporting a 90% interval would be much more robust.
Regarding the model a two small recommendations: see if non-centering mu speeds things up and consider fitting 1 / phi so that you can implement a prior that concentrates around a Poisson instead of one that concentrates around infinite overdispersion. For detail on the latter see Probabilistic Building Blocks.
6 Likes
Can you provide a reference for this? Google isn’t turning anything up.
Just an update on effective tools for debugging–I had a different model with some extra parameters that have the same problem. So I was unsurprised to see these errors:
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
But when I checked monitor the offending parameters had Rhat 1, Bulk ESS 4000:
result=monitor(fit,print=F)
as.data.frame(result)[2518,c("Bulk_ESS","Tail_ESS","Rhat")]
Bulk_ESS Tail_ESS Rhat
Rt_adj[1,1] 4000 4000 1
To find these parameters I used:
sims = as.array(fit)
rhats=vapply(1:dim(sims)[3], function(x) Rhat(sims[,,x]), 0)
2 Likes
So it looks like what you need to do is pass a data value turning off gq in the main fit, and then a data value turning them back on when you want gqs?
Sorry, there was a bug in monitor which was fixed Jan 10 2020 in monitor.R rstan repo. @flaxter, if you can, you could use the that fixed monitor.R.
Instead of thinking in terms of ESS, you can also think directly in terms of Monte Carlo error for the interval end points. Reference: Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner (2019). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. arXiv preprint arXiv:1903.08008. By default monitor() function is computing also MCSE_Q2.5, MCSE_Q25, MCSE_Q50, MCSE_Q75, and MCSE_Q97.5 even if it’s not printing those. Just do mon <- monitor(fit) and then mon has those columns, too. If you need other quantiles, you can use probs argument.
EDIT: typo fix
1 Like
I think you have the limits of the uniform backwards; it should be:
E_deaths0[1, m] += uniform_rng(1e-16, 1e-15);
Right?
Yes, check for a comprehensive set of changes pushed yesterday to our github. Still more fixes to make as suggested by Bob, but this fixes the infinite Rhats!
1 Like
Here’s the link: https://github.com/ImperialCollegeLondon/covid19model/blob/master/stan-models/base.stan
1 Like