Bayesian Methods Part 1: An Introduction to Bayesian and Likelihood Thinking

Have you ever made a decision and thought, “there’s a 5% chance I’m wrong”? In Frequentist statistics, that 5% error rate is calculated before we see the data. But once we’ve observed our data and made our call, are we actually wrong 5% of the time? Spoiler: not necessarily. In fact, we might be either certainly right or certainly wrong, and the 5% doesn’t help us one bit.

In this tutorial, we’ll walk through two hypothesis testing problems that expose this limitation and then introduce likelihood functions — a tool that lets us ask a much more natural question: given what I observed, which hypothesis does the data actually support? This is our first step toward Bayesian thinking.

Why Should We Care?

Before diving into the math, let’s think about why this matters.

Suppose we’re researchers trying to understand whether a new advising intervention reduces attrition among first-generation college students. We run our study, collect the data, and a standard hypothesis test returns \(p = 0.06\). Under the conventional \(\alpha = 0.05\) threshold, we “fail to reject” the null — and the intervention looks like it didn’t work. But should we really conclude that? Our sample was modest, prior literature suggests the intervention should help, and the effect size looks meaningful. The Frequentist framework gives us a binary verdict based on long-run error rates, but it can’t incorporate what we already knew going in, and it can’t tell us what to believe about this particular result.

This tension runs through much of education research. We study complex, multilevel systems — students nested in classrooms, classrooms nested in institutions — where sample sizes for specific subgroups (e.g., undocumented students, formerly incarcerated learners) can be small, prior knowledge from existing studies is abundant, and the stakes of getting it wrong are high. A Frequentist \(p\)-value tells us how often we’d see data this extreme if the null were true and we repeated the experiment forever. It says nothing about the probability that the intervention actually works for the students in front of us.

As Hoff (2009, Ch. 1) puts it, Frequentist inference focuses on what would happen across many datasets, not what our particular dataset tells us. Bayesian inference, by contrast, asks the more direct question: given this data and what we already know, what should we believe about the effect? For education researchers concerned with equity and student success, this shift is powerful — it means we can formally bring prior evidence into our analyses, make honest probability statements about effects for small subpopulations, and move beyond the blunt instrument of “significant or not.”

The examples below will make this distinction concrete.

1. Example 1: Earthquake Survival Kits (Binomial)

1.1 The Problem

Imagine we’re running quality control for earthquake survival kits in Los Angeles County. Our supplier is either providing kits where one-third are defective or where two-thirds are defective. Let \(\theta\) be the true defect rate. We want to test:

\[H_0 : \theta = \frac{1}{3} \qquad H_1 : \theta = \frac{2}{3}\]

We randomly sample 100 kits and count how many fail. Let \(X\) be the number of failures.

1.2 Step 1: Identify the Distribution of \(X\)

The first thing we need to figure out: what kind of random variable is \(X\)?

Think about it this way — each kit either fails or doesn’t (two outcomes). We’re testing 100 kits independently, each with the same defect probability \(\theta\). This is exactly the setup for a binomial distribution:

• Under \(H_0\) (\(\theta = 1/3\)): \(X \sim \text{Binomial}(100,\; 1/3)\)
• Under \(H_1\) (\(\theta = 2/3\)): \(X \sim \text{Binomial}(100,\; 2/3)\)

TipQuick Refresher: Binomial Distribution

If we flip a coin \(n\) times with probability \(p\) of heads, the total number of heads follows \(\text{Binomial}(n, p)\). The probability of observing exactly \(x\) heads is:

\[P(X = x) = \binom{n}{x} p^x (1-p)^{n-x}\]

In our problem, each kit is a “coin flip” with probability \(\theta\) of being defective, and we flip 100 times. See Hoff (2009, 2.4.1) for a detailed treatment of the binomial distribution.

Under \(H_0\), we’d expect about \(100 \times 1/3 \approx 33\) failures. Under \(H_1\), we’d expect about \(100 \times 2/3 \approx 67\) failures. So the data should look very different under the two hypotheses — the question is how we formalize the decision.

1.3 Step 2: Understand the Decision Rule

