I’m looking at mean differences between salary two groups (Male vs. Female) in a synthetic dataset. Gender was determined randomly, so the ground truth is that there is no gender difference.
At the end of the day, I’m looking for bayes factors assessing whether the evidence better supports the hypothesis that the difference between groups is centered on zero vs. centered on some alternative critical value (e.g., 10,000). The approach I was thinking:
- Create two priors: Null, which is normally distributed and centered on 0, and Alternative, which is normally distributed and centered on 10,000.
- Sample the posterior distributions and compute a Bayes factor to quantify how much more the data favor the starting assumption that the difference is centered over 0 vs 10,000.
I have already used an approach that samples the male and female data from a gamma distribution to generate a distribution of differences, which is normal and roughly zero:
with pm.Model() as diff_model:
M_data=df.loc[df['gender']=='m']['wealth'].to_numpy()
F_data=df.loc[df['gender']=='f']['wealth'].to_numpy()
#priors for means
mean_M=pm.Gamma("mean_M", mu=mu, sigma=sd)
mean_F=pm.Gamma("mean_F", mu=mu, sigma=sd)
#priors for variances
sigma_M=pm.HalfNormal("sigma_M", sd)
sigma_F=pm.HalfNormal("sigma_F", sd)
# Likelihoods
wealth_M= pm.Gamma('wealth_M', mu=mean_M, sigma=sigma_M, observed = M_data)
wealth_F= pm.Gamma('wealth_F', mu=mean_F, sigma=sigma_F, observed = F_data)
# Compute the difference and store it in the trace
mean_diff = pm.Deterministic('mean_diff', mean_M - mean_F)
sigma_diff = pm.Deterministic('sigma_diff', sigma_M - sigma_F)
#idata = pm.sample() #takes forever
n_tune=2000
n_chains=4
n_samples=2000
idata = pm.sample(init='ADVI', chains=n_chains, draws=n_samples, tune=n_tune)
This is fine for showing that the group differences overlap with zero, but really what I want to know here is how to get a Bayes Factor out of my sampling experiment so I can quantify the evidence for/against a hypothesis rather than eyeball the distribution plots. Back when I was exploring Bayesian statistics for the first time in R, I was using a package that included a Bayesian version of the t-test and ANOVA, and produced Bayes factors that were a good gateway statistic for students used to thinking about p-values.
So; here’s the thing. Bayes factors are kinda neat,but theyre also pretty unreliable for model selection. In general you should keep in mind that they dislike complex models and are super sensitive to the prior (remember, you have to compute the marginal likelihood over the prior). Issues with misspecification none withstanding.
What do you mean about p values? Bayes factors and p values are not terribly related; there are some choices of priors that produce equivalent confidence intervals in glms, but they do not have a nice relationship in general. What sort of hypothesis test are you after?
You have the posteriors for the two groups, so why would you not compare them directly?
I want to thank you @brontidon for both THAT you replied, and HOW you replied; I understand from some of my reading that using Bayes factors in this way is not discussed in polite company for some of the reasons you describe (e.g., sensitivity to the quality of the priors).
The reason I bring up p-values at all is that I am advocating for a Bayesian approach and trying to meet the audience where they’re at. The community I am addressing are not hardcore academic theorists or statistical enthusiasts. These are folks who are currently running t-tests because that’s more rigorous than just eyeballing a bar chart, and at the end of the day, they only care whether the p-value matches what they’re looking for. This is purely a data validation exercise, and the test outcomes aren’t being used to develop new hypotheses or drive future exploration. In the area of fairness testing, we’re actually hoping for null effects (i.e., the outcomes are equivalent across different groups). This is what drew me to explore Bayesian approaches in the first place, since we’ve all been trained that traditional NHST techniques don’t prove the null hypothesis; they just fail to reject. This speaks to why I’m particularly interested in Bayes factors as quantifying evidence for the null hypothesis.
What I want to be able to do is use a BF to indicate the strength of evidence that the outcomes for k groups differ by zero vs. some concerning threshold(s). It would be very beneficial to be able to explain to a stakeholder that the odds are 100:1 that one of the groups will be favored by no more than X. Having a concrete value to point at might also prove useful for monitoring as new data becomes available over time. So far I have only figured out how to generate posterior distributions of differences (as in my original code snippet), but pulling an actual Bayes Factor out of that has remained elusive.
Another point that I will add is that these tests are often done with really large samples – like thousands or even tens of thousands of data points. As a result, the folks doing the testing are often obtaining significant test statistics that require some additional work to explain why p=.00000001 doesn’t necessarily mean that there’s a problematic difference.
Okay, this makes more sense. But I might point you into a direction that I think will cause less of a potential headache on your audience’s part: It’s nice because you’re also going to provide required materials under FDA/NSF revisions for Bayesian modeling that started a couple of years ago.
-
Drop everything bayes factors related UNLESS you want to explicitly show scenarios where they break down and introduce more rigorous information theoretic tools we use to determine if our model is go/no go. Because your audience is not as steeped into the stats side of things; I do not feel they will grasp some of the nuances that may lead them to conclude the opposite of what you intend.
-
Really focus on posterior predictive checks. Because if you ever publish anything bayesian; you need to provide measures of calibration; and your audience can understand this more intuitively. Remember: People Love Pictures (PLP) Perhaps mention loo-cv, but I don’t think unless there is a strong appetite for math that you’re gonna get super far on explaining that. WAIC, which to me has a very funny informal name is also easy enough to grasp the first pass around.
In my mind, because you are trying to illustrate the overall viability of your model against the data generating process; these are what you are after. Any statistical model is basically correct in its little plane of existence; you wanna see if it’s useful in this one. This can be a nice mention for causal inference for your students and why this is a very important experimental design consideration. Bayesian CI is both very cool and very deep though; so maybe small bites for now.
-
You should perhaps mention model averaging in the case you have a set of closely performing models. Show that! You just use your already computed posteriors to build a new super posterior (the posterior heights are the weights of your support values 
). I , of course assume that you have a motivated a decision function and have determined what makes a model good enough lol
-
The best take down of p values is just…using p values haha. Under the null, p values are UNIFORM. This is the starting point for a lot of practical simulation you could run with a couple lines of code (just fix a model where the null is true and badda bing badda boom, your audience will see that you’re gonna get type one rejections no matter what ). The p value under the alternative is much more interesting and difficult to prove in any given scenario.
Also, just a note on your model; you could simplify sampling a bit by modeling the differences explicitly, or using a hierarchical set up. I saw your note on sampling taking a slow time, and that’s mostly because you are basically reusing work.