# Bayesian approach to report simulation studies?

I am running a simulation study where we want to estimate a proportion $$p$$. We are reporting the coverage of credible intervals with a uniform prior, and we are doing $$500$$ Monte Carlo simulations. We report the proportion of intervals that contain the true value of the parameter, and report a confidence interval for this proportion. I was wondering if I should also do a Bayesian analysis of the results from the simulation study or if the frequentist nature of the simulation study favours reporting confidence intervals.

Q. More in general, is it valid to use a Bayesian approach to report a simulation study?

• If your credible interval is computed without any approximations, then its average coverage will be equal to the nominal level if you take the average with respect to the your prior for $p$. So simulations may not be needed (if you think the average coverage in this sense is what is relevant). – Jarle Tufto Dec 11 '18 at 16:01
• But there is no true value of $p$, $p$ is a random variable. How can you "compare" the coverage probability for the credible interval where $p$ is random to the coverage probability for the confidence interval where $p$ is fixed? – AdamO Dec 11 '18 at 16:13
• @AdamO You compute the frequentist coverage for different values of $p$ and then compute the mean coverage with respect to the prior distribution of $p$ – Jarle Tufto Dec 11 '18 at 16:31
• @AdamO That brings in the interpretation of credible intervals, as a quantification of uncertainty. This is more philosophical, but one can do coverage analysis in Bayes as well. – Donald Dec 11 '18 at 16:41
• @AdamO What I mean by average coverage is $\int c(p)\pi(p)dp$ where $c(p)$ is frequentist coverage for a given value of $p$ and $\pi(p)$ is the prior. – Jarle Tufto Dec 11 '18 at 17:23

The following is a simple proof of the (perhaps obvious) fact that a Bayesian credible interval (or set) always has average frequentist coverage equal to the nominal level (assuming a proper prior). Some authors (e.g. Brown et al. 2001 and Fagerland, Lydersen & Laake 2017, section 1.4.2), in choosing between the large number of approximate frequentist confidence intervals, argue that such average coverage may be more relevant than ensuring that the coverage at least has the nominal level for all $$\theta$$. If one accepts this, an average weighted by prior belief in different values of $$\theta$$ is arguably a sensible way to assess overall coverage (as done below).

By definition, a Bayesian $$(1-\alpha)$$-credible interval (or set) $$C(x)$$ (some function of the data $$x$$) includes $$\theta$$ with probability $$P(\theta \in C(x)|x)=1-\alpha. \tag{1}$$ conditional on the data $$x$$. Thinking of both $$\theta$$ and $$x$$ as random variables, it is also true, by the law of total probablity, that $$P(\theta \in C(x))=\int P(\theta \in C(x)|x)f(x)dx=1-\alpha \tag{2}$$ where $$f(x)$$ is the marginal density of $$x$$.

For the frequentist coverage of the Bayesian credible interval, we instead condition on only the parameter $$\theta$$ such that $$x$$ and $$C(x)$$ are random. This gives the frequentist coverage $$P(\theta \in C(x)|\theta) \tag{3}$$ which in general will depend on $$\theta$$. But the average frequentist coverage, weighted by the prior $$\pi(\theta)$$, is $$\int P(\theta \in C(x)|\theta) \pi(\theta) d\theta=P(\theta \in C(x)) =1-\alpha, \tag{4}$$ by the law of total probability and (2), and so always perfectly matches the nominal level.

The following plot and R-code illustrates (3) and (4) for a binomial model with a uniform prior on $$p$$ and $$n=20$$:

R code

m <- 1e+5
p <- seq(0,1,len=m)
coverage <- numeric(m)
alpha <- .05
a <- 1
b <- 1
n <- 20
x <- 0:n # possible outcomes
lower <- qbeta(alpha/2, a+x, b+n-x) # lower and upper credible limits
upper <- qbeta(1-alpha/2, a+x, b+n-x) # for each value of x
for (i in 1:length(p)) {
px <- dbinom(x,size=n,prob=p[i]) # probability of possible outcomes given p
coverage[i] <- sum(px*ifelse(lower<=p[i] & p[i]<=upper, 1,0)) # coverage
}
plot(p,coverage,type="l",ylim=c(.85,1))
abline(h=1-alpha, lty=3)

> mean(coverage)
[1] 0.9499928


This is probably a bit opinion based, but I would normally expect that many of the usual arguments for using Bayesian methods may not apply. E.g. bringing in extra information to compensate for low sample size (usually one can just simulate more) and wanting to use prior information (firstly, exact prior information may be hard to come but and it may be less controversial to just simulate more) would not seem like strong arguments.

On the other hand, using a Bayesian method to get an easy to compute approximation to (exact) frequentist confidence intervals using percentiles of the beta distribution may be a reasonable idea - as long as you have enough simulations that your exact choice of prior will be uncontroversial. Putting too strong a prior on a desirable result (e.g. using a Beta(95,5) prior for the coverage probability of 95% credible intervals) might end up being hard to defend, while I would guess most people would not have much of a problem using vague priors such as e.g. Beta(1/3, 1/3), Beta(0.5, 0.5) or the uniform Beta(1, 1) prior, if you have a decent number of simulations.

The obvious frequentist alternative is something like an exact Clopper-Pearson confidence interval. I have never had anyone questioning (e.g. on the basis of methodological incoherence) me using one those for simulation studies, even when I was evaluating Bayesian methods. In a way, you are after all evaluating the frequentist properties of a method (repeat sampling performance). Around that I dimly recall some quote (perhaps by Rubin?) that you use the Bayesian approach to construct methods and evaluate these properties in a frequentist way (and that they usually turn out to have good frequentist properties).

On the specifics of your project, 500 simulations seems like a pretty low number for the coverage probability of credible intervals (unless they are e.g. 50% intervals, certainly quite low for 95% ones). Even if you observe exactly 95% of simulations with the CI covering your true value, your credible interval for the credible interval coverage under a Beta(0.5, 0.5) prior would be from 92.8% to 96.7% (2500 simulations would get you to 94.1% to 95.8% and 10,000 to 94.6 to 95.4%). I guess if that level of uncertainty is fine in your case that may not be an issue.

• Thanks, very interesting points. Of course, the last paragraph should account for the running times of the simulation. As if it requires 1 hour per iteration, you would think twice before running 10,000 simulations, and that level of accuracy may not be that bad. Are you the famous Referee that always want more simulations? (tongue in cheek). – Donald Dec 11 '18 at 15:42
• As a referree I have done that (although most often it was "that sounds nice in theory, but can you run at least some simulations to show that it does in finite sample size what you say it does."). Admittedly, I have access to a high performance cluster where I can run a couple of hundred jobs in parallel for a few weeks so that we can try out relatively complex Bayesian models using a decent number of simulation. – Björn Dec 11 '18 at 17:51