Week 5: Intro to Bayesian Inference

PS 818 - Data Analysis with Statistical Models

Anton Strezhnev

University of Wisconsin-Madison

October 5, 2026

Intro to Bayesian Inference

\[ \DeclareMathOperator*{\argmin}{arg\,min} \DeclareMathOperator*{\argmax}{arg\,max} \]

Previously

  • Wrapping up MLE
    • Latent variables and discrete choice: ordinal logit, multinomial logit, and the conditional logit
    • Duration models: modelling the hazard, censoring, and the Cox proportional hazards model
    • Estimated via Maximum Likelihood
      • \(\theta\) is a fixed, unknown constant
      • Uncertainty comes from repeated sampling - the sampling distribution of \(\hat{\theta}\)

This week

  • Bayesian inference
    • What if we treat \(\theta\) as having uncertainty
    • We want the posterior distribution conditioning on our observed data \(\theta | \mathbf{Y}\)
    • But Bayes’ rule tells us that the likelihood is only half of this.
      • We also need a prior on \(\theta\).
  • Defining a Bayesian Model
    • Writing down the data-generating process
    • “Plate notation” for concise model summaries

Intro to Bayesian Inference

Bayesian inference

  • In the frequentist paradigm, we took a particular approach to thinking about randomness when trying to learn about an unknown parameter \(\theta\) using observed data \(\mathbf{Y}\)

    • \(\theta\) is an unknown constant
    • \(\mathbf{Y}\) are random variables
    • \(\hat{\theta}\) is our estimator of \(\theta\) - it’s a function of random variables
    • Probabilities reflect behavior in repeated samples
  • Bayesian inference takes a different approach to the problem. Rather than treating \(\theta\) as a constant, we consider it to be random as well!

    • Probabilities reflect beliefs about a particular quantity
    • \(f(\theta)\) denotes our prior beliefs about the value of \(\theta\)
    • \(f(\theta | \mathbf{Y})\) is the posterior distribution given the observed data. This is our target of inference
    • \(\mathbf{Y}\) are still random variables…but the posterior distribution conditions on them
  • We derive the posterior using Bayes rule

    \[f(\theta | \mathbf{Y}) = \frac{f(\mathbf{Y} | \theta)}{f(\mathbf{Y})} f(\theta)\]

Bayesian inference

  • The denominator can be written as an integral over all of the possible values of \(\theta\) (marginalizing over \(\theta\))

    \[f(\mathbf{Y}) = \int f(\mathbf{Y},\theta)d\theta = \int f(\mathbf{Y}|\theta)f(\theta)d\theta\]

  • So often we’ll write the posterior distribution up to a proportionality constant as

    \[\underbrace{f(\theta | \mathbf{Y})}_{\text{posterior}} \propto \underbrace{f(\mathbf{Y} | \theta)}_{\text{likelihood}} \times \underbrace{f(\theta)}_{\text{prior}}\]

Bayesian inference

  • In addition to the posterior distribution of the parameters \(\theta\), we’ll often want to generate “forecasts” of an out-of-sample \(\tilde{Y}\) conditional on what we have observed
    • This involves integrating over the posterior values of the parameter \(\theta\)
    \[\begin{align*}f(\tilde{Y}|\mathbf{Y}) &= \int f(\tilde{Y}, \theta | \mathbf{Y}) d\theta\\ &= \int f(\tilde{Y}| \theta, \mathbf{Y}) f(\theta | \mathbf{Y}) d\theta\\ &= \int f(\tilde{Y}| \theta) f(\theta | \mathbf{Y}) d\theta\end{align*}\]

