How do Bayesians verify their methods using Monte Carlo simulation methods?

Question

Background: I have a PhD in social psychology, where theoretical statistics and math were barely covered in my quantitative coursework. Through undergrad and grad school, I was taught (much like many of you also in the social sciences, probably) through the "classical" frequentist framework. Now, I also love R and using simulation methods to verify that methods work makes way more sense to me than mathematical proofs (again: background in a quantitative social science, not theoretical statistics). Frequentist methods and simulation methods together make a ton of sense to me. Because frequentists see probability as long run odds (e.g., If I do this an arbitrarily large number of times, and it happens 50% of the time, then there is a 50% probability). We can simulate this long run with Monte Carlo methods!

Complications: Since undergrad, I have been very aware of Bayesian methods, and there has always been people in my life calling me to the Bayesian side, saying that the results were easier to interpret, that we get probability for a hypothesis instead of the data given a hypothesis, etc. I was really into this and took a Bayesian class, read some Bayesian books and papers, and now am pretty familiar with Stan and its associated R packages.

Enter Mayo: After thinking "Bayesian is probably the way of the future" for a while, I read Deborah Mayo's Statistical Inference as Severe Testing. She says she doesn't pick a side at the beginning of the book, but she does: She is a frequentist, and a lot of the book is defending frequentist methodologies. I don't want to necessarily get into a discussion of whether or not we think the way she sees evidence is valid, but this got me thinking: Is Bayes really all that it is advertised? I mean, the Bayes crowd is so fractured itself that I don't even know the "right" way to analyze data in a Bayesian framework often. Usually, I would just use rstanarm and present point estimates and credible intervals... which often line up closely with frequentist estimates and confidence intervals. I might do model comparisons, but I'm always afraid of describing Bayes factors as posterior probability comparisons, etc.

More Thinking: What I kept thinking through Mayo's book was: There is a way we can use computers to make sure our frequentist methods work, because probability is what we see in the long run and we can simulate that. Bayesians can't even agree on what probability really is, it seems, depending on the Bayesian school (default, subjective, etc.). Which leads me to my question:

Question: How do Bayesians verify that their methods define uncertainty properly (i.e., calculate valid credible intervals and posterior distributions) using Monte Carlo simulation methods, if probability is not defined as rates in the long run?

Example: I create a data generator. This is just going to simulate from a Bernoulli distribution with a probability of .5:

set.seed(1839)
p <- .50
n <- 100
gen_dat <- function(n, p) {
  rbinom(n, 1, p)
}

Now, let's say I want to make sure that confidence intervals in a logistic regression are actually valid. I can simulate a regression a large number of times and make sure that the actual population value falls in the 95% confidence interval 95% of the time. It's an intercept-only model, so I just want to make sure that I'm estimating p correctly:

set.seed(1839)
iter <- 10000
results <- sapply(seq_len(iter), function(zzz) {
  mod <- glm(gen_dat(n, p) ~ 1, binomial)
  conf <- suppressMessages(confint(mod))
  log(p / (1 - p)) < max(conf) & log(p / (1 - p)) > min(conf)
})
mean(results)

This takes a minutes to run, but we end up with the mean(results) call giving us 0.9416. This is about 95%, and I'm confident in saying that the glm command is describing uncertainty in a valid way. I'm sure it would have gotten closer to right on the nose at 95% if I upped iter and wanted to wait here at my laptop longer.

On the other hand, let's fit a Bayesian model for the same thing:

library(rstanarm)
set.seed(1839)
dat <- data.frame(y = gen_dat(n, p))
stan_mod <- stan_glm(y ~ 1, binomial, dat)
summary(stan_mod)

In part, this gives me:

Estimates:
                mean   sd    2.5%   25%   50%   75%   97.5%
(Intercept)    -0.1    0.2  -0.5   -0.2  -0.1   0.0   0.3  
mean_PPD        0.5    0.1   0.3    0.4   0.5   0.5   0.6  
log-posterior -73.0    0.7 -75.1  -73.1 -72.7 -72.5 -72.5  

Since Bayesians don't define probability as what we see in the long run, how can I use simulation methods to verify than stan_glm is accurately capturing uncertainty? That is, how could I trust that these credible intervals are valid, using simulation methods? And right now, I'm not even defining a prior—how does the inclusion of priors come into play here, since that will affect our measures of uncertainty?