Our decision rule is:

• Accept \(H_0\) if \(x < 50\) (fewer than 50 failures → looks more like \(\theta = 1/3\))
• Reject \(H_0\) if \(x > 50\) (more than 50 failures → looks more like \(\theta = 2/3\))
• Randomly choose \(H_0\) or \(H_1\) (50/50 coin flip) if \(x = 50\)

The boundary at 50 makes intuitive sense — it’s right between what we’d expect under \(H_0\) (about 33 failures) and \(H_1\) (about 67 failures). The coin-flip at \(x = 50\) is called a randomized decision rule, a theoretical device to handle the exact boundary case.

1.4 Step 3: Calculate Type I and Type II Errors

Now let’s compute the error probabilities. Let’s make sure we’re clear on what these mean:

• Type I error (\(\alpha\)): Rejecting \(H_0\) when \(H_0\) is actually true. This means concluding “two-thirds are defective” when really only one-third are. The consequence: we might unnecessarily switch suppliers, wasting money.
• Type II error (\(\beta\)): Accepting \(H_0\) when \(H_1\) is actually true. This means concluding “only one-third are defective” when actually two-thirds are. The consequence: defective earthquake kits go out to the public — potentially dangerous.

Computing \(\alpha\) (Type I error):

We reject \(H_0\) when \(X > 50\) or (with probability 0.5) when \(X = 50\). So under \(H_0\) (\(\theta = 1/3\)):

\[\alpha = P(X > 50 \mid \theta = 1/3) + 0.5 \cdot P(X = 50 \mid \theta = 1/3)\]

Computing \(\beta\) (Type II error):

We accept \(H_0\) when \(X < 50\) or (with probability 0.5) when \(X = 50\). So under \(H_1\) (\(\theta = 2/3\)):

\[\beta = P(X < 50 \mid \theta = 2/3) + 0.5 \cdot P(X = 50 \mid \theta = 2/3)\]

Let’s compute these in R:

# --- Type I error ---
# P(X > 50 | θ = 1/3): pbinom gives P(X ≤ 50), so subtract from 1
p_gt50_h0 <- 1 - pbinom(50, 100, 1/3)

# P(X = 50 | θ = 1/3): the exact point probability
p_eq50_h0 <- dbinom(50, 100, 1/3)

# Total: always reject when X > 50, reject half the time when X = 50
alpha <- p_gt50_h0 + 0.5 * p_eq50_h0

# --- Type II error ---
# P(X < 50 | θ = 2/3) = P(X ≤ 49 | θ = 2/3)
p_lt50_h1 <- pbinom(49, 100, 2/3)

# P(X = 50 | θ = 2/3)
p_eq50_h1 <- dbinom(50, 100, 2/3)

# Total: always accept H0 when X < 50, accept half the time when X = 50
beta <- p_lt50_h1 + 0.5 * p_eq50_h1

cat("Type I error (alpha):", alpha, "\n")

Type I error (alpha): 0.0003091369 

cat("Type II error (beta):", beta, "\n")

Type II error (beta): 0.0003091369 

TipR Functions: pbinom vs dbinom

These are our two key tools for binomial calculations:

• dbinom(x, n, p) = \(P(X = x)\) — the point probability (density). The “d” stands for density.
• pbinom(q, n, p) = \(P(X \leq q)\) — the cumulative probability (CDF). The “p” stands for probability.

So 1 - pbinom(50, 100, 1/3) gives \(P(X > 50)\) by subtracting \(P(X \leq 50)\) from 1. This pattern — 1 - pbinom(...) for upper-tail probabilities — is one we’ll use constantly.

Both error rates come out to approximately 0.00031. They’re equal. This is not a coincidence.

Why are \(\alpha\) and \(\beta\) the same? This happens because of a beautiful symmetry. If \(Y \sim \text{Bin}(100, 1/3)\), then \(W = 100 - Y \sim \text{Bin}(100, 2/3)\). So:

\[P(Y > 50 \mid \theta = 1/3) = P(100 - Y < 50 \mid \theta = 1/3) = P(W < 50 \mid \theta = 2/3)\]