Bayesian inference

  • We also might reason about two different values of \(\theta\): \(\{\theta_1, \theta_2\}\) and their relative odds. This has a nice expression in terms of the prior times a likelihood ratio

    \[\frac{f(\theta_1|\mathbf{Y})}{f(\theta_2|\mathbf{Y})} = \frac{f(\theta_1)}{f(\theta_2)} \times \frac{f(\mathbf{Y}|\theta_1)}{f(\mathbf{Y}|\theta_2)}\]

  • Inference on the posterior is often referred to as updating the prior and often can be conceptualized in terms of sequential changes to our beliefs about the parameter.

  • Suppose we have two samples \(\mathbf{Y}_1\) and \(\mathbf{Y}_2\). Our posterior can be written as:

    \[\begin{align*}f(\theta | \mathbf{Y}_1, \mathbf{Y}_2) &\propto f(\mathbf{Y}_2, \mathbf{Y}_1 | \theta) \times f(\theta) \\ &\propto f(\mathbf{Y}_2 | \mathbf{Y}_1, \theta) \times f(\mathbf{Y}_1 | \theta) \times f(\theta) \\ &\propto f(\mathbf{Y}_2 | \mathbf{Y}_1, \theta) \times f(\theta | \mathbf{Y}_1) \end{align*}\]

  • The posterior becomes the new “prior”!

Application: Absentee Ballot Rejection

  • A voter who returns an absentee ballot can still have it rejected

    • In Wisconsin, these decisions are made by municipal clerks
    • This has been the subject of some controversy when municipalities are overly strict
    • Reporting suggests that since 2024, officials in the city of Mequon have been using a particularly strict standard for adjudicating witness addresses (Shur 2026).
    • Let’s see if that is detectable in the data!
  • The EAC’s Election Administration and Voting Survey (EAVS) reports how many mail ballots were returned and rejected in every election jurisdiction (U.S. Election Assistance Commission 2025)

  • For each municipality \(i\), \(i \in \{1, 2, \dotsc, N\}\), we observe (November 2024 general election):

    • \(Y_i\): The number of returned absentee ballots that were rejected for discretionary reasons (i.e. not because the ballot arrived late or the voter died)
    • \(n_i\): The number of absentee ballots returned on time by living voters

Application: Absentee Ballot Rejection

  • We’ll start with a simple model
    • Assume \(Y_i\) is a draw from a binomial distribution with a total number of trials \(n_i\) and “rejection” probability \(\pi_i\).
    • \(\pi_i\) comes from a Beta distribution with hyperparameters \(\alpha_0\) and \(\beta_0\)
    \[Y_i \sim \text{Binomial}(n_i, \pi_i)\] \[\pi_i \sim \text{Beta}(\alpha_0, \beta_0)\]

Plate notation

  • A common method of writing statistical models is via plate notation
    • These encode independence and dependence assumptions across parameters and data.

Features of a Bayesian model

  • There are four types of variables in a Bayesian model
    • Observed Data: Variables that have a probability distribution but on whose observed values we condition - \(Y_i\)
    • Known Constants: Fixed quantities that do not have probability distributions (e.g. regressors or other features of the nodes) - \(n_i\), \(N\), \(\alpha_0\), \(\beta_0\)
    • Deterministic quantities: Transformations of other variables
    • Latent parameters: Variables that have a probability distribution that we do not observe - \(\pi_i\)

Features of a Bayesian model

  • Sometimes the known constants are actually known by us (e.g. \(N\) or \(n_i\)) and in other cases they are assumed to be known
    • These are typically parameters of the prior distribution which are called hyperparameters (here: \(\alpha_0\) and \(\beta_0\))
    • The hyperparameters govern the distribution of the prior
  • In Bayesian inference, we are interested in either obtaining the posterior of the latent parameters or integrating them out (in the latter case, they are sometimes referred to as “nuisance” parameters)
    • Here, we want to obtain \(f(\pi_i | \mathbf{Y})\)
    • When writing the posterior, we’ll often omit the implicit conditioning on the observed constants.

The Likelihood

  • Remember the posterior distribution is proportional to the likelihood times the prior

    \[\underbrace{f(\pi_i | \mathbf{Y})}_{\text{posterior}} \propto \underbrace{f(\mathbf{Y} | \pi_i)}_{\text{likelihood}} \times \underbrace{f(\pi_i)}_{\text{prior}}\]

  • First, let’s derive the likelihood \(f(\mathbf{Y} | \pi_i)\)

  • Since the \(Y_i\) are independent conditional on \(\pi_i\), we can write:

    \[f(\mathbf{Y} | \pi_i) = \prod_{j=1}^N f(Y_j | \pi_i)\]

  • Next, for \(j \neq i\), we have \(f(Y_j | \pi_i) = f(Y_j)\) since \(Y_j\) only depends on \(\pi_j\). That’s just a constant term and drops out of the posterior distribution, leaving

    \[f(\pi_i | \mathbf{Y}) = f(\pi_i | Y_i) \propto f(Y_i | \pi_i) f(\pi_i)\]

