Bayesian Inference Foundations

Core theorem: prior × likelihood → posterior

Bayesian inference rests on a single mechanistic rule. Before seeing data, the analyst encodes existing knowledge — from prior studies, expert opinion, theoretical constraints, or deliberate agnosticism — as a probability distribution p(θ) over the unknown parameter θ. This is the prior. The data D inform via the likelihood L(θ; D) = p(D | θ): the probability of observing exactly these data if θ were the true value. Bayes' theorem combines both:

posterior p(θ | D) ∝ prior p(θ) × likelihood p(D | θ)

The ∝ symbol ("proportional to") hides only a normalizing constant — a number that makes the posterior integrate to 1. That constant does not depend on θ and does not affect inference about θ. The posterior is the complete statement of post-data knowledge: it assigns a probability to every interval of plausible θ values.

In a finite conjugate update (e.g., the beta-binomial pair described below), the posterior is computed in closed form with no numerical integration. In all other cases — logistic or survival regressions with Bayesian priors, hierarchical multi-site models, models with correlated parameters — the posterior is approximated by Markov chain Monte Carlo (MCMC).

Credible interval vs confidence interval — the interpretation most people want

The single most practically important distinction between Bayesian and frequentist inference is in the interval that summarizes uncertainty. A 95% credible interval (CrI) is the interval [a, b] such that the posterior probability that θ ∈ [a, b] is exactly 0.95. The CrI is a direct probability statement about the parameter: "Given the prior and the data, I am 95% certain the true event rate is between 0.12 and 0.45."

This is precisely the interpretation most scientists, clinicians, and decision-makers apply when they see a confidence interval — and it is wrong for a frequentist confidence interval. A 95% CI is a statement about a procedure, not about a specific interval: if the study were repeated infinitely many times and a CI constructed each time, 95% of those intervals would contain the true fixed value of θ.

For any single realized interval, θ either is or is not inside it — there is no 95% probability. The parameter θ is fixed and unknown; only the interval is random in the frequentist framework.

When sample sizes are large and priors are weak (non-informative), the numerical difference between a 95% CrI and a 95% CI is negligible. When samples are small, priors are informative, or parameters are on the boundary of their support, the two can diverge substantially. In rare-disease RWE settings with n < 50, the CrI is the only interval that correctly propagates all sources of uncertainty — prior, likelihood, and model structure — into a single coherent statement.

Conjugate priors and the beta-binomial example

A conjugate prior is a prior distribution family that, combined with a specific likelihood, yields a posterior in the same family — enabling exact closed-form updating. The canonical example for proportions and event rates is the beta-binomial pair:

Prior: Beta(α₀, β₀) [event rate p has prior mean α₀/(α₀+β₀)] Data: k events in n Bernoulli trials [likelihood: Binomial(n, p)] Posterior: Beta(α₀ + k, β₀ + n − k) [exact; no MCMC required]

The shape parameters α and β carry the interpretation of pseudo-counts: α₀ is the number of prior "events" and β₀ is the number of prior "non-events," giving the prior an effective sample size (ESS) of α₀ + β₀. The posterior simply adds observed counts to the pseudo-counts. The Jeffreys non-informative prior for a binomial proportion is Beta(0.5, 0.5), which contributes only 1 pseudo-observation.

The link to the beta-distribution entry in this catalog: the beta-distribution entry explains how to fit Beta(α, β) from a published mean and SE (method of moments), which is exactly how informative priors are constructed from published literature.

Prior choice: informative, weakly informative, skeptical, and flat

The prior encodes what is known before the data are seen. Four classes arise in practice:

• Informative prior: encodes genuine external knowledge — a Beta(2, 8) prior on an event rate says "I believe the rate is around 0.20 and I have 10 prior pseudo-observations of evidence." Legitimate when external evidence is well-documented, transportable to the current patient population, and pre-specified before analysis.
• Weakly informative prior: rules out physically impossible or substantively absurd values while remaining flat over the scientifically plausible range. A Beta(1.5, 1.5) or a Normal(0, 1) prior on a log-odds are examples: they prevent the MCMC sampler from exploring extreme regions without pulling the posterior toward any particular interior value. Weakly informative priors are the default recommendation for hierarchical models in RWE.
• Skeptical prior: deliberately weights against large treatment effects, raising the evidentiary bar. A Normal(0, 0.355²) on log-OR places 95% of prior mass on OR between 0.5 and 2.0. Regulators and HTA reviewers increasingly require a prespecified skeptical prior to accompany the primary analysis in confirmatory Bayesian submissions.
• Flat / non-informative prior: places approximately equal mass across the parameter range. Truly flat priors are mathematically improper on unbounded parameters and are not invariant to reparameterization. In large samples a flat prior yields posteriors nearly identical to frequentist maximum-likelihood results — making it a natural bridge when the goal is Bayesian machinery with minimal prior impact.