The two hypotheses (\(\theta = 1/3\) and \(\theta = 2/3\)) are mirror images around the decision boundary (\(50/100 = 1/2\)). When both errors are equal, this single value serves as the overall Frequentist error probability for the test — no trade-off needed.

1.5 Step 4: The Conditional Error at \(x = 50\)

Here’s where things get interesting. If we observe exactly \(x = 50\), the rule flips a coin. What’s the actual probability we make an error?

• If \(H_0\) is true (\(\theta = 1/3\)): we reject with probability 0.5 → \(P(\text{error} \mid x = 50, H_0) = 0.5\)
• If \(H_1\) is true (\(\theta = 2/3\)): we accept \(H_0\) with probability 0.5 → \(P(\text{error} \mid x = 50, H_1) = 0.5\)

Either way, the conditional probability of error given \(x = 50\) is 0.5 — vastly different from the tiny unconditional error rate of ~0.031%.

NoteWhy Does This Matter?

The unconditional error rate of ~0.03% averages over all possible values of \(x\). For most values of \(x\) (like \(x = 30\) or \(x = 70\)), the decision is obviously correct. But when \(x = 50\), the data is genuinely ambiguous — and the 0.03% figure hides that ambiguity. This is the gap between pre-data guarantees and post-data reality.

2. Example 2: Uniform Distribution

Let’s try a different setup to see the same phenomenon even more clearly.

2.1 The Problem

Now \(X \sim U(\theta, \theta + 1)\) — uniformly distributed between \(\theta\) and \(\theta + 1\) — and we test:

\[H_0 : \theta = 0 \qquad H_1 : \theta = 0.9\]

Step 1: What does \(X\) look like under each hypothesis?

• Under \(H_0\) (\(\theta = 0\)): \(X \sim U(0, 1)\) — \(X\) is equally likely to be anywhere between 0 and 1
• Under \(H_1\) (\(\theta = 0.9\)): \(X \sim U(0.9, 1.9)\) — \(X\) is equally likely to be anywhere between 0.9 and 1.9

Notice the key feature: these two distributions overlap on the interval \((0.9, 1)\). If \(X\) falls in this region, the data are consistent with both hypotheses.

Step 2: The decision rule.

Accept \(H_0\) if \(x < 0.95\), reject \(H_0\) if \(x \geq 0.95\). The boundary 0.95 sits right in the middle of the overlap region.

2.2 Step 3: Type I and Type II Errors

For the uniform distribution, computing probabilities is about interval lengths. The CDF of \(X \sim U(a, b)\) is:

\[P(X \leq x) = \frac{x - a}{b - a}\]

This is just the fraction of the total interval that we’ve covered — a direct consequence of “uniform” meaning “equal probability everywhere.”

Type I error (reject \(H_0\) when \(H_0\) is true, i.e., \(X \sim U(0, 1)\)):

\[\alpha = P(X \geq 0.95 \mid \theta = 0) = \frac{1 - 0.95}{1 - 0} = \frac{0.05}{1} = 0.05\]

Type II error (accept \(H_0\) when \(H_1\) is true, i.e., \(X \sim U(0.9, 1.9)\)):

\[\beta = P(X < 0.95 \mid \theta = 0.9) = \frac{0.95 - 0.9}{1.9 - 0.9} = \frac{0.05}{1} = 0.05\]

Let’s verify in R:

# Type I error: P(X >= 0.95 | X ~ U(0,1))
# punif gives P(X <= x), so subtract from 1
alpha_q2 <- 1 - punif(0.95, 0, 1)
alpha_q2

[1] 0.05

# Type II error: P(X < 0.95 | X ~ U(0.9, 1.9))
beta_q2 <- punif(0.95, 0.9, 1.9)
beta_q2

[1] 0.05

TipR Function: punif

Just like pbinom is the CDF for the binomial, punif(x, min, max) is the CDF for the uniform distribution: it returns \(P(X \leq x)\) where \(X \sim U(\text{min}, \text{max})\).

Both errors equal 0.05. A nice, clean Frequentist result. But does that 5% actually describe how accurate our specific decision is?

