The Bayesian Workflow in R: All 5 Steps Explained

The Bayesian workflow is the disciplined loop that turns a single Bayesian model into a trustworthy answer: you set priors, fit the model, check whether the fit could have produced your real data, and revise when a check fails. It is iterative on purpose. This tutorial walks all five steps on one dataset in R, so by the end you can run the whole loop yourself, from prior to report.

What is the Bayesian workflow, and why is it a loop?

Most tutorials teach you to fit one Bayesian model and read off the numbers. Real analysis never works that way. Andrew Gelman and colleagues describe the honest practice as a loop: you propose a model, pit it against the data, then improve it wherever it falls short. The Bayesian workflow is that loop, done deliberately instead of by accident.

We will follow one question the whole way through. In the built-in mtcars data, does a car's weight affect its fuel economy differently for manual versus automatic transmissions? Let us start by centering weight at its average (so the intercept means "mpg for a typical-weight car") and looking at the raw slope of miles-per-gallon against weight inside each transmission group.

Read the two numbers. For automatic cars, every extra 1,000 pounds costs about 3.8 mpg. For manual cars, the same weight costs about 9.1 mpg, nearly two and a half times as much. That gap is the puzzle the whole workflow will test: is it a real difference the model must capture, or just noise from 32 cars? A single fit cannot answer that. The loop can.

Here are the five steps we will run, and how they connect.

Figure 1: The five workflow steps form a loop: a failed check sends you back to revise and refit.

1. Set priors and check them by simulating fake data before you look at the real data.
2. Fit the model to the data.
3. Diagnose the sampler to make sure the computation actually worked.
4. Check the fitted model against the data, looking for places it cannot reproduce.
5. Revise the model when a check fails, then run the loop again.

The arrow from step 5 back to step 2 is the whole point. A Bayesian analysis is finished not when the model runs, but when it survives the checks. Throughout, the model fitting itself uses the rstanarm package (which needs a local R session), while every prior and posterior check runs right here in your browser using base R.

Step 1: How do you set priors and check them before seeing data?

A prior is your belief about a parameter before you see the data, written as a probability distribution. A weakly-informative prior is one that rules out absurd values but stays open-minded about anything reasonable. For our model, we need priors for the intercept (average mpg at average weight) and the weight slopes, plus a prior for the size of the noise. Good starting choices are a Normal centered at 20 for the intercept, a Normal centered at 0 for each slope, and a positive prior for the noise size.

How do you know a prior is reasonable? You simulate. A prior predictive check draws parameter values from the priors alone, generates the fake datasets those priors imply, and asks a simple question: could this data plausibly exist? Let us draw 1,000 sets of parameters from weakly-informative priors and see what range of fuel economy they imply across the real weights.

The middle 90% of prior-implied mpg values runs from about -2 to 42. That is wide, which is the point of a weakly-informative prior: it does not pretend to know the answer. About 6% of the implied values dip below zero, which is impossible for real fuel economy. A small amount of leakage into impossible territory is acceptable here, because the data we are about to add will easily pull the estimates into a sensible range.

Now watch what a careless prior does. If we make the slope prior vague, a Normal centered at 0 with a standard deviation of 100, the implied data becomes nonsense.

Under the vague prior, the implied mpg swings from -127 to 167, and nearly a third of the values are physically impossible. A prior that expects cars getting -127 mpg is not "letting the data speak"; it is asserting a belief that most possible worlds are absurd. The prior predictive check catches this before any fitting happens, which is exactly when it is cheapest to fix.

Step 2: How do you fit the model and read it honestly?

We will start with the simplest model that could answer our question: mpg depends on weight and on transmission, but with a single common weight slope for both groups. In formula terms that is mpg ~ wt_c + am. We fit it with stan_glm(), which draws thousands of samples from the posterior distribution, the updated belief about each parameter after seeing the data.

$$ \text{mpg}_i = \beta_0 + \beta_{\text{wt}}\, \text{wt\_c}_i + \beta_{\text{am}}\, \text{manual}_i + \varepsilon_i $$

This additive form gives manual cars a different intercept (through the am term) but forces both transmission types to share the same weight slope. Let us fit it and read the summary.

library(rstanarm)
options(mc.cores = 1)