The Likelihood and Prior

  • What’s \(f(Y_i | \pi_i)\)? We’ve defined it as the binomial PMF

    \[f(Y_i | \pi_i) = {n_i \choose Y_i} \pi_i^{Y_i}(1- \pi_i)^{n_i - Y_i}\]

  • What about the prior? \(f(\pi_i)\)? We’ve chosen the beta distribution

    \[f(\pi_i) = \frac{\Gamma(\alpha_{0} + \beta_{0})}{\Gamma(\alpha_0)\Gamma(\beta_0)} \pi_i^{\alpha_0 - 1}(1- \pi_i)^{\beta_0 - 1}\]

  • The Beta distribution has some notable features.

    • Its mean: \(E[\pi_i] = \frac{\alpha_0}{\alpha_0 + \beta_0}\)
    • Its variance \(Var(\pi_i) = \frac{\alpha_0\beta_0}{(\alpha_0 + \beta_0)^2(\alpha_0 + \beta_0 + 1)}\)

The Beta Distribution

  • Let’s see how the parameters influence the shape of the beta. \(\text{Beta}(2,2)\):

The Beta Distribution

  • \(\text{Beta}(100,200)\):

Choosing a prior

  • There are two ways to think about choosing a prior distribution
    • Informative prior - Use the prior to encode our existing beliefs about the parameter
    • Uninformative/Diffuse prior - Pick a prior that will have the least impact on the posterior distribution
  • We also need to consider the shape of the prior distribution \(f(\pi_i)\)
    • Why did we pick the Beta distribution?
    • Because it has a special property when combined with a binomial likelihood.
    • It is a conjugate prior to the binomial likelihood.
  • Conjugate prior: A prior distribution is conjugate to a particular likelihood if the posterior distribution is of the same form as the prior
    • \(f(\theta)\) is beta; \(f(\mathbf{Y} | \theta)\) is binomial; \(f(\theta | \mathbf{Y})\) is beta

Choosing a prior

  • In this case, \(\pi_i\) must be between \(0\) and \(1\), so we already have some information.
  • One possible choice for an uninformative prior is the uniform distribution - each parameter value is equally likely: \[f(\pi_i) \propto 1\]
    • In settings where the parameter takes on values \((-\infty, \infty)\), we could still consider a “uniform” prior but it will be improper as it’s not a density that integrates to \(1\)
  • For \(\pi_i \in (0, 1)\), the uniform prior corresponds to the \(\text{Beta}(1,1)\) distribution
    • Originally proposed by Bayes.

Choosing a prior

  • Haldane proposed instead \(\text{Beta}(0,0)\) - notably this is improper \(\pi_i^{-1}(1-\pi_i)^{-1}\)

Choosing a prior

  • And Jeffreys proposed \(\text{Beta}(.5, .5)\)

The posterior density

  • We’ve mentioned that \(f(\pi_i | Y_i)\) is a beta distribution – let’s show that!

    \[f(\pi_i | \mathbf{Y}) \propto f(Y_i | \pi_i) f(\pi_i)\]

  • Plug in the densities (we’ll drop any multiplicative constants that don’t depend on \(\pi_i\))

    \[f(\pi_i | \mathbf{Y}) \propto \left[\pi_i^{Y_i}(1- \pi_i)^{n_i - Y_i}\right]\times\left[\pi_i^{\alpha_0 - 1}(1- \pi_i)^{\beta_0 - 1}\right]\]

  • Adding the exponents

    \[f(\pi_i | \mathbf{Y}) \propto \pi_i^{Y_i + \alpha_0 - 1}(1- \pi_i)^{n_i - Y_i + \beta_0 - 1}\]

  • And we can recognize this as the kernel of the beta distribution with parameters

    • \(\alpha = Y_i + \alpha_0\)
    • \(\beta = n_i - Y_i + \beta_0\)