2.3 Step 4: The Post-Data Reality

Now the crucial question: if we observe \(0.9 < x < 1\) — a value in the overlap region, consistent with both hypotheses — what is the conditional probability that our decision is wrong?

We need to condition on the event \(\{0.9 < X < 1\}\) and compute the error probability within this region. The decision boundary at \(x = 0.95\) splits the overlap region \((0.9, 1)\) into two halves of equal length (0.05 each), so:

Conditional Type I error (reject \(H_0\) when \(H_0\) is true, given \(0.9 < x < 1\)):

\[P(X > 0.95 \mid X \sim U(0,1), \; 0.9 < X < 1) = \frac{1 - 0.95}{1 - 0.9} = \frac{0.05}{0.1} = \frac{1}{2}\]

Conditional Type II error (accept \(H_0\) when \(H_1\) is true, given \(0.9 < x < 1\)):

\[P(X < 0.95 \mid X \sim U(0.9, 1.9), \; 0.9 < X < 1) = \frac{0.95 - 0.9}{1 - 0.9} = \frac{0.05}{0.1} = \frac{1}{2}\]

Either way, the conditional probability of error given \(0.9 < x < 1\) is 1/2 — vastly different from the unconditional 5%.

Notice the parallel with Example 1: there, the conditional error at \(x = 50\) was also 1/2, compared to the unconditional ~0.031%. Both examples tell the same story — the pre-data error rate can be dramatically different from the post-data error probability once we condition on what we actually observed.

The point is: given we observe \(0.9 < x < 1\), do we continue to regard the probability of making an error to be 0.05, or do we regard it to be \(\frac{1}{2}\)? The latter typifies the conditional perspective that motivates Bayesian inference.

ImportantThe Frequentist Limitation

The pre-data error rate of 0.05 describes long-run performance across many experiments. It does not describe the accuracy of any particular decision after the data are observed. This disconnect is one of the key motivations for Bayesian inference, where we can compute \(P(\text{hypothesis} \mid \text{data})\) directly.

As Hoff (2009, §1.1) describes, Bayesian learning starts with a prior distribution \(p(\theta)\) representing our beliefs before seeing data, then uses Bayes’ rule to update:

\[p(\theta \mid y) = \frac{p(y \mid \theta) \, p(\theta)}{p(y)}\]

This gives us something Frequentist inference cannot: a direct probability statement about \(\theta\) given our observed data \(y\). Bayes’ rule “does not tell us what our beliefs should be, it tells us how they should change after seeing new information” (Hoff, 2009, §1.1).

3. Enter the Likelihood Function

So if Frequentist error rates don’t tell us about a specific observation, what does? This is where likelihood functions come in — the bridge between Frequentist and Bayesian statistics.

3.1 What Is a Likelihood Function and Why Should We Care?

In our intro stats course, we learned to think about \(p(x \mid \theta)\) — the probability of data \(x\) given a parameter \(\theta\). We fix \(\theta\) and ask “how likely is this data?”

The likelihood function flips this perspective. We fix the data \(x\) (because we’ve already observed it.) and ask: “how well does each possible value of \(\theta\) explain what we saw?”

\[L(\theta \mid x) = p(x \mid \theta), \quad \text{viewed as a function of } \theta\]

Same formula, completely different interpretation. This shift — from “what data would I expect?” to “what parameter does my data support?” — is fundamental.

NoteFormal Definition (Hoff, 2009)

Following Hoff’s framework (Ch. 1-3), the components of Bayesian inference are:

1. Prior distribution \(p(\theta)\): our beliefs about \(\theta\) before seeing data
2. Sampling model \(p(y \mid \theta)\): how data would be generated given \(\theta\) — this becomes the likelihood once data are observed
3. Posterior distribution \(p(\theta \mid y)\): our updated beliefs after seeing data

The likelihood \(p(y \mid \theta)\) is the ingredient that connects data to parameters. It’s the piece that both Frequentists and Bayesians agree on — the disagreement is about what else we bring to the table (nothing vs. a prior).

