Bayesian Statistics in R

Reviewed & published by Brayan K

Bayesian statistics updates beliefs with evidence: start from a prior, combine it with the likelihood of the data, and get a posterior distribution that quantifies your remaining uncertainty.

Part of the free R course at LearnCodingFast — hands-on lessons with worked examples and the output they print, plus practice exercises and a quick quiz.

By the end of this lesson you'll understand priors, likelihoods, and posteriors, fit Bayesian regressions with brms and rstanarm via MCMC, and read credible intervals and posterior predictive checks.

What You'll Learn in This Lesson

• What priors, likelihoods, and posteriors are
• How Bayes' theorem updates belief with data
• How to fit Bayesian regression with brms / rstanarm
• What MCMC sampling does and why it is needed
• How credible intervals differ from confidence intervals
• How posterior predictive checks assess model fit

1️⃣ Bayes' Theorem: Prior x Likelihood

The posterior is proportional to the prior times the likelihood. With a conjugate prior the math is clean: a Beta prior plus Binomial data gives a Beta posterior, no sampling required.

# Bayes' theorem in action: updating belief that a coin is fair.
# Posterior  ~  Prior  x  Likelihood

# Suppose a coin lands heads 7 of 10 flips.
# A Beta(2, 2) prior says "probably near 0.5, but unsure".
# The Beta is conjugate to the Binomial, so the posterior is also Beta:
#   posterior = Beta(prior_a + heads, prior_b + tails)

heads <- 7; tails <- 3
post_a <- 2 + heads   # 9
post_b <- 2 + tails   # 5

# Posterior mean for p (probability of heads):
post_a / (post_a + post_b)   # -> [1] 0.6428571

2️⃣ Credible Intervals

A credible interval is a direct probability statement about a parameter. From the Beta(9, 5) posterior, qbeta() gives the 2.5% and 97.5% quantiles — the 95% credible interval.

# A 95% CREDIBLE interval: a range that, given the model,
# contains the parameter with 95% posterior probability.

# Using the Beta(9, 5) posterior from before:
qbeta(c(0.025, 0.975), 9, 5)
# -> [1] 0.3857 0.8614

# Interpretation: there is a 95% probability that p lies in
# roughly [0.39, 0.86]. That is a DIRECT statement about p --
# unlike a frequentist confidence interval, which is about the
# long-run behaviour of the procedure, not this one interval.

3️⃣ Bayesian Regression with brms

brms writes Stan code for you. You give a formula and priors; it runs MCMC chains and returns posterior means, credible intervals, and the Rhat convergence diagnostic.

library(brms)

# Fit a Bayesian linear regression with brms (built on Stan).
# brms draws from the posterior with MCMC sampling.
fit <- brm(
  mpg ~ wt + hp,        # outcome ~ predictors
  data    = mtcars,
  prior   = prior(normal(0, 5), class = "b"),  # weakly informative prior on slopes
  chains  = 4,          # 4 independent MCMC chains
  iter    = 2000,       # iterations per chain
  seed    = 1
)

summary(fit)   # posterior means, credible intervals, and Rhat

4️⃣ Posterior Predictive Checks

A fitted model should reproduce the data it was trained on. pp_check() simulates datasets from the posterior and overlays them on the observed data so you can spot misfit at a glance.

library(brms)

# Posterior predictive check: simulate datasets from the fitted model
# and overlay them on the observed data. If the model fits well,
# the simulated curves (yrep) should look like the observed one (y).
pp_check(fit)

# rstanarm offers a similar workflow with stan_glm():
#   library(rstanarm)
#   fit2 <- stan_glm(mpg ~ wt + hp, data = mtcars)
#   posterior_interval(fit2, prob = 0.95)   # credible intervals

Your turn. Fill in the ___ blank, run it, and compare with the expected output.

# 🎯 YOUR TURN — replace each ___ then run it.

# A drug helped 8 of 12 patients. Start from a uniform Beta(1, 1) prior.
heads <- 8        # successes
tails <- 12 - 8   # failures

# 1) Beta is conjugate to the Binomial, so update the parameters:
post_a <- 1 + heads
post_b <- 1 + ___   # 👉 replace ___ with tails

# 2) Posterior mean for the success rate
cat("Posterior mean:", round(post_a / (post_a + post_b), 3), "\n")

