This page fits models by Bayesian inference, using the JAGS sampler driven from R: a single proportion, a logistic regression, a survival model, and a random-effects model, and then averages over linear regression models (Dobson and Barnett 2018, chap. 14). It uses the priors of the Bayesian Inference page and the sampling and convergence checks of the Markov Chain Monte Carlo page.
1 A first example: a single proportion
From here on, the models are fit by Bayesian inference with JAGS (“Just Another Gibbs Sampler”), a program that builds an MCMC sampler from a text description of a model (Dobson and Barnett 2018, chap. 13), driven from R through the rjags package. We begin with the simplest possible model, a single Bernoulli probability, whose exact posterior is known from the Beta-Bernoulli example, so the output can be checked. The example shows the mechanics that every later analysis reuses: specifying a model, supplying data, running a burn-in, monitoring parameters, and summarizing and checking the draws. In JAGS, dnorm(mean, precision) is parameterized by the precision, the reciprocal of the variance.
Example 1 (A Bernoulli model fitted with JAGS) The data are \(r = 55\) successes in \(n = 91\) trials, and the prior is uniform, as in the Beta-Bernoulli example. Each chain gets its own random-number generator seed, so the output is reproducible:
y <-rep(c(1L, 0L), times =c(55L, 36L))bernoulli_model <-"model { for (i in 1:N) { y[i] ~ dbern(p) } p ~ dbeta(1, 1)}"inits <-list(list(.RNG.name ="base::Mersenne-Twister", .RNG.seed =1),list(.RNG.name ="base::Mersenne-Twister", .RNG.seed =2))fit <- rjags::jags.model(textConnection(bernoulli_model),data =list(y = y, N =length(y)),inits = inits,n.chains =2,quiet =TRUE)stats::update(fit, n.iter =1000, progress.bar ="none") # burn-indraws <- rjags::coda.samples( fit,variable.names ="p",n.iter =4000,progress.bar ="none")p_draws <-as.matrix(draws)[, "p"]c(mean =mean(p_draws), stats::quantile(p_draws, c(0.025, 0.975)),Rhat = coda::gelman.diag(draws)$psrf[1, "Point est."]) |>round(3)#> mean 2.5% 97.5% Rhat #> 0.601 0.499 0.698 1.000
For a binary outcome \(Y_i \sim \operatorname{Bernoulli}(\pi_i)\) with \(\operatorname{logit}(\pi_i) = {\tilde{x}_i}^{\top}\tilde{\beta}\), the logistic regression model, a Bayesian analysis places a prior on the coefficient vector \(\tilde{\beta}\) and samples the posterior \(\operatorname{p}(\tilde{\beta}\mid \tilde{y})\) by MCMC (Dobson and Barnett 2018, chap. 14, p. 318). Each \(\beta_j\) is summarized by the mean and quantiles of its draws. No large-sample Gaussian approximation is needed: a credible interval is read directly from the posterior quantiles, and the posterior of an odds ratio \(e^{\beta_j}\), or of any other function of \(\tilde{\beta}\), is obtained by transforming the draws.
Example 2 (Logistic regression fitted by Bayesian inference) We simulate \(N = 200\) observations from a logistic model with one continuous predictor, intercept \(\beta_0 = -0.5\) and slope \(\beta_1 = 1.2\), and place a diffuse \(\operatorname{N}\mathopen{}\left(0, 10^2\right)\mathclose{}\) prior on each coefficient (precision \(0.01\) in JAGS). The JAGS program also monitors the odds ratio \(e^{\beta_1}\) for a one-unit increase in \(x\).
Each 95% credible interval contains its coefficient’s true value, and every \(\hat{R}\) is near 1. As comparing MCMC estimates to maximum likelihood leads us to expect, the posterior means are close to the maximum likelihood estimates:
stats::glm(y ~ x, family = stats::binomial()) |> stats::coef() |>round(3)#> (Intercept) x #> -0.647 1.183
3 Survival analysis
Parametric survival models, such as those with exponential or Weibull event times, can be fit by Bayesian inference: priors are placed on the baseline-hazard parameters and the regression coefficients, and the posterior is sampled by MCMC (Dobson and Barnett 2018, chap. 14, p. 327). Censoring enters through the likelihood, exactly as in maximum likelihood estimation with the likelihood with censoring: an observed event at time \(t\) contributes its density \({\lambda}(t)\,\operatorname{S}(t)\), and an observation right-censored at time \(t\) contributes its survival probability \(\operatorname{S}(t)\), where \({\lambda}\) is the hazard and \(\operatorname{S}\) the survival function.
When a model’s log-likelihood contribution \(\ell_i\) for observation \(i\) is not one of the distributions built into JAGS, the zeros trick supplies it. An observed value of 0 from a Poisson distribution with mean \(\phi_i\) has probability \(e^{-\phi_i}\), and with \(\phi_i \stackrel{\text{def}}{=}C - \ell_i\) for a constant \(C\),
\[
\begin{aligned}
e^{-\phi_i}
&= e^{-(C - \ell_i)}
&& \text{(definition of $\phi_i$)}\\
&= e^{-C} e^{\ell_i}
&& \text{(splitting the exponent)}\\
&\propto e^{\ell_i}
&& \text{($C$ does not depend on the parameters)},
\end{aligned}
\]
which is observation \(i\)’s likelihood contribution. So declaring data zeros[i] = 0 with zeros[i] ~ dpois(C - loglik[i]) multiplies the likelihood by \(e^{\ell_i}\). The constant \(C\) only has to be large enough to keep every \(\phi_i\) positive.
Example 3 (Exponential survival regression fitted by Bayesian inference) We simulate \(N = 200\) right-censored exponential survival times, with a binary covariate \(x_i\) (say, treatment) and hazard \({\lambda}_i = \operatorname{exp}\mathopen{}\left\{\beta_0 + \beta_1 x_i\right\}\mathclose{}\), where the log baseline hazard is \(\beta_0 = \log 0.05\) and the log hazard ratio is \(\beta_1 = -0.7\). With event indicator \(\delta_i\) (1 for an event, 0 for a censored time) and follow-up time \(t_i\), the exponential density is \({\lambda}_i e^{-{\lambda}_i t}\) and its survival function is \(e^{-{\lambda}_i t}\), so observation \(i\) contributes the log-likelihood
About 50% of the follow-up times end in an event. The 95% credible interval for \(\beta_1\) contains the true \(-0.7\), though the posterior mean is some distance from it. The maximum likelihood estimates, from the same log-likelihood, show that the gap is sampling variation in this data set rather than an effect of the prior:
The monitored \(e^{\beta_1}\) gives the hazard ratio and its credible interval directly.
4 Random effects
The hierarchical model of the two-level Gaussian example is the random-effects model structure, with group-level parameters \(\theta_j \sim \operatorname{N}\mathopen{}\left(\mu, \tau^2\right)\mathclose{}\) drawn from a common distribution. What changes here from a maximum likelihood fit is the inference method: the Bayesian analysis places priors on the hyperparameters \(\mu\) and \(\tau\) and samples their joint posterior together with the group-level parameters (Dobson and Barnett 2018, chap. 14, p. 329). The posterior shrinks each group’s estimate toward the overall mean, by an amount the data determine through \(\tau\), and the posterior for \(\tau\) carries the uncertainty about the between-group spread into every group-level summary.
Example 4 (A random-intercept model fitted by Bayesian inference) We simulate \(J = 8\) groups of 12 observations each, with group means drawn from \(\operatorname{N}\mathopen{}\left(\mu, \tau^2\right)\mathclose{}\) (\(\mu = 5\), \(\tau = 1.5\)) and within-group standard deviation \(\sigma = 2\). The priors are a diffuse \(\operatorname{N}\mathopen{}\left(0, 100^2\right)\mathclose{}\) for \(\mu\) and flat priors on \((0, 100)\) for the standard deviations \(\tau\) and \(\sigma\). Being flat on the standard-deviation scale is a choice: as a flat prior on the log-odds shows for a probability, it is not flat on another scale, such as the variance.
set.seed(2025)n_groups <-8Ln_per_group <-12Ltheta_true <- stats::rnorm(n_groups, mean =5, sd =1.5)group <-rep(seq_len(n_groups), each = n_per_group)y <- stats::rnorm(n_groups * n_per_group, mean = theta_true[group], sd =2)random_intercept_model <-"model { for (i in 1:N) { y[i] ~ dnorm(theta[group[i]], 1 / sigma^2) } for (j in 1:J) { theta[j] ~ dnorm(mu, 1 / tau^2) } mu ~ dnorm(0, 1.0E-4) sigma ~ dunif(0, 100) tau ~ dunif(0, 100)}"inits <-list(list(.RNG.name ="base::Mersenne-Twister", .RNG.seed =1),list(.RNG.name ="base::Mersenne-Twister", .RNG.seed =2))fit <- rjags::jags.model(textConnection(random_intercept_model),data =list(y = y, group = group, N =length(y), J = n_groups),inits = inits,n.chains =2,quiet =TRUE)stats::update(fit, n.iter =2000, progress.bar ="none")draws <- rjags::coda.samples( fit,variable.names =c("mu", "tau", "sigma", "theta"),n.iter =10000,progress.bar ="none")draws_matrix <-as.matrix(draws)posterior_summary <-cbind(mean =colMeans(draws_matrix),t(apply(draws_matrix, 2, stats::quantile, probs =c(0.025, 0.975))),Rhat = coda::gelman.diag(draws, multivariate =FALSE)$psrf[, "Point est."])round(posterior_summary[c("mu", "tau", "sigma"), ], 3)#> mean 2.5% 97.5% Rhat#> mu 5.572 4.733 6.405 1.002#> tau 0.868 0.125 2.002 1.004#> sigma 2.161 1.867 2.520 1.000
The 95% credible intervals contain the true \(\mu = 5\), \(\tau = 1.5\) and \(\sigma = 2\). The interval for \(\tau\) is wide: eight groups say little about how spread out group means are, and the posterior reports that uncertainty instead of a single estimate of the variance component.
Each group’s posterior mean lies between its sample mean and the overall sample mean:
This is the borrowing of strength described in distributions and hierarchies: group 4, with the largest sample mean, is pulled furthest toward the overall mean.
5 Bayesian model averaging
When several candidate models are plausible, committing to a single “best” one ignores the uncertainty about which model is right. Bayesian model averaging instead averages over the models, weighting each by its posterior probability (Dobson and Barnett 2018, chap. 14, p. 339). This approach carries model uncertainty, not just parameter uncertainty, into the final inference. It is an alternative to choosing a single model by predictor selection.
Definition 1 (Posterior model probability) Let \(M_1, \ldots, M_K\) be candidate models with prior probabilities \(\operatorname{p}(M_k)\) summing to 1, and let \(\operatorname{p}(\tilde{y}\mid M_k)\) be the marginal likelihood of the data under model \(M_k\). The posterior model probability of \(M_k\) is
Example 5 (Posterior probabilities of two models) Two models with equal prior probabilities \(\operatorname{p}(M_1) = \operatorname{p}(M_2) = 0.5\) and marginal likelihoods \(\operatorname{p}(\tilde{y}\mid M_1) = 0.02\) and \(\operatorname{p}(\tilde{y}\mid M_2) = 0.06\) have posterior probabilities \(\operatorname{p}(M_1 \mid \tilde{y}) = 0.01 / (0.01 + 0.03) = 0.25\) and \(\operatorname{p}(M_2 \mid \tilde{y}) = 0.75\).
Definition 2 (Bayesian model averaging) With the posterior model probabilities\(\operatorname{p}(M_k \mid \tilde{y})\) of candidate models \(M_1, \ldots, M_K\), Bayesian model averaging estimates a quantity \(\Delta\) that has the same meaning in every model by its posterior distribution averaged over the models:
Definition 3 (Posterior inclusion probability) When the candidate models of Definition 2 differ in which predictors they include, the posterior inclusion probability of a predictor is the total posterior probability of the models that include it.
Example 6 (A BIC approximation to Bayesian model averaging) Marginal likelihoods are hard to compute, but with equal prior probabilities for the models, the Bayesian information criterion gives a large-sample approximation to the posterior model probabilities (Schwarz 1978):
We simulate data in which only \(x_1\) and \(x_2\) affect the outcome, fit a linear regression for every subset of the three predictors \(x_1\), \(x_2\) and \(x_3\), and compute these approximate posterior model probabilities:
Subtracting the smallest BIC before exponentiating leaves the normalized weights unchanged and avoids numerical underflow. The posterior inclusion probability of each predictor sums the weights of the models that contain it:
includes <-function(v) vapply(subsets, function(s) v %in% s, logical(1))vapply(predictors, function(v) sum(weight[includes(v)]), numeric(1)) |>round(3)#> x1 x2 x3 #> 1.000 1.000 0.152
For the model-averaged slope of \(x_1\), we approximate its posterior mean within each model by the least-squares estimate, taken as \(0\) in models that exclude \(x_1\), and average with the posterior model probabilities as weights:
The predictors that affect the outcome, \(x_1\) and \(x_2\), have inclusion probabilities near 1, and the irrelevant \(x_3\) a much smaller one. The models without \(x_1\) get almost no weight, so the model-averaged slope of \(x_1\) is close to its least-squares estimate in the model with \(x_1\) and \(x_2\) only. Both differ from the true slope \(0.8\) by sampling variation in this data set, not because of the averaging.
References
Dobson, Annette J, and Adrian G Barnett. 2018. An Introduction to Generalized Linear Models. 4th ed. CRC press. https://doi.org/10.1201/9781315182780.