Why should we care? Because the likelihood function lets the data speak for themselves. Instead of asking “would my procedure be right in the long run?” (Frequentist), we can ask “which parameter value best explains this data?” — and that’s the question we usually actually want to answer.

3.2 Likelihood for the Binomial Problem (Example 1)

For the binomial model \(X \sim \text{Bin}(100, \theta)\), the likelihood given observation \(X = x\) is:

\[L(\theta \mid x) = \binom{100}{x} \theta^x (1-\theta)^{100-x}\]

Since we only have two possible values of \(\theta\), we just evaluate this at each:

\[L\!\left(\frac{1}{3} \;\middle|\; x\right) = \binom{100}{x} \left(\frac{1}{3}\right)^x \left(\frac{2}{3}\right)^{100-x}\]

\[L\!\left(\frac{2}{3} \;\middle|\; x\right) = \binom{100}{x} \left(\frac{2}{3}\right)^x \left(\frac{1}{3}\right)^{100-x}\]

Let’s see how these compare for a few values of \(x\):

# Likelihood under each hypothesis for different x values
x_vals <- c(30, 40, 50, 60, 70)

for (x in x_vals) {
  L_h0 <- dbinom(x, 100, 1/3)   # L(θ=1/3 | x)
  L_h1 <- dbinom(x, 100, 2/3)   # L(θ=2/3 | x)
  ratio <- L_h0 / L_h1           # likelihood ratio
  cat(sprintf("x = %d: L(1/3|x) = %.2e, L(2/3|x) = %.2e, ratio = %.2e\n",
              x, L_h0, L_h1, ratio))
}

x = 30: L(1/3|x) = 6.73e-02, L(2/3|x) = 6.12e-14, ratio = 1.10e+12
x = 40: L(1/3|x) = 3.08e-02, L(2/3|x) = 2.93e-08, ratio = 1.05e+06
x = 50: L(1/3|x) = 2.20e-04, L(2/3|x) = 2.20e-04, ratio = 1.00e+00
x = 60: L(1/3|x) = 2.93e-08, L(2/3|x) = 3.08e-02, ratio = 9.54e-07
x = 70: L(1/3|x) = 6.12e-14, L(2/3|x) = 6.73e-02, ratio = 9.09e-13

TipLikelihood Ratio Intuition

Comparing \(L(1/3 \mid x)\) to \(L(2/3 \mid x)\) tells us which hypothesis the data support more strongly:

• \(x = 30\): Ratio is huge → data overwhelmingly favor \(\theta = 1/3\)
• \(x = 50\): Ratio is closer to 1 → data are ambiguous
• \(x = 70\): Ratio is tiny → data overwhelmingly favor \(\theta = 2/3\)

The likelihood ratio is the core quantity in Bayesian hypothesis comparison. As we’ll see in later chapters of Hoff (2009), the Bayes factor — the ratio of likelihoods integrated over priors — is the Bayesian analog of the p-value for comparing hypotheses.

3.3 Likelihood for the Uniform Problem (Example 2)

For the uniform model \(X \sim U(\theta, \theta + 1)\), the density is either 1 (if \(x\) is in the support) or 0 (if not). So the likelihood is simple:

\[L(\theta \mid x) = f(x \mid \theta) = \begin{cases} 1 & \text{if } \theta \leq x \leq \theta + 1 \\ 0 & \text{otherwise} \end{cases}\]

Evaluating at our two hypotheses:

\[L(0 \mid x) = \begin{cases} 1 & \text{if } 0 \leq x \leq 1 \\ 0 & \text{otherwise} \end{cases} \qquad L(0.9 \mid x) = \begin{cases} 1 & \text{if } 0.9 \leq x \leq 1.9 \\ 0 & \text{otherwise} \end{cases}\]

3.4 Step-by-Step Interpretation: What the Likelihood Tells Us

This is the key payoff. The likelihood function partitions the sample space into three regions. Let’s walk through each one carefully:

Region 1: \(0 \leq x < 0.9\)

