Does coverage matter?

In response to Andrew Gelman’s extended April Fool’s diatribe on Objections to Bayesian Statistics, Larry Wasserman commented regarding physicists who want  guaranteed frequentist coverage for their confidence intervals that  “Their desire for frequentist coverage seems well justified. Someday, we can count how many of their intervals trapped the true parameter values and assess the coverage. The 95 percent frequentist intervals will live up to their advertised coverage claims. A trail of Bayesian intervals will, in general, not have this property”.

One thing to note about this statement is that it’s just not true.  Confidence intervals produced in actual scientific research are notorious for not covering the true value, even when they are produced using frequentist recipes.  This is why high-energy physicists insist on such absurdly high confidence levels (or absurdly low p-values) before declaring discoveries — what they call “five sigma” evidence, which corresponds to a p-value of less than 10^-6. If taken seriously, quoting such a small p-value would be pointless, since any reader would surely assign a higher probability than that to the possibility that the “discovery” results from fraud or gross incompetence. The high confidence levels demanded are just an ad hoc way of trying to compensate for possible inadequacies in the statistical model used, which can easily make the true coverage probability be much less than advertised (or the true Type I error rate much higher than advertised).

Let’s ignore this, though, since discussions of theory omitting messy practical issues can be valuable.  The next thing to ask, then, is how it is possible that a 90% Bayesian probability interval — which purports to contain the true value with 90% probability —  can contain the true value less than 90% of the time.  A simple example will show how this can happen, and provide insight into whether we should care.

One can trivially get an example of low or zero coverage using a Bayesian method in which some parameter values are assigned low or zero prior probability, but that’s not very interesting, since if you use such a prior, you presumably want such values to be considered low probability (and hence excluded from posterior intervals) even if they’re not especially disfavoured by the data.  But here I’ll give an example in which the coverage is zero for a parameter value that has just as high prior probability as all the other parameter values.

Suppose we have an unknown parameter, θ, with ten possible values, 0, 1, …, 9.  We obtain data, x, which  has nine possible values, 1, …, 9.  If θ=0, the observation x is equally likely to be any of the values 1, …, 9 —  ie, the probabiity of  each of these values is 1/9.  If θ>0, we will observe x=θ with probability one.  Suppose we use a uniform prior for θ, so each of the values 0, 1, … 9 has prior probability 1/10.

Once we observe x, only two values of θ will remain possible — θ=x and θ=0.  By assumption, the ratio of prior probabilities for these two values is one.  The ratio of likelihoods is 1 over 1/9, or 9, in favour of θ=x.  The posterior odds are the product of the prior odds and the likelihood ratio, so we find that the posterior probabilities are 9/10 for θ=x, 1/10 for θ=0, and zero for any other value of θ.  A 90% Bayesian probability region is therefore easy to define — it’s just the set consisting only of the observed value of x.

What’s the coverage of this Bayesian probability region?  Well, if θ>0, the coverage is 100%, since when θ>0 we are guaranteed to observe x=θ.  But when θ=0, the coverage is 0%, since we never produce a posterior probability region that includes 0.  Frequentist coverage is the minimum probability, for any true θ, that the region will include the true θ.  So the coverage for these Bayesian probability regions is zero.

Should zero coverage be cause for worry?  That depends first of all on whether we actually believe the model and the prior that were used to obtain the Bayesian regions, and secondly on whether a posterior probability region actually provides the information that we need.  It’s certainly possible that the answers to both these questions are “yes”, in which case the zero coverage should not be a cause for worry.  But I think that often Bayesian posterior regions are not what we want.

Recall that one way of looking at a 90% frequentist confidence region is as the set of all parameter values for which a hypothesis test would produce a p-value greater than 0.1 — ie, the set of parameter values that are consistent with the data, according to a hypothesis test using a 10% significance level.  From this point of view, the confidence region is a way of telling theorists what theories to discard — namely, all those theories that predict a value for the parameter outside the confidence region.

Does this work for the Bayesian posterior regions in the example above?  If we observe x=9, we construct the 90% posterior region {9}, and conclude that all theories that predict any value for the parameter other than 9 should be discarded (at least if we think 90% is high enough).  This is certainly the right thing to do for theories that predict that the parameter is 1, 2, …, 8, since those values are excluded by the data with certainty.  But what about parameter value 0?  It’s outside the 90% posterior region, but note that its posterior probability of 1/10 is exactly the same as its prior probability.  The observation of x=9 has not reduced the probability that θ=0 at all, so it certainly seems strange to say that θ=0 is now excluded by the data!

I think part of the problem is that reports of experimental results should not be aimed at presenting conclusions, as may seem most natural from a Bayesian viewpoint, but rather at providing the information with which the readers may draw conclusions.  This may be the source of some objections to the prior distribution in Bayesian analysis, which can be seen as corrupting the objective presentation of the experimental results, even though frequentist methods like p-values are not suitable presentations either.  In simple examples like the one above, the experimental results can be communicated fully and objectively (assuming the model is uncontroversial) by reporting the likelihood function — which in this example is 1 for Θ equal to the observed x, 1/9 for Θ equal to 0, and zero for any other Θ.  There is no need for either frequentist confidence regions or Bayesian posterior regions.   Readers can use the likelihood function to produce posterior regions based on whatever priors they like.  (Frequentists using methods that violate the Likelihood Principle are out of luck, but you can’t please everyone.)

In more complex problems with a high-dimensional parameter, however, just presenting the likelihood function is both infeasible and unenlightening.  One solution is to present a marginal likelihood function for just a smaller set of parameters of interest, integrating with respect to a prior distribution for the other “nuisance” parameters.  This only works if the prior for nuisance parameters is relatively uncontroversial, but that may often be the case (or at least, the experimenters may be the best people to formulate this prior).

If you think you must instead produce something like a “Bayesian confidence region”, perhaps the best method is to report the region for which the ratio of posterior probability (or density) to prior probability (or density) is not less than 0.1 (for an analogue of a 90% region).  This can also be thought of as the region containing all  Θ values for which the Bayes Factor is at least 0.1 for a comparison of the model in which this Θ has prior probability one to the model in which Θ has whatever prior you have set up as the “default”.  This can also be seen as the region of Θ where a theorist who assigns prior probability 1/2 to their pet theory predicting that value of Θ would not abandon their theory (using a 90% confidence level) if the they use your “default” prior for parameter values conditional on their theory being false.  For the example above, the regions produced in this way (with uniform default prior) would consist of the observed x plus 0.  Coverage would be 100%.

I hasten to add that this suggestion assumes a context in which various theories predict specific values for the parameter, and our interest is in which of these theories is true.  This is quite unlike many applications of confidence intervals, in which the parameters are things such as the average daily caloric intake of Canadian adults, or the regression coefficient of income on years of education.  For such parameters, the likelihood function (or marginal likelihood function) should be reported.  If this likelihood happens to be approximately normal, and concentrated in a region where any reasonable prior would be nearly uniform, you could go ahead and present a 90% posterior region if you’re so inclined.