fit_add <- stan_glm(mpg ~ wt_c + am, data = mt,
                    prior = normal(0, 10), prior_intercept = normal(20, 10),
                    prior_aux = exponential(0.2),
                    seed = 2027, refresh = 0)
print(fit_add)
#> stan_glm
#>  family:       gaussian [identity]
#>  formula:      mpg ~ wt_c + am
#>  observations: 32
#>  predictors:   3
#> ------
#>             Median MAD_SD
#> (Intercept) 20.1    0.8
#> wt_c        -5.3    0.8
#> ammanual     0.0    1.6
#>
#> Auxiliary parameter(s):
#>       Median MAD_SD
#> sigma 3.1    0.4

Read this honestly. The weight slope is about -5.3 mpg per 1,000 pounds, tightly estimated. But look at the transmission term, ammanual: its posterior median is 0.0 with a spread of 1.6. Once the model already knows a car's weight, adding "is it manual?" does essentially nothing. The coefficient table shows that once weight is in the model, transmission does not matter.

A single number is not a Bayesian answer. The whole appeal of the Bayesian approach is that it gives you a full distribution for each parameter, so you can state a credible interval: a range that contains the parameter with a stated probability. Let us draw from the fitted posterior and read off the 90% credible interval for the transmission effect. The four numbers in mu_a are the posterior means of the three coefficients plus log sigma, and S_a is their covariance matrix; both are summaries of the local fit above.

The 90% credible interval for the transmission effect runs from about -2.5 to +2.5 mpg. It straddles zero, and it is roughly symmetric around zero. In plain language: after accounting for weight, this model sees no evidence that transmission changes fuel economy in either direction. If you stopped here, you would report "transmission does not matter" and move on. Hold that thought, because it is wrong, and the workflow is about to show you why.

Step 3: Did the sampler actually converge?

Before you trust a single number from a fitted model, you have to know the computation worked. stan_glm() does not solve the posterior with algebra; it explores it with a Markov chain that wanders through parameter space. If that wandering has not settled down, the summaries are garbage no matter how sensible they look.

Two numbers tell you whether it settled. R-hat compares the several chains the sampler runs in parallel: if they have all converged to the same distribution, R-hat sits at 1.00, and anything above about 1.01 is a warning. The effective sample size (n_eff) estimates how many truly independent draws you have; you want it in the hundreds or thousands, not the tens. Here they are for our additive fit.

diag <- summary(fit_add)[c("(Intercept)", "wt_c", "ammanual", "sigma"),
                         c("Rhat", "n_eff")]
round(diag, 2)
#>             Rhat n_eff
#> (Intercept)    1  2223
#> wt_c           1  2133
#> ammanual       1  2268
#> sigma          1  2413

Every R-hat is 1.00 and every effective sample size is over 2,000. The sampler converged and gave us plenty of independent information about each parameter. This is the step people skip, and skipping it is how wrong conclusions get published. For the full battery of convergence diagnostics, including trace plots that show the chains overlapping like a fuzzy caterpillar, see the companion tutorial on MCMC diagnostics.

Step 4: Does the model reproduce the data?

Here is the heart of the workflow. A model can converge perfectly and still be wrong, because converging only means the computation matched the model you specified, not that the model matches reality. A posterior predictive check closes that gap. It asks: if this fitted model is true, what data would it generate, and does that fake data look like the data we actually saw?

The recipe is simple. Draw parameters from the posterior, use them to simulate a replicated dataset the same size as the real one, compute some summary statistic on that replicate, and repeat thousands of times. Then compare the real data's statistic to the cloud of replicated statistics. Pick a statistic that targets what you care about. Our question is about group-specific slopes, so the natural statistic is the gap between the manual weight slope and the automatic weight slope.

$$ T(y) = \text{slope}_{\text{manual}} - \text{slope}_{\text{automatic}} $$

We already know the observed gap is about -5.3 (manuals lose mpg roughly 5.3 faster per 1,000 pounds). Let us simulate that gap from the additive model many times and see where the observed value falls.

Look at what happened. The additive model replicates slope gaps centered on zero, because it forces both groups to share one slope, so any gap it produces is pure noise. The observed gap of -5.3 sits far out in the left tail. The fraction of replicates as extreme as the data, a posterior predictive p-value, is 0.002. That is a clear failure: the model almost never generates data with a slope gap like the one we actually see.

Figure 2: A posterior predictive check compares the data to what the fitted model would generate, and the verdict decides whether you revise.