Prior sensitivity analysis is mandatory reporting. Every Bayesian submission should include results under at least two or three prior specifications spanning from non-informative to the analyst's preferred informative prior. Robustness of the conclusion across this grid means the prior is not load-bearing; strong sensitivity means the data are insufficient to overcome prior beliefs — which is itself a substantive finding that must be disclosed. Regulatory agencies (FDA 2023 Bayesian guidance) and HTA bodies (NICE DSU Technical Support Documents) require prior sensitivity analyses for Bayesian submissions.

MCMC and convergence diagnostics

Most practical Bayesian models do not have conjugate closed forms. Markov chain Monte Carlo (MCMC) — Metropolis-Hastings, Gibbs sampling, and the modern No-U-Turn Sampler (NUTS) used in Stan, PyMC, and brms — approximates the posterior by constructing a Markov chain whose stationary distribution is the target posterior. After a burn-in (warmup) phase where the chain finds the high-probability region and tunes its step size, the retained draws are treated as samples from the posterior.

Convergence diagnostics exist not to prove convergence (that cannot be proven) but to detect non-convergence: the R-hat (potential scale reduction factor) should be < 1.01 for every parameter; Effective Sample Size (ESS) should exceed 400 in the bulk and tail. Trace plots should look like well-mixed "hairy caterpillars" — no trends, no stuck chains, no chains drifting apart. Always report these diagnostics. A posterior from a non-converged chain is meaningless regardless of how sophisticated the model. Use Stan/brms (R) or PyMC (Python) — they implement NUTS and report diagnostics automatically.

Where Bayes earns its keep in RWE/HEOR

Bayesian inference is not universally superior to frequentist inference; it earns its additional modeling and computational cost in specific settings:

• Rare-disease historical borrowing: the beta-binomial conjugate update (or a MAP prior) lets external control data informatively augment a small concurrent trial with a quantifiable, pre-specified discount. See the Bayesian borrowing entry in this catalog.
• Probabilistic sensitivity analysis in cost-effectiveness: PSA is already implicitly Bayesian — it propagates parameter uncertainty through a decision model by drawing from distributions. Embedding the PSA in a formal Bayesian framework (parameters as posterior draws from a hierarchical model) gives coherent posterior probabilities over cost-effectiveness thresholds.
• Sequential safety monitoring: Bayesian adaptive designs accumulate posterior probability of harm continuously and can trigger early stopping when P(harm exceeds threshold | data) crosses a pre-specified gate, without paying a frequentist multiplicity penalty for interim looks.
• Hierarchical shrinkage across sites: a Bayesian hierarchical model partially pools site-specific estimates toward the overall mean, shrinking noisy small-site estimates and reducing false-discovery rates vs independently estimated site effects — mathematically equivalent to empirical Bayes / random-effects meta-analysis but with explicit priors on the between-site variance.
• Bias priors in probabilistic bias analysis: unmeasured confounding can be modeled by placing a prior over bias parameters (sensitivity and specificity of confounder measurement, prevalence in exposed vs unexposed), propagating the resulting uncertainty into adjusted effect estimates. See the probabilistic bias analysis entry in this catalog.

Frequentist-Bayes pragmatism (calibrated Bayes)

The tribal framing — "Bayesian vs frequentist" — is not useful in applied RWE/HEOR. Most production analyses use both paradigms at different layers. A Cox proportional hazards model is frequentist; adding Normal(0, σ²) priors on the regression coefficients and sampling via MCMC makes it Bayesian. Calibrated Bayes (Rubin 1984) is the pragmatic synthesis: use Bayesian machinery where it offers genuine advantages (informative priors, hierarchical structure, exact small-sample uncertainty, direct probability statements), and frequentist methods where they are adequate (large samples, no genuine prior, regulatory conventions requiring frequentist operating characteristics).