• \(L(0 \mid x) = 1\) ✓ (x falls within [0, 1])
• \(L(0.9 \mid x) = 0\) ✗ (x does NOT fall within [0.9, 1.9])
• Interpretation: Only \(H_0\) has nonzero likelihood. The data are impossible under \(H_1\), so we know with certainty that \(\theta = 0\). No ambiguity whatsoever.

Region 2: \(0.9 \leq x \leq 1\) (the overlap region)

• \(L(0 \mid x) = 1\) ✓
• \(L(0.9 \mid x) = 1\) ✓
• Interpretation: Both hypotheses have equal likelihood. The data are perfectly consistent with either \(\theta = 0\) or \(\theta = 0.9\). The data alone provide no information to distinguish the hypotheses. Any decision we make in this region must rely on something other than the data — this is exactly where a prior distribution would help!

Region 3: \(1 < x \leq 1.9\)

• \(L(0 \mid x) = 0\) ✗ (x does NOT fall within [0, 1])
• \(L(0.9 \mid x) = 1\) ✓
• Interpretation: Only \(H_1\) has nonzero likelihood. The data are impossible under \(H_0\), so we know \(\theta = 0.9\) with certainty.

Here’s a summary table:

Region	\(L(0 \mid x)\)	\(L(0.9 \mid x)\)	What the data say
\(0 \leq x < 0.9\)	1	0	Definitely \(H_0\)
\(0.9 \leq x \leq 1\)	1	1	Can’t tell — need more info
\(1 < x \leq 1.9\)	0	1	Definitely \(H_1\)

NoteThe Honest Answer

Notice something remarkable: in the overlap region, the likelihood honestly tells us “the data cannot distinguish these hypotheses.” No decision rule can change that fact — it’s a property of the data itself. The Frequentist 5% error rate masks this uncertainty by averaging over all possible data, including the unambiguous cases.

This is where Bayesian inference shines. If we had a prior — say, \(P(H_0) = 0.5\) and \(P(H_1) = 0.5\) — then even in the overlap region, we could compute a meaningful posterior probability for each hypothesis via Bayes’ rule. The prior encodes our beliefs before seeing the data, and the posterior tells us what we should believe after. We’ll explore this in future posts.

4. The Big Picture: From Frequentist to Bayesian

Let’s step back and see the full arc of what we’ve learned:

	Frequentist View	Likelihood / Bayesian View
What’s fixed?	\(\theta\) is fixed, data varies	Data is fixed, \(\theta\) varies
Key question	“How often will I be wrong?”	“Given this data, what should I believe?”
Error assessment	Pre-data, long-run average	Post-data, specific to what was observed
When data are ambiguous	Still applies the 5% rate	Honestly says “data can’t tell”
Needs a prior?	No	Yes (for full Bayesian inference)

The likelihood function is the bridge between these two worlds. As Hoff (2009, §1.1) describes the idealized form of Bayesian learning:

1. Start with a prior distribution \(p(\theta)\) — our beliefs before seeing data
2. Specify a sampling model \(p(y \mid \theta)\) — this becomes our likelihood once data are observed
3. Compute the posterior distribution \(p(\theta \mid y)\) via Bayes’ rule — our updated beliefs

\[p(\theta \mid y) = \frac{p(y \mid \theta) \, p(\theta)}{\int p(y \mid \tilde{\theta}) \, p(\tilde{\theta}) \, d\tilde{\theta}}\]

Hoff emphasizes a pragmatic view: even when our prior \(p(\theta)\) only approximates our true beliefs, the posterior is still useful. And in many complex problems, “there are no obvious non-Bayesian methods of estimation or inference” — Bayes’ rule provides a principled way to generate procedures, which can then be evaluated even by Frequentist criteria (Hoff, 2009, §1.2).

TipKey Takeaways

1. Frequentist error rates (\(\alpha\), \(\beta\)) describe long-run performance, not the accuracy of any specific decision.
2. Conditional on the observed data, we may be certainly right or certainly wrong — the pre-data error rate doesn’t apply.
3. The likelihood function lets us see what the data actually say about each hypothesis value, without reference to long-run error rates.
4. Bayesian inference combines the likelihood with a prior to give direct probability statements about \(\theta\) given the data — the question we usually want to answer.