This is the decision point of the loop. The coefficient table in step 2 said transmission does not matter. The check says the model that assumes transmission does not matter cannot reproduce the data. When the summary and the check disagree, the check wins, because the check is looking at the actual data-generating behavior of the model, not just one parameter in isolation.

Step 5: How do you revise the model when a check fails?

A failed check is not a dead end; it is an instruction. The check told us the additive model cannot produce different slopes for different transmissions, so we revise the model to allow exactly that. We add an interaction between weight and transmission, written mpg ~ wt_c * am. The interaction term lets each group have its own weight slope.

$$ \text{mpg}_i = \beta_0 + \beta_{\text{wt}}\, \text{wt\_c}_i + \beta_{\text{am}}\, \text{manual}_i + \beta_{\text{wt:am}}\, (\text{wt\_c}_i \times \text{manual}_i) + \varepsilon_i $$

The new coefficient, the one multiplying weight-times-manual, is the extra slope manual cars get on top of the automatic slope. Let us refit and read it.

fit_int <- stan_glm(mpg ~ wt_c * am, data = mt,
                    prior = normal(0, 10), prior_intercept = normal(20, 10),
                    prior_aux = exponential(0.2),
                    seed = 2027, refresh = 0)
print(fit_int)
#> stan_glm
#>  family:       gaussian [identity]
#>  formula:      mpg ~ wt_c * am
#>  observations: 32
#>  predictors:   4
#> ------
#>               Median MAD_SD
#> (Intercept)   19.2    0.7
#> wt_c          -3.8    0.8
#> ammanual      -2.1    1.4
#> wt_c:ammanual -5.2    1.5
#>
#> Auxiliary parameter(s):
#>       Median MAD_SD
#> sigma 2.6    0.4

Now the story changes. The automatic weight slope is about -3.8, and the interaction term wt_c:ammanual is about -5.2, well separated from zero given its spread of 1.5. Manual cars lose roughly 5.2 mpg more per 1,000 pounds than automatics, which recovers the raw -3.8-versus-9.1 gap we saw at the very start. The transmission effect the additive model missed was never about the intercept; it was about the slope.

But a smaller residual noise (sigma dropped from 3.1 to 2.6) and a plausible-looking coefficient are not proof. We revised because a check failed, so we are only done when the same check passes. Let us rerun the identical posterior predictive check on the interaction model.

The replicated slope gaps now center on -5.26, right where the observed -5.3 lives, and the posterior predictive p-value is 0.489, almost exactly one half. The observed data is now completely typical of what the model generates. The revision worked, and we can prove it with the same check that rejected the old model. That is the loop closing: a failed check drove a specific change, and the change fixed the specific failure.

How sensitive are your conclusions to the prior?

You chose those priors. A careful reader will wonder whether your conclusion is a fact about the data or an artifact of that choice. A sensitivity analysis answers the question directly: refit the model under a much wider prior and a much tighter one, and see whether the conclusion moves. If the answer barely changes, the data is doing the talking. If it swings, you must report that your result depends on your assumptions.

We refit the interaction model three ways: a skeptical prior that strongly doubts large effects (Normal centered at 0 with standard deviation 2), our weakly-informative baseline (standard deviation 10), and a nearly flat prior (standard deviation 50). Then we compare the posterior for the interaction term.

fit_tight <- stan_glm(mpg ~ wt_c * am, data = mt, prior = normal(0, 2),
                      prior_intercept = normal(20, 10), prior_aux = exponential(0.2),
                      seed = 2027, refresh = 0)
fit_wide  <- stan_glm(mpg ~ wt_c * am, data = mt, prior = normal(0, 50),
                      prior_intercept = normal(20, 10), prior_aux = exponential(0.2),
                      seed = 2027, refresh = 0)

pick <- function(f) {
  x <- as.matrix(f)[, "wt_c:ammanual"]
  round(c(median = median(x), q05 = quantile(x, 0.05), q95 = quantile(x, 0.95)), 2)
}
rbind(tight = pick(fit_tight), base = pick(fit_int), wide = pick(fit_wide))
#>       median q05.5% q95.95%
#> tight  -3.67  -5.60   -1.64
#> base   -5.24  -7.73   -2.86
#> wide   -5.28  -7.73   -2.78