Quantities of interest

  • Once we have a posterior distribution, we typically will report summaries

  • Point summaries

    • Posterior Mean: \(\hat{\theta} = E[\theta | \mathbf{Y}]\)
    • Posterior Mode: \(\hat{\theta} = \argmax_{\theta} f(\theta | \mathbf{Y})\)
  • Credible interval: A \(95\%\) credible interval \((l_{95}, h_{95})\) is a range of values that contains \(95\%\) of the posterior probability

    \[\int_{l_{95}}^{h_{95}} f(\theta | \mathbf{Y})d\theta = .95\]

  • There is no one unique credible interval! As a result, there are a few common choices for how to construct a credible interval

    • Highest Density interval (HDI) - no values outside of the interval have higher density than values inside the interval
    • Equal-tailed interval (ETI) - the probability of being below the lower limit is equal to the probability above the upper limit

Application

  • Let’s take a look at the EAVS data for Wisconsin. Start by loading it in (see data/README.md for the details)
# EAVS 2024, Wisconsin municipalities (C1b = mail ballots returned by voters, C9a = mail ballots rejected,
# C9b = rejected because late, C9m = rejected because the voter died)
wi <- read_csv("data/eavs_2024_wi.csv")
wi <- wi %>% mutate(n = C1b - C9b - C9m, y = C9a - C9b - C9m) %>% # Discretionary rejections only
  filter(n > 0) %>% # Drop municipalities where no absentee ballots were returned
  transmute(FIPSCode, municipality = Jurisdiction_Name, n, y, mle = y/n)
head(wi)
# A tibble: 6 × 5
  FIPSCode municipality                               n     y    mle
  <chr>    <chr>                                  <dbl> <dbl>  <dbl>
1 00100    CITY OF ABBOTSFORD - MULTIPLE COUNTIES    59     0 0     
2 00175    TOWN OF ABRAMS - OCONTO COUNTY           185     0 0     
3 00225    TOWN OF ACKLEY - LANGLADE COUNTY          30     2 0.0667
4 00275    CITY OF ADAMS - ADAMS COUNTY              70     0 0     
5 00300    TOWN OF ADAMS - ADAMS COUNTY             127     0 0     
6 00325    TOWN OF ADAMS - GREEN COUNTY              30     0 0     
  • The statewide rejection rate
state_rate <- sum(wi$y)/sum(wi$n)
state_rate
[1] 0.00208

Application

  • Our baseline approach – just take the mean everywhere. Which municipalities reject the most ballots?
wi %>% arrange(desc(mle)) %>% select(municipality, n, y, mle) %>% head(5)
# A tibble: 5 × 4
  municipality                            n     y   mle
  <chr>                               <dbl> <dbl> <dbl>
1 TOWN OF ALVIN - FOREST COUNTY           4     2 0.5  
2 VILLAGE OF LOWELL - DODGE COUNTY       18     3 0.167
3 TOWN OF VANCE CREEK - BARRON COUNTY     7     1 0.143
4 TOWN OF RIETBROCK - MARATHON COUNTY    29     4 0.138
5 TOWN OF PERSHING - TAYLOR COUNTY       15     2 0.133
  • All tiny towns!
  • And at the other end, 79% of municipalities have a rejection rate of exactly zero.

Application

  • The most extreme rejection rates come from the municipalities with the fewest ballots

Application

  • Consider the “worst” municipality in the state: the Town of Alvin in Forest County
alvin <- wi %>% filter(municipality == "TOWN OF ALVIN - FOREST COUNTY")
alvin
# A tibble: 1 × 5
  FIPSCode municipality                      n     y   mle
  <chr>    <chr>                         <dbl> <dbl> <dbl>
1 01625    TOWN OF ALVIN - FOREST COUNTY     4     2   0.5
  • 2 of 4 returned ballots were rejected
    • Do we really believe that the Town of Alvin rejects 50% of its absentee ballots?

Application

  • Let’s get the posterior distribution for \(\pi_i\) for the Town of Alvin under an uninformative uniform prior
posterior_alpha <- alvin$y + 1
posterior_beta <- alvin$n - alvin$y + 1
# Plot the posterior
alv_posterior <- ggplot() + xlim(0, 1)  + geom_function(fun=dbeta, args=list(shape1=posterior_alpha, shape2=posterior_beta), n=1001) + theme_bw() + xlab(expression(pi)) + 
  ylab("Density")