# 3) 95% credible interval
print(round(qbeta(c(0.025, 0.975), post_a, post_b), 3))

# ✅ Expected output:
# Posterior mean: 0.643
# [1] 0.386 0.861

Mini-Challenge: Prior to Posterior

Update a Beta prior with spam-filter data and report the posterior mean and a 90% credible interval. The conjugate Beta-Binomial pair keeps the arithmetic simple.

# 🎯 MINI-CHALLENGE: prior to posterior, by hand.
#
# An email filter flags 30 of 50 messages as spam.
# 1. Start from a Beta(2, 2) prior on the spam rate.
# 2. Update to the posterior Beta(2 + spam, 2 + ham).
# 3. Print the posterior mean.
# 4. Print the 90% credible interval with
#    qbeta(c(0.05, 0.95), post_a, post_b).
#
# ✅ Example output:
#    Posterior mean: 0.593
#    90% CI: 0.481 0.700

# your code here

Common Errors (and the fix)

• Misreading a confidence interval as a credible one — only the Bayesian credible interval makes a direct probability statement about the parameter.
• Ignoring high Rhat — an Rhat well above 1.0 means the chains did not converge; run more iterations or rethink the model.
• Overly tight priors — a strong prior can dominate small datasets. Prefer weakly informative priors unless you truly have prior knowledge.
• "could not find function 'brm'" — load the package with library(brms) (and ensure Stan is installed).

Pro Tips

• 💡 Start weakly informative: priors like normal(0, 5) on slopes regularize without overpowering the data.
• 💡 Always check chains: look at Rhat and effective sample size before trusting estimates.
• 💡 Visualize the posterior: plot(fit) and pp_check(fit) tell you more than a single number.

📋 Quick Reference — Bayesian R

🎉 Lesson Complete!

• ✅ Posterior ∝ prior × likelihood (Bayes' theorem)
• ✅ Conjugate priors give a posterior in closed form
• ✅ brms and rstanarm fit Bayesian regression via MCMC
• ✅ Credible intervals make direct probability statements
• ✅ pp_check() and Rhat assess fit and convergence

Practice quiz

Continue this course

• Previous: Machine Learning with caret
• Next: Cross-Validation & ROC — k-fold CV, overfitting, and the ROC curve / AUC for honest model evaluation
• Quick reference: R cheat sheet

Frequently asked questions

What is the core idea of Bayesian statistics?

Bayesian statistics treats unknown quantities as having probability distributions that express our uncertainty. You start with a prior distribution — what you believe before seeing the data — then update it with the likelihood, which measures how well each parameter value explains the observed data. The result is the posterior distribution, your updated belief. Bayes' theorem is the engine: the posterior is proportional to the prior times the likelihood. Everything you report — point estimates, intervals, predictions — comes from summarizing that posterior.

How does a credible interval differ from a confidence interval?

A 95% credible interval is a range that, according to your model, has a 95% posterior probability of containing the parameter — a direct probability statement about the unknown quantity. A 95% confidence interval is subtler: it is a property of the procedure, meaning that if you repeated the experiment many times, about 95% of the intervals constructed this way would contain the true value. People often interpret confidence intervals as if they were credible intervals; Bayesian inference lets you make that intuitive statement legitimately, at the cost of specifying a prior.

What are brms and rstanarm, and why use MCMC?

brms and rstanarm are R packages that fit Bayesian regression models using Stan, a probabilistic programming engine. You write a familiar formula like mpg ~ wt + hp, and they translate it into a Stan model and sample from the posterior. Most posteriors have no closed-form solution, so these tools use Markov chain Monte Carlo (MCMC) — specifically Hamiltonian Monte Carlo — to draw samples whose distribution approximates the posterior. From those draws you compute posterior means, credible intervals, and predictions.

How do I know the model fit and the sampler worked?

Two checks matter. First, sampler diagnostics: the Rhat statistic should be close to 1.00 for every parameter, indicating the independent chains converged to the same distribution, and the effective sample size should be reasonably large. Warnings about divergent transitions signal trouble. Second, posterior predictive checks: simulate new datasets from the fitted model with pp_check() and compare them to the observed data. If the simulated data look like the real data, the model captures its important features; systematic mismatches point to a model that needs revising.