Read the three rows. Under the skeptical prior the interaction shrinks toward zero, from -5.2 down to -3.7, exactly what a skeptical prior is supposed to do. But even then its 90% credible interval, -5.6 to -1.6, stays entirely below zero. The wide prior is almost identical to the baseline. Across an aggressive range of prior beliefs, the sign and the significance of the effect never change; only the magnitude wobbles a little under the most skeptical prior. That is a conclusion that holds up, and being able to say so honestly is worth more than any single point estimate.

How do you report a Bayesian analysis honestly?

An honest Bayesian write-up shows its work at every step of the loop. It states the priors and why they are reasonable, confirms the sampler converged, presents the posterior predictive checks that the final model passes, reports a sensitivity analysis, and gives uncertainty as intervals rather than bare point estimates. A reader should be able to see not just your answer but every place the model could have failed and did not.

One more tool belongs in the report: a fair comparison of the models you considered. Leave-one-out cross-validation (LOO) estimates how well each model predicts data it has not seen, which is a better basis for comparison than in-sample fit. It is a large topic; here we use it only to confirm that the revision earned its extra complexity.

loo_add <- loo(fit_add)
loo_int <- loo(fit_int)
loo_compare(loo_add, loo_int)
#>         elpd_diff se_diff
#> fit_int  0.0       0.0
#> fit_add -5.5       2.4

The interaction model (fit_int) sits at the top, and the additive model predicts about 5.5 units worse, with a standard error of 2.4 on that difference. The gap is a little over two standard errors, so LOO agrees with the predictive check: the interaction model is the better description, and its extra parameter pays for itself. With only 32 cars this comparison is suggestive rather than decisive, which is itself worth stating in a report.

The complete workflow in one script

Here is the entire loop in one place: fit the simple model, check it, revise when the check fails, then confirm the fix works. This block uses the real posterior draws from rstanarm (through posterior_predict()), so it complements the browser-friendly approximations above and gives the same verdict. Run it in a local R session with rstanarm installed.

library(rstanarm)
cars <- mtcars
cars$am   <- factor(cars$am, labels = c("automatic", "manual"))
cars$wt_c <- cars$wt - mean(cars$wt)

gap <- function(y, d) {
  a <- coef(lm(y[d$am == "automatic"] ~ d$wt_c[d$am == "automatic"]))[2]
  m <- coef(lm(y[d$am == "manual"]    ~ d$wt_c[d$am == "manual"]))[2]
  as.numeric(m - a)
}
obs <- gap(cars$mpg, cars)

# Step 1-2: priors + fit the simple additive model
m_add <- stan_glm(mpg ~ wt_c + am, data = cars, prior = normal(0, 10),
                  prior_intercept = normal(20, 10), prior_aux = exponential(0.2),
                  seed = 2027, refresh = 0)

# Step 4: posterior predictive check on the group-slope gap
set.seed(5)
rep_add <- apply(posterior_predict(m_add, draws = 2000), 1, gap, d = cars)
cat("additive model check p-value:", round(mean(rep_add <= obs), 2), "\n")
#> additive model check p-value: 0

# Step 5: revise (add interaction) and re-check
m_int <- stan_glm(mpg ~ wt_c * am, data = cars, prior = normal(0, 10),
                  prior_intercept = normal(20, 10), prior_aux = exponential(0.2),
                  seed = 2027, refresh = 0)
set.seed(5)
rep_int <- apply(posterior_predict(m_int, draws = 2000), 1, gap, d = cars)
cat("revised model check p-value:", round(mean(rep_int <= obs), 2), "\n")
#> revised model check p-value: 0.49

The simple model fails the check with a p-value of essentially 0; the revised model passes with a p-value near 0.5. That is the workflow in miniature: propose, check, revise, confirm.

Practice Exercises

These combine several steps of the workflow. Each runs in your browser using the posterior draws you built above (tha for the additive model, thi for the interaction model). Use fresh variable names so you do not overwrite the tutorial's objects.

Exercise 1: Summarize and probe the interaction effect

Using the interaction draws in thi, report the 90% credible interval for the extra weight penalty manual cars pay (column 4, wt_c:ammanual), and the posterior probability that this penalty is steeper than -3 mpg per 1,000 pounds.

Exercise 2: Build a reusable posterior predictive check