alv_posterior

Application

  • We know the form of the beta mean, so our posterior mean estimate is
posterior_mean <- posterior_alpha/(posterior_alpha + posterior_beta)
posterior_mean
[1] 0.5
  • And we can obtain an equal-tailed 95% credible interval simply via the quantile function
posterior_ci <- c(qbeta(.025, shape1=posterior_alpha, shape2=posterior_beta), qbeta(.975, shape1=posterior_alpha, shape2=posterior_beta))
posterior_ci
[1] 0.147 0.853

Application

  • But we have a lot more information from the other 1846 municipalities in the state
others_rate <- wi %>% filter(municipality != alvin$municipality) %>% summarize(rate = sum(y)/sum(n)) %>% pull(rate)
others_rate
[1] 0.00208
  • Our posterior mean under the uniform prior (0.5) is about 240 times the rate in the rest of the state!
    • The “uninformative” prior isn’t actually uninformative
    • It expects that municipalities reject half of all ballots!
  • What if we instead constructed our prior such that its mean was centered on the rejection rate in the rest of the state?

Application

  • We can interpret the hyperparameters as “previously observed” counts
    • \(\frac{\alpha_0}{\alpha_0 + \beta_0}\) is the prior mean – set it to the rejection rate in the rest of the state
    • \(\alpha_0 + \beta_0\) controls the strength of the prior – set it to \(500\): the prior is “worth” 500 ballots
prior_mean <- others_rate
prior_strength <- 500
# Convert these to alpha/beta 
prior_alpha <- prior_mean*prior_strength
prior_beta <- (1 - prior_mean)*prior_strength
c(prior_alpha, prior_beta)
[1]   1.04 498.96

Application

  • Let’s plug-in the new prior parameters
posterior2_alpha <- alvin$y + prior_alpha
posterior2_beta <- alvin$n - alvin$y + prior_beta
# Plot the posterior (cutting off the y-axis - the new posterior is much more concentrated)
alv_posterior2 <- ggplot() + xlim(0, 1)  +  geom_function(fun=dbeta, args=list(shape1=posterior_alpha, shape2=posterior_beta), n=1001, lty=2, lwd=1.5, col="dodgerblue") + 
  theme_bw() + xlab(expression(pi)) + 
  ylab("Density") + geom_vline(xintercept=prior_mean, lty=2, lwd=1.5, col="indianred") +
  geom_function(fun=dbeta, args=list(shape1=posterior2_alpha, shape2=posterior2_beta), n=1001,lwd=1.5) +
  coord_cartesian(ylim=c(0, 10))

Application

  • Uniform prior posterior (blue, dashed), prior mean (red, dashed), and statewide prior posterior (black)
alv_posterior2

Application

  • Our new posterior mean
posterior2_mean <- posterior2_alpha/(posterior2_alpha + posterior2_beta)
posterior2_mean
[1] 0.00603
  • And 95% credible interval
posterior2_ci <- c(qbeta(.025, shape1=posterior2_alpha, shape2=posterior2_beta), qbeta(.975, shape1=posterior2_alpha, shape2=posterior2_beta))
posterior2_ci
[1] 0.00127 0.01442
  • With only 4 ballots, the prior dominates. Alvin’s posterior mean is much closer to the statewide rate than to the MLE: 50%.

Application

  • Let’s do this for every municipality: pull each estimate toward the statewide rate
wi <- wi %>% mutate(post_mean = (y + prior_strength*state_rate)/(n + prior_strength))
  • Now who rejects the most ballots?
wi %>% arrange(desc(post_mean)) %>% select(municipality, n, y, mle, post_mean) %>% head(5)
# A tibble: 5 × 5
  municipality                          n     y    mle post_mean
  <chr>                             <dbl> <dbl>  <dbl>     <dbl>
1 TOWN OF DOVER - RACINE COUNTY       297     9 0.0303    0.0126
2 TOWN OF UNION - EAU CLAIRE COUNTY   313     9 0.0288    0.0124
3 TOWN OF CHETEK - BARRON COUNTY      205     7 0.0341    0.0114
4 CITY OF PARK FALLS - PRICE COUNTY   139     6 0.0432    0.0110
5 TOWN OF ALMA - JACKSON COUNTY        66     5 0.0758    0.0107
  • The top of the list now includes places with more ballots – a high rate on a lot of ballots is more convincing than 2 out of 4