The key sanity check: a Bayesian posterior interval should have approximately correct frequentist coverage when the model is well specified. When Bayesian and frequentist results agree, the prior is not carrying load. When they disagree, the prior is doing real work — and must be justified explicitly.

Pros, cons, and trade-offs

Pros of Bayesian inference:

• Direct probability statements about parameters: "There is a 94% posterior probability the drug reduces 12-month mortality by more than 5 percentage points."
• Natural, principled incorporation of prior information in small-sample and rare-disease settings — the prior ESS is a transparent, reportable quantity.
• Coherent uncertainty propagation through complex hierarchical models without ad hoc variance adjustments.
• No frequentist multiplicity penalty for pre-specified interim looks in adaptive designs.
• The 95% credible interval has the interpretation practitioners routinely and incorrectly apply to confidence intervals.

Cons and limitations:

• Prior specification adds a modeling choice that must be pre-specified, justified, and sensitivity-analyzed; a poorly chosen or post-hoc prior can bias inference in opaque ways.
• Computational cost: MCMC can be slow for large RWE datasets or models with many parameters; convergence must be verified and can fail.
• Frequentist operating characteristics (type-I error, power) must still be verified by simulation for regulatory submissions — the Bayesian framework does not guarantee error control unless the design is pre-specified.
• Regulatory reviewers trained in frequentist methods may require additional justification; acceptance is increasing but uneven across jurisdictions.

When NOT to use — and when it is actively misleading

• Do not choose the prior to reach a desired result. Prior specification must be locked in the analysis plan before unblinding. Post-hoc prior tuning to achieve statistical significance is data dredging with additional obfuscation. Pre-specify the prior ESS and sensitivity grid.
• Do not skip prior sensitivity analysis. Reporting only the result under one prior — even a nominally "non-informative" flat prior — hides sensitivity. Every Bayesian submission should include at least three prior specifications from non-informative to the planned informative prior.
• Do not interpret a high posterior probability as a causal claim. A 98% posterior probability that a treatment coefficient is positive is an associational statement about the model, not proof of causation. Causal claims require the same identification assumptions regardless of inferential framework.
• Do not use flat priors on variance components in hierarchical models. A Uniform(0, ∞) prior on a between-site standard deviation is improper and causes pathological posterior behavior with few groups. Use a weakly informative Half-Normal or Half-Cauchy prior.
• Bayesian machinery does not fix a confounded design. Sophisticated MCMC sampling of a misspecified or confounded model yields a precise posterior of a wrong quantity.

Interpreting the output

From the worked example: prior Beta(2,8) (mean 0.2, equivalent to 10 pseudo-observations); observe 6 events in 20 trials; posterior Beta(8,22), mean 8/30 ≈ 0.267.

(1) Formal interpretation. The posterior Beta(8,22) is a probability distribution over the unknown event rate p. The posterior mean 8/30 ≈ 0.267 is a weighted compromise between the prior mean (2/10 = 0.2) and the data MLE (6/20 = 0.3), with weights proportional to the prior information (10 pseudo-observations) and the data (20 observations). The total information in the posterior is 8+22 = 30 pseudo-observations.

A 95% credible interval from Beta(8,22) spans approximately [0.12, 0.45]: the posterior probability that the true event rate p lies in this interval is exactly 0.95 — a direct probability statement, not a repeated-sampling claim. This posterior is exact because of conjugacy; no MCMC is required for this model.

(2) Practical interpretation. Before seeing the new data, we believed the event rate was around 20%, with the strength of 10 observations. We then observed 30% in 20 patients. Our best updated estimate is approximately 26.7% — not 30%, because the prior pulled us back toward 20%; not 20%, because the new data pulled us upward.

A decision-maker should note that the prior is doing real work: if its 10 pseudo-observations came from a comparable population, the shrinkage is appropriate. If the prior population differed (different disease stage, earlier era), the posterior is biased toward the prior, and a sensitivity analysis under Beta(0.5, 0.5) — contributing only 1 pseudo-observation — should accompany the primary result to show how much the conclusion depends on the informative prior.

Decision diagram