When I was trying to write a beta regression with a hurdle model component in Stan from scratch once, I had someone recommend to me: "Simulate data. Do it a bunch of times, and the true estimates should be in the credible interval about 95% of the time." But to me, that goes against the very thing that Bayesians believe in! That relies on frequentist understandings of probability! So how would a Bayesian convince me that the credible interval I'm getting from the summary() call to my model is accurately describing uncertainty, using simulation methods?

Purpose of Question: This is a trivial example, but many times clients provide me with difficult problems. And I try things that I am unfamiliar with, so I often run a simulation study to make sure that what I am doing is valid. If I were to write a custom model in Stan, how would I know that what I am doing is legit? How could I use simulation methods to verify that what I'm doing in Stan is actually going to tell me what I want to know?

Answer 1

I think you see the logical problem in your question. In the frequentist paradigm, it is fine to presume a population truth, generate data, and see if the estimates have good coverage, because that is what they are suppose to do. In the Bayesian paradigm, however, there is no ground truth to generate data from! Bayesians ask the probability of such truths given data, so in simulation we need different truths that give rise to data and then condition on the data. In practice, one ends up simulating the law of conditional probability, which, fortunately, holds always by definition. I take up this exact issue in Rouder, 2014, Psychonomic Bulletin and Review. https://dx.doi.org/10.3758/s13423-014-0595-4

Answer 2

How do Bayesians verify that their methods define uncertainty properly (i.e., calculate valid credible intervals and posterior distributions) using Monte Carlo simulation methods, if probability is not defined as rates in the long run?

I believe the confusion here is about the purpose of simulation methods in Bayesian statistics. The only purpose of Markov Chain Monte Carlo methods such as Gibbs Sampling or Hamiltonian Monte Carlo is to calculate the denominator of Bayes rule.

Of course, there are often other methods available which would make MCMC needless. Some models can be expressed using conjugacy, others through applying a fine grid over the parameter space, yet others can be solved with acceptance-rejection testing. Where MCMC comes in handy is when the integral is ill-behaved.

While I would love to avoid math, that really cannot be avoided. In looking at Bayes rule $$\pi(\theta|x)=\frac{f(X|\theta)\pi(\theta)}{\int_{\theta\in\Theta}f(X|\theta)\pi(\theta)\mathrm{d}\theta},$$ the numerator is made up of $f(X|\theta)$ and $\pi(\theta)$. $f(X|\theta)$ is a likelihood and not a probability, so it doesn’t sum to one except by chance. The denominator assures us that $\pi(\theta|X)$ sums to one. The goal of MCMC is to determine the bottom number. Note that the bottom number is a constant. It is the expectated likelihood.

The accuracy of that number does determine some but not all parameter estimates. If you were using the maximum a posteriori estimator, then MCMC is an unnecessary step. You should build a hill climbing algorithm instead. On the other hand, it is necessary to determine the posterior mean or an interval. That is because the 95% interval has to be 95% of something and the denominator determines what the scale of that something is.

The goal of MCMC in Bayesian methodologies is to get the Markov chains to converge to the posterior density. That is it. It doesn’t test the validity of anything. It is just an attempt to determine a fixed point value. It is a form of numerical integration. As there is no way to know without letting the algorithm to run to infinity whether all dense regions have been covered, there is some human judgment. The algorithm will have a cutoff when it believes it is done, but that does not mean it is actually done.

In Frequentist methodologies, MCMC is often used to test the reasonableness of a model or to numerically approximate a solution when an analytic one is not available. It serves no similar purpose here.

If I were to write a custom model in Stan, how would I know that what I am doing is legit? How could I use simulation methods to verify that what I'm doing in Stan is actually going to tell me what I want to know?

This question is far more difficult. Stan is a fast algorithm, which means it trades speed for an added risk of inaccuracy. Stan, by construction, will more often be correct than incorrect. There are other algorithms that are designed to search the parameter space widely for local maximums which may be more accurate, but which will be very slow.

What you should do, before using a particular algorithm, is read the literature on that algorithm and look at its functional limitations. Unfortunately, that is usually mathematical work as the only real goal of any non-conjugate method is to estimate $$\int_{\theta\in\Theta}f(X|\theta)\pi(\theta)\mathrm{d}\theta.$$

The second thing you can do is to validate it with an alternative algorithm. The numbers will never match, but if you deem them close enough, then you are fine.