Application

  • The MLE vs. the posterior mean for every municipality (note the different y-axes)
  • Small municipalities get pulled almost all the way to the statewide rate while big cities barely move.

Does pooling help?

  • How do we know whether this is a good idea? We also have EAVS data from the November 2022 election (U.S. Election Assistance Commission 2023)
    • Use each municipality’s 2024 data to predict its 2022 rejections
    • The prior only matters when there’s very little data, so focus on the municipalities with fewer than 200 returned ballots in 2024
wi_2022 <- read_csv("data/eavs_2022_wi.csv") %>%
  mutate(n_2022 = C1b - C9b - C9m, y_2022 = C9a - C9b - C9m) %>% filter(n_2022 > 0) %>% # Discretionary rejections only
  transmute(FIPSCode, n_2022, y_2022)
oos <- wi %>% filter(n < 200) %>% inner_join(wi_2022, by="FIPSCode") # Municipalities in both years
scale_2022 <- (sum(wi_2022$y_2022)/sum(wi_2022$n_2022))/state_rate # Statewide 2022 rate relative to 2024
nrow(oos)
[1] 1418

Does pooling help?

  • Score each set of 2024 estimates by the log-likelihood of the 2022 data (higher is better)
    • Rejections were 2.5 times as common in 2022, so scale every prediction up by that ratio
# Truncate predictions at .0001 so that a predicted probability of 0 doesn't give log(0)
pred_loglik <- function(p) sum(dbinom(oos$y_2022, oos$n_2022, prob = pmin(pmax(p*scale_2022, .0001), .9999), log=T))
strengths <- c(50, 200, 500, 1000)
pull_loglik <- sapply(strengths, function(k) pred_loglik((oos$y + k*state_rate)/(oos$n + k)))
best_k <- strengths[which.max(pull_loglik)] # The best-predicting prior strength (we'll map it below)
tibble(estimator = c("MLE", "Uniform prior", "Statewide rate only", paste("Pull to statewide rate, prior strength", strengths)),
       loglik = c(pred_loglik(oos$mle), pred_loglik((oos$y + 1)/(oos$n + 2)), pred_loglik(rep(state_rate, nrow(oos))), pull_loglik))
# A tibble: 7 × 2
  estimator                                   loglik
  <chr>                                        <dbl>
1 MLE                                         -2490.
2 Uniform prior                               -2824.
3 Statewide rate only                         -1249.
4 Pull to statewide rate, prior strength 50   -1392.
5 Pull to statewide rate, prior strength 200  -1222.
6 Pull to statewide rate, prior strength 500  -1202.
7 Pull to statewide rate, prior strength 1000 -1210.

Does pooling help?

  • How did the municipalities with the most extreme 2024 rates do in 2022?
oos %>% mutate(group = case_when(rank(-mle, ties.method = "first") <= 50 ~ "Top 50 in 2024",
                                 y == 0 ~ "No rejections in 2024", TRUE ~ "Everyone else")) %>%
  group_by(group) %>%
  summarize(municipalities = n(), times_state_2024 = (sum(y)/sum(n))/state_rate,
            times_state_2022 = (sum(y_2022)/sum(n_2022))/(state_rate*scale_2022))
# A tibble: 3 × 4
  group                 municipalities times_state_2024 times_state_2022
  <chr>                          <int>            <dbl>            <dbl>
1 Everyone else                    113             6.19             2.60
2 No rejections in 2024           1255             0                1.38
3 Top 50 in 2024                    50            23.9              6.00
  • And the City of Mequon? 4.1 times the statewide rate in 2024, but 0.5 times the statewide rate in 2022.