Write a function that runs a posterior predictive check for any test statistic, then use it to check whether the interaction model reproduces the maximum mpg in the data. A good model should generate replicates whose maximum is similar to the observed maximum.

Exercise 3: Interpret the effect at an average-weight car

The interaction changed the story about transmission. Using thi, compute the posterior probability that a manual car gets higher mpg than an automatic car of average weight (weight centered at zero). At that point the difference is just the ammanual term, column 3. Then say in one sentence what this means for the old claim that "manuals get better mileage."

Frequently Asked Questions

Is the Bayesian workflow the same as just fitting a Bayesian model?

No. Fitting is only step 2. The workflow is the full loop of setting and checking priors, fitting, diagnosing the sampler, checking predictions against data, and revising. The discipline lives in the checking and revising, not in the fit.

What is the difference between a prior predictive check and a posterior predictive check?

A prior predictive check simulates data from the priors alone, before fitting, to catch absurd assumptions early. A posterior predictive check simulates data from the fitted model, after seeing the data, to find where the model fails to reproduce reality. One guards the input, the other tests the output.

Is a posterior predictive p-value the same as a frequentist p-value?

No. A posterior predictive p-value measures how typical your observed data is among datasets the fitted model would generate. It is a model-checking tool rather than a hypothesis test, and it has no fixed 0.05 threshold. Values near 0 or 1 flag a model that cannot reproduce the chosen statistic; values in the middle mean the model reproduces it well.

How many times should I go around the loop?

Until the checks you care about pass and a sensitivity analysis shows your conclusions are stable. In practice that is often two or three iterations. Stop when new revisions stop changing the checks that matter, not when the model runs without error.

Do I always need rstanarm, or can I use brms or Stan?

Any engine works. This tutorial uses rstanarm because it fits standard regressions with one function call, but the same five steps apply to brms and to raw Stan. The workflow is a discipline, not a package.

The tutorial fit models with rstanarm but ran checks in base R. Why?

The model fitting needs a local R session with rstanarm installed. The checks are written in base R so you can run them here in your browser, drawing from a compact summary of the fitted posterior. For linear models like this one, that summary behaves just like the full posterior, which is why the browser checks and the local posterior_predict() version give the same verdict.

Summary

The Bayesian workflow turns a single model into a defensible conclusion by looping through five steps and revising whenever a check fails. On the mtcars question, the loop moved us from a model that wrongly declared transmission irrelevant to one that correctly captured a weight-by-transmission interaction, and it did so because a posterior predictive check forced the revision.

Figure 3: Each step and the tool it produces, at a glance.

The lesson that outlasts this example: a coefficient table can hide a model's failures, and only checking predictions reveals them. Set priors you can defend, always diagnose before interpreting, let the checks drive your revisions, and report the sensitivity of your conclusions. That loop is what makes a Bayesian analysis trustworthy.

Continue learning

• Prior Predictive Checks in R: go deeper on step 1, including how to elicit and calibrate priors before fitting.
• Posterior Predictive Checks in R: more test statistics and graphical checks for step 4.
• MCMC Diagnostics in R: the full set of convergence checks behind step 3, including trace plots and divergences.

References

1. Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C., and Modrák, M. Bayesian Workflow. arXiv:2011.01808 (2020). Link - the paper that named and organized the five-step loop this tutorial teaches.
2. McElreath, R. Statistical Rethinking: A Bayesian Course with Examples in R and Stan, 2nd Edition. CRC Press (2020). Link - a ground-up Bayesian course that builds the same modeling habits with worked R and Stan examples.
3. rstanarm documentation: How to Use the rstanarm Package. Link - the package reference for stan_glm() and the prior arguments used throughout.
4. Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A. Visualization in Bayesian workflow. Journal of the Royal Statistical Society A (2019). Link - shows the graphical prior and posterior predictive checks behind steps 1 and 4.
5. Vehtari, A., Gelman, A., and Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing (2017). Link - the method behind loo() and loo_compare() and how to read elpd_diff.
6. Muth, C., Oravecz, Z., and Gabry, J. User-friendly Bayesian regression modeling: A tutorial with rstanarm and shinystan. The Quantitative Methods for Psychology (2018). Link - a gentle applied walkthrough of Bayesian regression with rstanarm.
7. Stan Development Team. Posterior and Prior Predictive Checks (Stan User's Guide). Link - the reference explanation of prior and posterior predictive checks.