flowchart TD
  Prior[Prior p(theta)\nencodes existing knowledge\ne.g. Beta_2_8 mean=0.20] --> Bayes[Bayes' theorem\nposterior prop prior x likelihood]
  Data[Observed data\ne.g. 6 events in 20 trials\ndata MLE = 0.30] --> Bayes
  Bayes --> Post[Posterior p(theta|data)\nBeta_8_22 mean approx 0.267]
  Post --> CrI[95% credible interval\ndirect probability statement\nP(theta in CrI) = 0.95]
  Post --> Prob[Direct probability queries\nP(rate > 0.25 | data) = ?]
  Post --> Sens[Prior sensitivity analysis\nrepeat under flat and skeptical priors\nmandatory reporting]
  Sens --> Report[Report: posterior is robust or\nprior is load-bearing - say which]

Worked example

Scenario

A health outcomes analyst is studying an adverse event rate for a drug used in a rare pediatric disease. From published literature, the event rate is believed to be around 20%, based on evidence roughly equivalent to 10 patients. The analyst encodes this as a prior Beta(2,8) — 2 prior events and 8 prior non-events. A new registry study then observes 6 adverse events among 20 patients. The analyst wants to compute the updated (posterior) event rate using the exact beta-binomial conjugate formula, and confirm it lies between the prior belief (20%) and the data estimate (30%).

Dataset

Summary counts feeding the beta-binomial conjugate update. Each row represents one source of information: prior belief (encoded as pseudo-counts) and new observed data.

Steps

Result

Posterior Beta(8,22): prior mean 2/10 = 0.2, data MLE 6/20 = 0.3, posterior mean 8/30 ≈ 0.267. The posterior ESS = 8+22 = 30, composed of 10 prior pseudo-observations plus 20 real data observations. The 95% credible interval from Beta(8,22) is approximately [0.12, 0.45] — the posterior probability the true event rate lies in this range is exactly 0.95. Because the posterior mean (0.267) falls between the prior (0.2) and the data (0.3), it is correctly shrunk toward the prior in proportion to how informative the prior was relative to the new data.

Trade-offs

Runnable example

import numpy as np
from scipy import stats

# ── Beta-binomial conjugate update (exact posterior, no MCMC) ─────────────────
alpha_0, beta_0 = 2, 8          # prior Beta(2,8): mean = 2/(2+8) = 0.2; ESS = 10
k, n            = 6, 20         # observe 6 events in 20 trials

alpha_post = alpha_0 + k        # 2 + 6 = 8
beta_post  = beta_0 + (n - k)  # 8 + (20 - 6) = 22

prior     = stats.beta(alpha_0, beta_0)
posterior = stats.beta(alpha_post, beta_post)

prior_mean = alpha_0 / (alpha_0 + beta_0)             # 2/10 = 0.200  (exact)
data_mle   = k / n                                     # 6/20 = 0.300  (exact)
post_mean  = alpha_post / (alpha_post + beta_post)     # 8/30 ≈ 0.267

print(f"Prior     Beta({alpha_0}, {beta_0}):  mean = {prior_mean:.4f}")
print(f"Data MLE: {k}/{n} = {data_mle:.4f}")
print(f"Posterior Beta({alpha_post}, {beta_post}): mean = {post_mean:.4f}")

# 95% credible interval — direct posterior probability statement
cri_lo, cri_hi = posterior.ppf(0.025), posterior.ppf(0.975)
print(f"95% credible interval: [{cri_lo:.3f}, {cri_hi:.3f}]")
print("Interpretation: P(event rate in [{:.3f}, {:.3f}] | data) = 0.95".format(cri_lo, cri_hi))

# Direct posterior probability query: P(p > 0.25 | data)
prob_q = 1 - posterior.cdf(0.25)
print(f"P(event rate > 0.25 | prior + data): {prob_q:.3f}")

# Prior sensitivity: Jeffreys Beta(0.5, 0.5) (non-informative)
a_jeff = 0.5 + k
b_jeff = 0.5 + (n - k)
post_jeff = stats.beta(a_jeff, b_jeff)
print(f"\nSensitivity — Jeffreys prior Beta(0.5,0.5):")
print(f"  Posterior Beta({a_jeff}, {b_jeff}): mean = {a_jeff/(a_jeff+b_jeff):.4f}")
print(f"  95% CrI: [{post_jeff.ppf(0.025):.3f}, {post_jeff.ppf(0.975):.3f}]")

# ── PyMC sketch: Bayesian logistic regression (non-conjugate) ─────────────────
# import pymc as pm
# import arviz as az
# with pm.Model():
#     # Weakly informative priors on log-odds scale
#     intercept = pm.Normal("intercept", mu=0, sigma=2)
#     beta_treat = pm.Normal("beta_treat", mu=0, sigma=1)   # skeptical prior
#     beta_age   = pm.Normal("beta_age",   mu=0, sigma=1)
#     p = pm.math.invlogit(intercept + beta_treat * treat + beta_age * age)
#     obs = pm.Bernoulli("obs", p=p, observed=y)
#     idata = pm.sample(2000, tune=2000, target_accept=0.9,
#                       chains=4, random_seed=1)
# # Convergence diagnostics — MANDATORY before interpreting posterior
# assert az.summary(idata)["r_hat"].max() < 1.01, "R-hat >= 1.01: chain not converged"
# print(az.summary(idata, var_names=["beta_treat"], hdi_prob=0.95))
# # Direct posterior probability: P(treatment increases rate)
# print("P(beta_treat > 0 | data):", float((idata.posterior["beta_treat"] > 0).mean()))

# ── Beta-binomial conjugate update (exact posterior, base R) ──────────────────
alpha_0 <- 2;  beta_0 <- 8    # prior Beta(2,8): mean = 2/(2+8) = 0.2; ESS = 10
k <- 6;  n <- 20              # observe 6 events in 20 trials

alpha_post <- alpha_0 + k              # 2 + 6 = 8
beta_post  <- beta_0 + (n - k)        # 8 + 14 = 22

prior_mean <- alpha_0 / (alpha_0 + beta_0)            # 2/10 = 0.200
data_mle   <- k / n                                    # 6/20 = 0.300
post_mean  <- alpha_post / (alpha_post + beta_post)   # 8/30 ≈ 0.267

cat(sprintf("Prior     Beta(%d,%d): mean = %.4f\n", alpha_0, beta_0, prior_mean))
cat(sprintf("Data MLE: %d/%d = %.4f\n", k, n, data_mle))
cat(sprintf("Posterior Beta(%d,%d): mean = %.4f\n", alpha_post, beta_post, post_mean))

# 95% credible interval — P(theta in CrI | data) = 0.95 exactly
cri <- qbeta(c(0.025, 0.975), alpha_post, beta_post)
cat(sprintf("95%% credible interval: [%.3f, %.3f]\n", cri[1], cri[2]))

# Direct posterior probability: P(p > 0.25 | data)
cat(sprintf("P(event rate > 0.25 | data): %.3f\n",
            1 - pbeta(0.25, alpha_post, beta_post)))

# Prior sensitivity: Jeffreys Beta(0.5, 0.5) — non-informative
a_jeff <- 0.5 + k;  b_jeff <- 0.5 + (n - k)
cat(sprintf("\nSensitivity check — Jeffreys prior Beta(0.5, 0.5):\n"))
cat(sprintf("  Posterior Beta(%.1f, %.1f): mean = %.4f\n",
            a_jeff, b_jeff, a_jeff / (a_jeff + b_jeff)))
cat(sprintf("  95%% CrI: [%.3f, %.3f]\n",
            qbeta(0.025, a_jeff, b_jeff), qbeta(0.975, a_jeff, b_jeff)))

# ── brms sketch: Bayesian logistic regression via Stan ────────────────────────
# library(brms)
#
# fit <- brm(
#   outcome ~ treatment + age + cci,
#   data   = rwe_df,
#   family = bernoulli(link = "logit"),
#   prior  = c(
#     prior(normal(0, 2), class = Intercept),
#     prior(normal(0, 1), class = b, coef = treatment),  # skeptical prior
#     prior(normal(0, 1), class = b)                      # weakly informative
#   ),
#   chains = 4, iter = 2000, warmup = 1000, seed = 1,
#   backend = "cmdstanr"   # faster; requires cmdstanr + CmdStan
# )
#
# # MANDATORY convergence check before interpreting results
# stopifnot(max(rhat(fit)) < 1.01)
# cat("All R-hat < 1.01: chains have converged.\n")
# summary(fit)  # posterior means, 95% CrI, Bulk ESS, Tail ESS for each parameter
#
# # Direct posterior probability: P(treatment coefficient > 0 | data)
# h <- hypothesis(fit, "treatment > 0")
# print(h)   # P+ = posterior probability treatment increases log-odds

/* ── Beta-binomial conjugate update (exact, DATA step) ─────────────────────── */
data work.bayes_update;
  /* Prior Beta(2,8): prior mean = 2/(2+8) = 0.200, ESS = 10 pseudo-observations */
  alpha_0   = 2;   beta_0    = 8;
  k         = 6;   n         = 20;

  /* Conjugate posterior                                                          */
  alpha_post = alpha_0 + k;             /* 2 + 6 = 8                             */
  beta_post  = beta_0  + (n - k);      /* 8 + 14 = 22                            */
  post_total = alpha_post + beta_post;  /* 8 + 22 = 30                           */

  prior_mean = alpha_0 / (alpha_0 + beta_0);        /* 2/10 = 0.200              */
  data_mle   = k / n;                                /* 6/20 = 0.300              */
  post_mean  = alpha_post / post_total;              /* 8/30 approx 0.267         */

  /* 95% credible interval from Beta(8,22)                                        */
  cri_lo = quantile('BETA', 0.025, alpha_post, beta_post);
  cri_hi = quantile('BETA', 0.975, alpha_post, beta_post);

  /* P(event rate > 0.25 | prior + data)                                          */
  prob_gt_025 = 1 - cdf('BETA', 0.25, alpha_post, beta_post);

  /* Sensitivity: Jeffreys Beta(0.5, 0.5) non-informative prior                   */
  a_jeff = 0.5 + k;   b_jeff = 0.5 + (n - k);
  post_mean_jeff = a_jeff / (a_jeff + b_jeff);
run;

proc print data=work.bayes_update label noobs;
  var prior_mean data_mle post_mean cri_lo cri_hi prob_gt_025 post_mean_jeff;
  format prior_mean data_mle post_mean cri_lo cri_hi prob_gt_025 post_mean_jeff 8.4;
  label prior_mean     = "Prior mean (Beta 2,8)"
        data_mle       = "Data MLE"
        post_mean      = "Posterior mean"
        cri_lo         = "95% CrI lower"
        cri_hi         = "95% CrI upper"
        prob_gt_025    = "P(rate>0.25|data)"
        post_mean_jeff = "Posterior mean (Jeffreys prior)";
  title "Beta-binomial conjugate update: prior Beta(2,8) + 6 events in 20";
run;

/* ── PROC MCMC: Bayesian logistic regression for binary RWE outcome ─────────── */
/* Requires work.rwe_data with variables: outcome (0/1), treat, age, cci.        */
/*
proc mcmc data=work.rwe_data nbi=2000 nmc=10000 seed=1 thin=2
          monitor=(b_int b_treat b_age b_cci)
          diagnostics=all;   * prints R-hat (PSRF) and ESS — check PSRF < 1.01;
  parms b_int 0  b_treat 0  b_age 0  b_cci 0;
  * Weakly informative priors (normal on log-odds scale);
  prior b_int   ~ normal(0, var=4);   * Normal(0, 2^2);
  prior b_treat ~ normal(0, var=1);   * Normal(0, 1^2) skeptical;
  prior b_age   ~ normal(0, var=1);
  prior b_cci   ~ normal(0, var=1);
  p_event = logistic(b_int + b_treat*treat + b_age*age + b_cci*cci);
  model outcome ~ binomial(n=1, p=p_event);
  * After run: verify PSRF < 1.01 and ESS > 400 for all parameters;
run;
*/

/* ── PROC GENMOD BAYES: faster alternative for standard GLMs ────────────────── */
/* The BAYES statement adds MCMC sampling to a standard GENMOD call.             */
/*
proc genmod data=work.rwe_data;
  class treat;
  model outcome = treat age cci / dist=binomial link=logit;
  bayes seed=1 nmc=10000
        coeffprior=normal(mean=0 var=1)
        diagnostics=all   * check convergence before trusting posterior summaries;
        plots=none;
run;
*/

Citations

1. [1]Bland JM, Altman DG. Statistics notes: Bayesians and frequentists. BMJ. 1998;317(7166):1151-1160.
2. [2]Greenland S. Bayesian perspectives for epidemiological research: I. Foundations and basic methods. International Journal of Epidemiology. 2006;35(3):765-775.

1. [3]Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Statistical Science. 1992;7(4):457-472.