Does pooling help?

  • The MLE does terribly for forecasting small areas - the extreme rates of 2024 were mostly noise

    • The “worst” municipalities of 2024 rejected about 6 times the statewide rate in 2022, not 24 times
    • The municipalities with no rejections in 2024 weren’t at zero in 2022
  • The uniform prior does even worse - it pulls toward rejecting half of all ballots!

  • The statewide rate alone – literally just guessing the mean – does much better

  • But a combination of the statewide rate and the local rate does best

    • There’s some additional information in a municipality’s own history: the top 50 still rejected well above the statewide rate in 2022
  • We directly manipulated the strength of the prior. Could we instead estimate it from the data?

    • Yes! This is the idea behind hierarchical (or multilevel) models - we’ll cover this in a few weeks.

Mapping rejection rates

  • Let’s map every municipality: the MLE vs. the posterior mean under our “best” prior (prior strength 500)
library(sf)
wi_map <- st_read("data/wi_municipalities_2024.gpkg", quiet = T)
wi_est <- wi %>% mutate(post_best = (y + best_k*state_rate)/(n + best_k))
wi_map <- wi_map %>% left_join(wi_est, by = c("COUSUBFP" = "FIPSCode")) %>%
  pivot_longer(c(mle, post_best), names_to = "estimate", values_to = "rate") %>%
  mutate(estimate = factor(estimate, levels = c("mle", "post_best"),
                           labels = c("MLE", paste0("Posterior mean (prior strength ", best_k, ")"))))

Mapping rejection rates

  • Note that maps of shrunken estimates can mislead too - the highest adjusted rates end up in more populated places (Gelman and Price 1999)

Application

  • The posterior mean can be thought of as a “weighted average” of the prior mean and the MLE

    • The weights are controlled by the variance of the prior.

    \[E[\pi_i | Y_i] = \frac{Y_i + \alpha_0}{n_i + \alpha_0 + \beta_0} = \underbrace{\frac{n_i}{n_i + \alpha_0 + \beta_0}}_{\text{weight on data}} \times \frac{Y_i}{n_i} + \underbrace{\frac{\alpha_0 + \beta_0}{n_i + \alpha_0 + \beta_0}}_{\text{weight on prior}} \times \frac{\alpha_0}{\alpha_0 + \beta_0}\]

    • Can interpret the hyper-parameters for this case - \(\alpha_0\) and \(\beta_0\) - as the number of “previously observed” counts.
  • Stronger priors \(\leadsto\) narrower credible intervals

    • But it takes a lot more data to move the posterior distribution from the prior.
  • Often, the prior serves to regularize our estimates

    • We want our estimates to be “pulled” towards a particular value if there’s very little data.
    • A common type of “regularizing” prior is designed to attenuate our estimates to \(0\).

Next week

  • Estimating Bayesian models
    • Most posteriors don’t have a nice closed form
      • Beta-binomial is the classic conjugate prior.
      • Normal-Normal is also common (you’ll see a variation on this on the problem set)
    • Markov Chain Monte Carlo (MCMC): drawing samples from a posterior that we can only evaluate up to a proportionality constant
    • The Metropolis-Hastings algorithm and Gibbs sampling
  • Bayesian regression
    • Putting priors on regression coefficients and the error variance
  • Introduction to Stan
    • Write down the model – let Stan implement the MCMC sampler
  • Readings: Chapters 3 and 9 of McElreath (2020)

References

Gelman, Andrew, and Phillip N. Price. 1999. “All Maps of Parameter Estimates Are Misleading.” Statistics in Medicine 18 (23): 3221–34. https://doi.org/10.1002/(SICI)1097-0258(19991215)18:23<3221::AID-SIM312>3.0.CO;2-M.
McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. 2nd ed. Chapman and Hall/CRC. https://xcelab.net/rm/.
Shur, Alexander. 2026. How Inconsistent Standards Led to Dozens of Disenfranchised Voters in Mequon. Votebeat; republished by Wisconsin Watch. https://wisconsinwatch.org/2026/07/wisconsin-election-mequon-disenfranchised-voters-absentee-ballot-envelope-witness-addresses/.
U.S. Election Assistance Commission. 2023. Election Administration and Voting Survey 2022 (Public Release Dataset, V1). Https://www.eac.gov/research-and-data/datasets-codebooks-and-surveys.
U.S. Election Assistance Commission. 2025. Election Administration and Voting Survey 2024 (Public Release Dataset, V1). Https://www.eac.gov/research-and-data/datasets-codebooks-and-surveys.