In the next post, we’ll explore how prior distributions work and see the full Bayesian machinery in action using conjugate models from Hoff Ch. 3.

5. Bayesian Methods in Education Research: Where the Field Is Heading

The ideas we’ve developed in this tutorial — likelihood functions, the limitations of Frequentist error rates, and the promise of updating beliefs with data — are not just theoretical curiosities. They are increasingly shaping how researchers study real educational problems: student retention, achievement gaps, graduation rates, and school effectiveness.

Paper 1: Crisp, G., Doran, E., & Salis Reyes, N. A. (2018). Predicting graduation rates at 4-year broad access institutions using a Bayesian modeling approach. Research in Higher Education, 59(2), 133–155.

• Research question: What institutional characteristics predict six-year graduation rates at four-year broad access institutions (BAIs), and do the same predictors generalize across overall, African American, and Latina/o graduation rates?
• Data: Two panels of IPEDS data — Panel 1: 412 BAIs using 2007–2008 data; Panel 2: 402 BAIs using 2014–2015 data — with 24 candidate institutional-level predictors (e.g., enrollment size, % full-time students, SES composition, institutional expenditures, religious affiliation).
• Bayesian methods:
  • Bayesian model averaging (BMA) implemented via the BMS package in R, using Markov Chain Monte Carlo (MCMC) model composition to explore the space of \(2^{24}\) candidate regression models.
  • Each candidate model is weighted by its posterior model probability (PMP) given the data; predictions are averaged across all models rather than selecting a single “best” specification.
  • Posterior inclusion probabilities (PIPs) are computed for each predictor and evaluated using Raftery’s (1995) criteria: PIP > 0.50 (weak evidence), > 0.75 (positive), > 0.95 (strong), > 0.99 (very strong).
  • A uniform model prior (each model equally likely a priori) and the unit information prior (UIP) for coefficients are used as defaults.
• Key findings: For overall graduation rates, key predictors with high PIPs included religious affiliation, % full-time students, SES composition, enrollment size, and institutional revenue/expenditures. Fewer variables reliably predicted race-specific graduation rates — e.g., only 2–3 predictors exceeded PIP > 0.50 for African American and Latina/o rates, compared to 5–7 for overall rates. The BMA approach produced more robust predictions than conventional stepwise model selection, which is especially important for BAIs where subpopulation samples are small and overfitting risk is high.

Paper 2: Matheny, K. T., Thompson, M. E., Townley-Flores, C., & Reardon, S. F. (2023). Uneven progress: Recent trends in academic performance among U.S. school districts. American Educational Research Journal, 60(3), 447–485.

• Research question: How have district-level academic achievement and achievement disparities (by race/ethnicity and socioeconomic status) changed across U.S. school districts from 2009 to 2019, and are improvements in achievement and reductions in disparities synergistic or in tension?
• Data: Stanford Education Data Archive (SEDA), based on approximately 430 million Grades 3–8 math and reading test scores from over 90,000 U.S. public schools aggregated into ~12,800 geographic school districts, supplemented by district-level covariates from the Common Core of Data (CCD), the Civil Rights Data Collection (CRDC), and the American Community Survey (ACS).
• Bayesian methods:
  • Precision-weighted random effects (hierarchical linear) models to estimate district-level achievement trends, where SEDA’s empirical Bayes (EB) shrinkage stabilizes noisy estimates for smaller districts by pulling them toward the grand mean.
  • Within-district model: \(\hat{\mu}_{dgyb} = \beta_{0d} + \beta_{1d}(grade) + \beta_{2d}(cohort) + \beta_{3d}(math) + e_{dgyb} + \varepsilon_{dgyb}\), where \(\varepsilon_{dgyb}\) is the SEDA-provided estimation error treated as known, and \(\beta_{2d}\) (the cohort trend) is the parameter of interest.
  • Between-district model relates estimated trends to district covariates via precision-weighted regressions, downweighting low-precision estimates.
  • Fitted using HLM v7 software via maximum likelihood, with the EB framework ensuring that districts with limited data borrow strength from the national distribution.
