This post is by Bob
There are several reasons I prefer central intervals.
First, you can’t tell how much of the probability is above or below a highest posterior density interval. For example, if I look at the HPD for an exponential distribution, it will run from 0 to the inverse cdf of the interval probability. If I look at the HPD for a normal, it matches the central interval.
Second, the range is going to be parameterization dependent. If I parameterize a probability in terms of log odds in (-infinity, infinity) vs. a log probability in (-infinity, 0) versus a probability (0, 1), I get different HPDs. Similarly, you’ll get different intervals for standard deviations vs. variance (even accounting for the squaring). In contrast, quantile-based intervals like central intervals are invariant to changes of variables.
Third, you don’t even get an interval if the highest density regions are disconnected, as in something like a Beta(0.5, 0.5) distribution. There, the highest density regions are from 0 up and from 1 down, so any highest density region is going to exclude the central interval of lowest density.
There is a lot of confusion about probabilities vs. densities in the wild. For instance, here’s a document from RStudio: Highest Posterior Density (HPD) Intervals. The document defines quantile-based (e.g., central) credible intervals in the usual way. Then it makes the common mistake of conflating density and probability, with the claim, “However, this definition does not guarantee the interval contains the highest probability regions.” Both the 90% HPD and 90% central interval contain exactly the same amount of probability mass, namely 90%. What this means is that the probability that a random draw falls in either is the same—90%. What the authors presumably meant is that central intervals don’t necessarily contain the highest density points. I think this confusion between probability and density is a lot of what’s motivating HPD. And also, the highest density points are going to depend on the parameterization (my second objection above).
On the other hand, Andrew’s the one who was pushing for replacing central intervals with HPD intervals in Stan, so I’m curious what his motivation for that change was given that he’s definitely not confused about probability versus density.
Bob:
I have a few things to say on this:
1. Usually I’m fine with central intervals.
2. If the marginal posterior distribution has two large and widely separated modes, then, yeah, the highest probability density region is two disconnected intervals. In some settings, that’s awkward; in other settings, that’s actually a good summary, if the data and model really do show that a certain parameter could be near A or near B but not in between.
3. If the marginal posterior distribution has two large and widely separated modes, then I doubt that a central probability interval is a good summary anyway. It’s easy to compute, and we compute it in Stan, fine, but it’s not much of a summary of what’s going on.
4. If the marginal posterior density is unimodal, then the “highest posterior intervals” are also “shortest posterior intervals.” I prefer the “shortest” formulation as that seems like a better motivation for me. For the reasons you say, I don’t particularly care about the posterior density function, but there is an appeal in taking the shortest interval. Yes, this depends on parameterization, but that’s fine: if there’s a desire for a shortest interval, that will indeed depend on how the problem is parameterized.
5. For that matter, we’re talking about marginal summaries, and any of these marginal intervals will depend on parameterization once the model has multiple parameters.
6. Another virtue of the “shortest posterior interval” formulation is that it immediately suggests a computational approach based on taking the shortest interval including the specified percentage of simulation draws. The Spin paper (with Ying Liu and Tian Zheng) uses this as a starting point. Ultimately I’m not thrilled with the algorithm in that paper—it’s just too hacky for my taste—but the general idea seems reasonable to me.
7. Finally, my motivating example: what led me to work on this problem in the first place. It’s the marginal posterior distribution for the hierarchical scale parameter in the 8 schools problem. Here it is, directly from BDA:
The central 95% interval for tau is awkward, because its left endpoint comes close to zero but doesn’t quite get there. For this particular problem, if I need to give an interval summary, I’d rather give the shortest probability interval, which will clearly have a left endpoint at zero.
If we switch to log(tau), it’s another story. For log(tau), the central interval seems to me to be a reasonable summary, and the shortest interval will of course not go to -∞. And that’s fine with me. Another way of saying it is that, in this example, the data are consistent with values of tau all the way down to zero, and tau has a flat prior, so the 95% shortest probability interval for tau includes zero. It’s also true that the data are consistent with values of log(tau) all the way down to -∞, but the prior for log(tau) rules out extremely low values, hence the 95% shortest probability interval for log(tau) does not go to -∞. The dependence of the shortest probability interval on the parameterization is directly connected to the dependence of the prior on the parameterization.
8. To return to point 1, usually I’m fine with central intervals. There are cases where I don’t like the central interval, as in the 8 schools example, but even there it’s not horrible. And there are bimodal posteriors where I don’t like any interval summary.
9. Stepping back for a moment, we can ask: Why summarize using intervals at all? And I don’t have any great answer to that question! As I’ve written elsewhere, I think that confidence intervals would be better termed “uncertainty intervals.”
Can uncertainty about a fitted model be usefully summarized by supplying an interval for each parameter in the model? Sometimes, sometimes not. I think it’s good that Stan supplies these intervals; we should just recognize that they don’t always work. No two-number summary does.
When it comes to scale parameters I usually utilize a Gamma that goes to zero density at x=0 it rarely makes sense to say a distribution is a delta function as implied by scale = 0 so I like my priors to avoid zero. I also like to use dimensionless parameters where ive already incorporated a factor that brings my parameter to scale around 1, So something like Gamma(2,1) or Gamma(3,1/2) is a reasonable first pass for a scale parameter.
Daniel:
Yes, we propose zero-avoiding gamma priors in our 2013 paper, whose results I announced in this 2011 post.
Thanks, Andrew. My complaint was the opposite of (7), that the HPD interval is awkward in these cases because as a summary, I usually want to know not just an interval but how much probability is above and below it.
And I completely agree the bimodal case is problematic—a single interval summary isn’t good.
I. think the right question is (9), why summarize with intervals at all? I think the answer is that it’s easy to write in a paper, whereas plots take up a lot of space and mental energy trying to figure out what their shapes mean.
I also agree that (9) is a problem—the whole idea of summarizing with marginal intervals for parameters is a mess if we think about inference, when correlation among parameters becomes important. I find this often clouds discussions of identifiability in Bayesian statistics.
I really like HPD for summarizing (marginal) bivariate posteriors, eg using the ggdensity package: Interpretable Bivariate Density Visualization with ggplot2 • ggdensity https://share.google/lhnHrhP7dcteMSRLb
Not sure there is a good “central region” analog here?
Thanks for bringing up the multivariate case. We should probably be paying more attention to this.
We used a notion of points whose log densities were in the central alpha interval (alpha in (0, 1)):
I(Theta, alpha) = { theta : log p(theta) in central alpha interval of log p(Theta)}.
We introduced this because it feels better than a typical set, which doesn’t have specified coverage. It’s defined for epsilon > 0 as:
T(Theta, epsilon) = { theta : abs(E[log p(Theta)] – log p(theta)) < epsilon}
> We used a notion of points whose log densities were in the central alpha interval (alpha in (0, 1))
That seems a lot of work to avoid the obvious high-density region.
For a two-dimentional standard Gaussian the 95% interval will be a circular ring between r=0.22 and r=2.72.
I don’t think that it’s any better than the circular region with r=2.45. However, if skipping the high probability (sorry, high density) center that accounts for 2.5% of the mass is desirable you may want to do the same in the one-dimensional case. Forget the central interval [-1.96 1.96] and use [-2.24 -0.03] ∪ [0.03 2.24] instead!
> First, you can’t tell how much of the probability is above or below a highest posterior density interval.
Someone could prefer left tails for the same reason. Why do you care about the center anyway?
> Second, the range is going to be parameterization dependent. […] In contrast, quantile-based intervals like central intervals are invariant to changes of variables.
To some changes of variables at least.
> Third, you don’t even get an interval if the highest density regions are disconnected, as in something like a Beta(0.5, 0.5) distribution. There, the highest density regions are from 0 up and from 1 down, so any highest density region is going to exclude the central interval of lowest density.
That seems a good thing? It may be better to localize 80% of the distribution in [0 0.35] and [0.65 1] rather than in [0.02 0.98]. It also seems that for an exponential [0 1.6] may be a better 80% interval than [0.1 2.3].
> Both the 90% HPD and 90% central interval contain exactly the same amount of probability mass, namely 90%. What this means is that the probability that a random draw falls in either is the same—90%.
I would expect that most people have the right intuition about all 90% probability regions containing a similar amount of probability mass. I don’t think that their use of “high probability region” is particularly confusing. (For example, you wrote some time ago that “the typical set *is* the high probability *mass* set by definition” and that was much more confusing! https://statmodeling.stat.columbia.edu/2017/03/15/ensemble-methods-doomed-fail-high-dimensions/ )
Anyway, one of the reasons to prefer high density intervals is that the generalization to high density regions is trivial.
It’s not that I care about the center per se. If all I was doing was testing calibration it wouldn’t matter. But I thought the intervals were supposed to be summaries. My main gripe is that with an HPD interval, you just don’t know how much is above or below the interval you specified. That feels like a serious deficiency for a two-number summary. That is, you can give me a 90% HPD interval and I don’t know if 10% fall above or below or a mix. So let’s say I give you a survival interval for a patient as an HPD. They don’t know how much probability there is of dying before the interval starts or surviving after the interval ends. You get both from central intervals.
I shouldn’t have brought up multimodal as central intervals are also terrible summaries there. I don’t know if it’s any worse, but a lot of people seem to think so.
I think people will have the right intuition that each 90% interval has the same probability. At least people who haven’t been drilled into focusing on interval tail statistics to calculate p-values. I think the write-up from RStudio was confusing probability mass and density.
The discussion of typical sets has been very confusing over the course of Stan because Michael Betancourt was using the term informally. Sorry for contributing to that! Technically, an epsilon-typical set is the points whose densities fall within epsilon of the expected log density (and it’s used in a concentration example in increasing dimensions, to make matters even more confused). Michael was using it to mean the part of the probability mass we sample with MCMC. With MCMC, that tends to be central values because the deep tails are rare by construction (i.e., they are atypical in the sense of “typical” intended for “typical set”).
> But I thought the intervals were supposed to be summaries.
You can also have summaries that won’t tell you how much is above or below some level, right?
As Ken Rice mentioned you can define different intervals based on different loss functions (minimising the regret which is the “expected excess loss” – pretty similar by the way to the “average excess loss” used in other contexts).
You can also have single-value “summaries” which don’t have to tell you how much is above or below that level. A median (or a tail quantile) may tell you that but a mean or a mode won’t. They have other desirable properties though.
> The discussion of typical sets has been very confusing […] Sorry for contributing to that!
I think the way it’s discussed in the new book is an improvement. Still, we’re told it’s important but it’s not clear why. We get an example where the so-called typical set covers some points where the probability is several quadrillion times higher than in other points (sorry, points where the probability density is several quadrillion times higher than in other points). That’s a big epsilon! (We don’t want our simulations to be in the areas of parameter space where 𝑙𝑝 log 𝑝(𝜃|𝑦) is near its expectation. We want them to be everywhere, according to the probability density of each area.)
Reading a bit more to see if the concept was relevant at some point I didn’t found anything, but I noticed the following: “The Monte Carlo standard error
of the average of 10 independent draws of 𝜃 is 1/√10 times the standard deviation of 𝜃in the target distribution; thus inference about the posterior mean based on 10 simulations is √1.1 1.05 times as uncertain as from infinite draws.”
Maybe I misunderstand the reference to “inference about the posterior mean” but with infinite draws the value of the posterior mean is estimated with infinite precision. I think that you mean that drawing from a perfectly defined posterior distribution there is some amount of uncertainty and drawing for a not-so-perfectly defined posterior distribution introduces just a bit of incremental uncertainty, but referring to that as “inference about the posterior mean” seems misleading.
Carlos:
Good catch. We’ll change “inference about the posterior mean” in that sentence on page 198 to “inference about theta.”
With an unlimited number of ways to construct a posterior interval/region with a given level of support, it really helps to be clear about what criteria we’re using to pick between them. One way of formalizing this is here, where a “regret” argument motivates central and HPD intervals, and some other well-known choices.
Thinking about the problem in this way may also be a way to wean people off defaulting to 95 percent coverage, which is almost as deeply-ingrained as alpha=0.05, for no compelling reason of which I’m aware. Working backward from a given coverage level, the corresponding “regret” level tells you how bad an estimate (relative to the best possible choice) you are willing to tolerate. Regret is measured on the same scale as theta itself, which may be easier to connect with the underlying science than percentage coverage.
I haven’t explored it, but you could even have the whole discussion on the scale of theta, picking an interval based on the regret one can tolerate. The posterior would still quantify overall uncertainty, but while the resulting interval’s coverage could be evaluated it wouldn’t be fixed at any specific value.
That seems like a differential notion of “regret” than the standard one used in reinforcement learning. The paper’s paywalled, so I could only read the abstract (blame the rent-seekers at the American Statistical Association).
I would have thought that thinking probabilistically, any interval with the same coverage would be interchangeable.
Free full text is here. Reinforcement learning’s not my area but I’m aware it gets used elsewhere, I hope the “regret” terms are at least close.
… I’m aware “regret” gets used elsewhere, I hope the terms are at least close
Years ago I did investigate reparameterization invariant intervals and regions based on maximizing predictive utility (minimizing predictive regret?) https://arxiv.org/abs/2109.10843. These are theoretically nice, but eventually I concluded that the practical value is low as most of the time we shouldn’t report just one interval or if we do report it doesn’t matter much how we construct that interval, and then it’s better to use easy to compute approaches. It’s rare that we would need to report an interval and the decision based on that intervals would be such that it would be rally important to get that interval to have specific properties (like maximizing utility).
we have been using credibility intervals, see for example 3.8.3 https://www.springerprofessional.de/en/modern-statistics/23516892