Third, most of the prebuilt packages provide warnings that something may be amiss. If a warning comes up, use something else after investigating the source of the problem, so you do not recreate it in another algorithm.

Fourth, look at your prior density. Imagine you had a prior density of $\Pr(\mu)=\mathcal{N}(7,2^2)$ with $\sigma^2$ known just to simplify it and a likelihood of $\mathcal{N}(25,.1^2)$. At the least, you should be going wow, either I was wrong, the sample was bad, or there is something else going on that I should investigate.

Fifth, and you should do this before you start Stan in the first place, graph out your marginal likelihoods in one or two dimensions. Are there surprises anywhere that may interfere with the algorithm?

Since Bayesians don't define probability as what we see in the long run, how can I use simulation methods to verify than stan_glm is accurately capturing uncertainty? That is, how could I trust that these credible intervals are valid, using simulation methods? And right now, I'm not even defining a prior—how does the inclusion of priors come into play here, since that will affect our measures of uncertainty?

If you do not define a prior, then your model is not valid. If you are not defining a reasonable prior density, then why would you use a Bayesian model? Frequentist models minimize the risk of the maximum loss that could happen from gathering a bad sample. They are very pessimistic and it often takes more information to produce the same result a Bayesian method would.

Nonetheless, that is of no use without using a good prior density. The prior density allows the Bayesian method to minimize the average loss from choosing a bad sample. The information in the prior acts as a weighting scheme so that if some extreme sample is chosen by unfortunate chance, the prior weakens the role that the data plays.

EDIT I realized I didn't provide one specific answer. It was to the question

How could I use simulation methods to verify that what I'm doing in Stan is actually going to tell me what I want to know?

What makes this question challenging is that in the Bayesian paradigm the fixed points are ,$X$, the sample. In Frequentist methods, the parameters are fixed and thousands of unseen samples are created. On the Bayesian side of the coin, it is the sample which is fixed. You need to simulate thousands of parallel universes.

To see what that may be like, imagine all density functions of a coin toss with an unknown probability $p$ of being heads and $1-p$ of being tails. You observe six heads and two tails. Imagine a small parameter space where $p\in\{1/3,1/2,2/3\}$. Your simulation would consider all the cases where six heads could be obtained over the three objective binomial distributions. The posterior would be the weighted average of each parameter being the true value. Your predictive distribution would be the sum of the weighted binomial distributions.

Of importance to you, it is impossible for the Bayesian prediction to ever be the true distribution. One of the three distributions is the true distribution. The Bayesian methods weight their probability based on the observed value and the prior. The posterior can never be the true distribution, nor the predictive density.

It is asking "what is the probability of seeing six heads and two tails over the set of all possible explanations (parameters, models, etc)."

The Frequentist would assert one of the three choices was the true value by making it the null. Six heads and two tails would falsify $H_0:p=1/3,$ but not the others. If, by chance, you chose the correct one of the three distributions, then you are perfectly correct. Otherwise, you will be wrong.

If you would use simulations to hold a sample fixed, you would find that Stan would perform admirably as Bayes theorem is a mathematical theorem. It is ex-post optimal. All you would find is that the algorithm correctly implemented Bayes theorem up to the natural error level in estimating the denominator.

There are three things you can do. First, you can use model scoring methods for out-of-sample data. Second, you can use a Bayesian model selection or model averaging process. Third, you can treat it as a Frequentist problem and construct the sampling distribution of estimators.

For the first, scoring methods are an entire literature unto itself. You should research them. Bayesian model selection and model averaging treat models as parameters. For model selection, the probability of the models being true is calculated. For model averaging the probability each model is true is calculated and that serves as weighting over the model space. Finally, you can treat it as a Frequentist model.

The last one will be a problem in many standard cases because of the prior. For models with three or more dimensions and a normal distribution, the posterior density will not integrate to unity if the prior density isn't a proper density. In other words, you have to bite the bullet and choose a prior for any model with any real complexity.

The presence of a correctly centered proper prior forces the case where the Bayesian method will be superior to the corresponding Frequentist method due to the improved information. The Bayesian method will win under any reasonable standard. That isn't due to a flaw in the Frequentist method, but the Bayesian method assumes exterior information. The Frequentist method, by only considering the information in the sample, will have less information in it if you have a real prior.

Again, if you do not have a real prior, then why are you using a Bayesian method?