• Key findings: Average test scores improved very modestly (~0.001 SD/year) in the typical district, but trends varied considerably across districts. Socioeconomic (nonpoor–poor) disparities widened by about 0.05 SD over the decade; White-Black disparities grew by ~0.03 SD; White-Hispanic disparities narrowed by ~0.05 SD. Trends in overall achievement and trends in disparities were essentially uncorrelated (correlations near zero), indicating no systematic synergy or trade-off between performance gains and equity gains. The strongest predictors of increasing inequality were within-district racial and socioeconomic segregation and differential access to certified teachers.

Paper 3: Bowman, N. A., Logel, C., Lacosse, J., Canning, E. A., Emerson, K. T. U., & Murphy, M. C. (2023). The role of minoritized student representation in promoting achievement and equity within college STEM courses. AERA Open, 9, 1–16.

• Research question: To what extent does the course-level representation of underrepresented racial minority (URM) and first-generation college students predict individual grades in postsecondary STEM courses, and do these relationships vary by students’ own minoritized or privileged identities?
• Data: 87,027 course grades received by 11,868 STEM-interested undergraduates in 8,468 STEM courses at 20 colleges and universities, drawn from the College Transition Collaborative’s social-belonging dataset (fall 2015 and fall 2016 cohorts). Student-level measures include race/ethnicity, first-generation status, sex, and ACT/SAT scores; course-level measures include proportion of URM and first-generation students, class size, STEM discipline, and academic term.
• Bayesian methods:
  • Cross-classified multilevel models (using the lmer package in R) to handle the non-hierarchical nesting structure where students and courses are crossed (each student takes multiple courses, and each course enrolls multiple students), with institutions at level 3.
  • While estimated via maximum likelihood rather than full Bayesian MCMC, the cross-classified multilevel framework embodies the Bayesian logic of partial pooling: estimates for individual courses and students are shrunk toward group means, borrowing strength across units — the same principle underlying empirical Bayes estimation.
  • Student fixed effects analyses (via the xtreg package in Stata) as robustness checks, removing all between-student variation to address self-selection concerns.
  • Benjamini-Hochberg procedure for false discovery rate correction across multiple comparisons.
• Key findings: Higher proportions of URM and first-generation students in a STEM course were positively associated with grades among all students, and these positive relationships were significantly stronger for students who themselves held minoritized identities. Specifically, the URM grade gap was 27% smaller in high- versus low-URM-representation courses, and the first-generation gap was reduced by over half (56%). Results were robust to student fixed effects and multiple comparison corrections, and were generally consistent across STEM disciplines, class sizes, and student subgroups. These findings support the importance of classroom-level representation as a contextual factor shaping equity in STEM outcomes.

TipKey Takeaway for Education Researchers

Bayesian and empirical Bayes methods are not just a different way to get \(p\)-values. They offer a fundamentally different relationship with uncertainty — one where prior knowledge and population-level information are features rather than nuisances, where small samples are handled gracefully through partial pooling, and where the output can be a direct probability statement about the quantity we care about. For researchers studying equity, access, and student success in higher education, these are not abstract advantages; they are practical tools for producing more honest and useful evidence.

References

• Bowman, N. A., Logel, C., Lacosse, J., Canning, E. A., Emerson, K. T. U., & Murphy, M. C. (2023). The role of minoritized student representation in promoting achievement and equity within college STEM courses. AERA Open, 9, 1–16.
• Crisp, G., Doran, E., & Salis Reyes, N. A. (2018). Predicting graduation rates at 4-year broad access institutions using a Bayesian modeling approach. Research in Higher Education, 59(2), 133–155.
• Hoff, P. (2009). A First Course in Bayesian Statistical Methods. Springer. (Chapters 1-3 for the foundations discussed here.)
• Matheny, K. T., Thompson, M. E., Townley-Flores, C., & Reardon, S. F. (2023). Uneven progress: Recent trends in academic performance among U.S. school districts. American Educational Research Journal, 60(3), 447–485.

This tutorial is based on UCLA STATS C216: Applied Bayesian Social Statistics taught by Prof Mark Handcock.