Components of a Bayesian Analysis

In contrast to the classical statistical approach (a.k.a., frequentist) which conducts inferences using only the existing data, Bayesian analysis combines information in the data, through the likelihood, with prior information about the model parameters to obtain a combined assessment of uncertainty of these unknown quantities called the posterior (Table 1). Our examples are regression examples, and the parameters are the regression coefficients and the error variance for linear regression. Through an application of the well-known Bayes theorem, these components are linked as follows: P(parameters|data) ∝L (data|parameters) X P (parameters). The likelihood, L (data|parameters), and prior distribution of the parameters, P(parameters), are defined by the user. Computing the posterior distribution of the parameters, P (parameters|data), or distribution of the parameters given the data, is the analysis goal and analogous to estimating the regression coefficients and computing confidence intervals in frequentist analysis.

Table 1. Glossary of terms

Likelihood: The likelihood is related to the model, or distribution, of the data (i.e., outcome) in the study, which is the case in both classical and Bayesian approaches. In linear regression, the standard model of the data (or outcome) for participant i is assumed to be a normal distribution with a mean equal to β ₀ + β ₁ X _1i + ⋯ + β _k X _ki , where the X _ji ’s are the variables of interest measured on participant i, and a variance of σ _e ². Given the independence assumption, the model for all participants is the product of these normal distributions. The likelihood function is the same mathematical function except now the data are fixed at the observed data and the function is studied as a function of the parameters. Thus, the likelihood function characterizes how “likely” the parameters are for a given set of data. In logistic regression, the likelihood is based on a binomial distribution where a logit transformation of the probability of success (p) is linked to variables of interest as follows: $\log \left({p \over 1-p}\right)=\beta _{0}+\beta _{1}X_{1i}+\cdots +\beta _{k}X_{ki}$ . By selecting a linear or logistic regression analysis, the likelihood is typically specified by the package/procedure. The analyst selects the variables of interest when specifying the model in the statistical code.

Prior Distribution on Parameters: Prior distributions on the parameters quantify a belief in the parameter values before observing any study data, with the spread of the distribution quantifying the strength of those beliefs. In linear and logistic regression, the main parameters are the beta coefficients and typical distributions for the prior are normal distributions. Normal distributions are chosen for mathematical convenience but also allow for a full range of the parameters (−∞,+∞) to be accommodated in the prior. In linear regression, there is a prior distribution chosen for the inverse of the model error variance (model precision), which is typically a gamma distribution to take on only non-negative values. Other choices include truncated normal, uniform, half T, or half Cauchy. There is vast literature on defining prior distributions and eliciting prior distributions that capture a scientist’s beliefs about parameter values [Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin19–Reference Johnson, Tomlinson, Hawker, Granton and Feldman21]. In practice, we recommend C&T teams discuss different forms of prior knowledge, including a clinician’s perspective from past research or clinical care and information published from other studies. We also recommend analysis with several different priors identified with this approach to formally understand how conclusions may change based on a range of a priori assumptions C&T researchers may have about parameters. In this work, we specifically investigate various types, or “strengths,” of priors (Fig. 1). These ranged from vague with a wide, flatter distribution to skeptically and optimistically informative depending on the value of the mean in the priors for the betas. Priors with means set equal to 0 with more prior probability of no treatment difference (created by a smaller or tighter variance) compared to a vague prior have been labeled skeptical and quantify the strength of a clinician’s prior belief that the treatment will not have an effect. Priors with means not equal to 0 (negative or positive depend on clinical context) with more prior probability of a treatment difference (again, created with a smaller variance) compared to a vague prior have been labeled optimistic toward associations and quantify the strength of a clinician’s prior belief that the treatment will have an effect. As discussed earlier, the value of the mean can be chosen based on clinician input, historical data (e.g., previous trials), or biological plausibility. The labeling of priors as vague, skeptical, and optimistic depends directly on the scale of the outcome. For example, in one problem a variance of 1000-units² could be quite large and vague while another outcome with large values and highly variable measurements may need a much larger variance to be vague.

Figure 1. Panel A). The four priors on the treatment effect in the linear regression. The black solid line is the vague prior and shows an even weight for the largest range of the values. The gray dashed line is the “pseudo” vague prior, which has more weight around zero over a fairly large range of treatment effect values (but does not cover the range of values consistent with distribution of the outcome and thus, informative for the intercept). The red medium dashed line is the skeptical prior with much more weight centered around small treatment effects and the blue dotted line is the optimistic prior with nearly all its weight on treatment effects less than zero. Panels B and C). The linear regression coefficient values of the treatment group of the prior (green solid line), likelihood (solid blue), and posterior (orange dashed). Panel B plots values from the vague prior scenario showing that this prior specification does not pull the coefficient from the likelihood, as the two density curves are nearly identical. Panel C plots the informative prior scenario showing that an informative prior can influence the posterior from the likelihood, and the posterior is a combination of the prior and likelihood curves.

Posterior distribution: The posterior distribution is the distribution of the parameters conditioned on the data and expresses the uncertainty in the parameters after observing the data. As noted earlier, the posterior is proportional to the likelihood times the prior emphasizing that, on a log scale, the information about the parameters in the posterior is essentially the sum of information in the prior and the data. It is common to summarize the posterior by its mean, median, or mode and credible intervals. For example, a 95% equal-tailed probability or credible interval is computed by finding the 2.5^th and 97.5^th percentiles of the posterior distribution. This interval is interpreted as the 95% posterior probability that the parameter is between the bounds given the observed data. We note that this interpretation is commonly misapplied to the classical confidence interval (i.e., if we repeated sampling infinitely, the resulting 95% confidence intervals would contain the true population value 95% of the time; therefore, the correct interpretation is we are 95% confident that the bounds contain the parameter value) [Reference Hoekstra, Morey, Rouder and Wagenmakers22–Reference Andrade25]. Another commonly used credible interval is the highest posterior density (HPD) interval, which is the interval with the smallest interval width among all credible intervals. Using the posterior distribution, one can also compute posterior probabilities of ranges of parameter values aligning with clinically meaningful differences between the treatment arms. An example might be the posterior probability that the treatment has a higher (or lower) mean compared to the control group.

Computing the Posterior

In Bayesian regression analysis, the posterior is rarely a known distribution. Instead, Markov Chain Monte Carlo (MCMC) algorithms are used to simulate samples from the posterior. MCMC is a set of algorithms for sampling from a distribution even when the actual distribution cannot be mathematically derived [Reference Hastings26–Reference Martin, Frazier and Robert28]. Gibbs sampling and Metropolis-Hastings (M-H) are two common algorithms among many other approaches [Reference Geman and Geman29–Reference Green, Łatuszyński, Pereyra and Robert31]. Detailed descriptions are available elsewhere [Reference Martin, Frazier and Robert28]. The beauty of software development in the past decade is that, in SAS, R, or STATA, the user only needs to specify the variables in the regression model, similar to classical analysis, and distribution of the priors. The programs assemble the math for the MCMC algorithms behind the scenes for the user making Bayesian implementation highly feasible for most analysts. Code for each package and software are publicly available at https://github.com/nichole-carlson/BayesianClinicalResearch.

The output of MCMC algorithms is a rectangular dataset where columns are a sample from the posterior distribution of a model parameter and each row is a single iteration from the algorithm. The set of samples is often called “chains” reflecting that the samples are derived from a Markov chain, the mathematical system behind MCMC. By design, chain iterations are serially correlated (called autocorrelation). Assessing the autocorrelation strength is useful for investigating algorithm performance including likely convergence to the posterior and full posterior exploration. In practice, some run the MCMC algorithm 2–3 times, each a unique chain of sampling with different random starting values to further assess convergence. However, once adequate algorithm performance is determined, a single longer chain is used for final posterior estimation [Reference Geyer32].

To investigate these concepts, there are several graphical diagnostics to be assessed for each column. The most common diagnostic plots are trace, autocorrelation, and density plots (Fig. 2). A trace plot has parameter values on the y-axis for each iteration (x-axis). The autocorrelation plot is the pairwise correlation (y-axis) between parameter values for different lags between iterations (x-axis). When multiple chains are run, we first assess whether, after a suitable number of iterations, the trace plots start to overlap. Lack of overlap raises the possibility that convergence has not been achieved in the number of iterations selected. Assuming the trace plots eventually converge, further graphical assessment (and posterior summary measures) should not be visually influenced by a particular starting value. Thus, a portion of the initial algorithm iterations (called burn-in) is removed. The burn-in is often ∼ 10% of the iterations with no firm recommendation. The trace plot looks like random noise when convergence is likely (Fig. 2, top-left panel) versus having strong patterns or long stretches of wandering (Fig. 2, top-right panel). The autocorrelation plot exhibits a steep decline as the lags increase (Fig. 2, bottom-left, both panels). Strong patterns in the trace plots and long lags with a high correlation indicate poor mixing and a lower possibility of convergence or limited exploration of the posterior. The results may not be valid. Density (or histogram) plots of the parameters should have a smooth shape with a single mode (Fig. 2, bottom-right, both panels). In typical regression analyses, other bimodal or unusual patterns may also indicate lack of convergence. There are other diagnostics, e.g., the Gelman-Rubin statistic, for assessing convergence [Reference Gelman and Rubin33,Reference Vats and Knudson34]. For simplicity, we focused on visual diagnostics.

Figure 2. Diagnostic plots for various scenarios. The left panel indicates convergence is likely and the right where convergence is less likely and the MCMC algorithm is modified. The top figure in each panel is a trace plot. The bottom-left figure is an autocorrelation plot, and the bottom-right figure is a posterior density plot. These were generated by SAS PROC MCMC. Similar graphics are available for the other software.

After verifying the simulation has likely converged to the posterior distribution, the posteriors are summarized using the mean and measures of variability described by credible intervals and/or posterior standard errors and standard deviations.

Application of Bayesian Analysis to a Clinical Trial

Here we draft a formal methods section for our illustrative application as an example for the reader to follow when writing their own sections in clinical journals. Key components to consider are in Table 2. Our anesthesiologist clinical research collaborators were interested in evaluating how sublingual sufentanil, a novel opioid medication for moderate to severe pain, performed relative to the existing standard of care therapy of intravenous (IV) fentanyl. The study details are published elsewhere [Reference Berg, Habeck, Strigenz, Pearson, Kaizer and Hutchin35]. In brief, 75 patients were randomized to two study arms, and 66 were included in the per-protocol illustrative analysis. The exposure of interest was drug treatment with either sublingual sufentanil or IV fentanyl (the referent group). The primary continuous outcome was time to readiness for discharge after arrival in post-anesthesia care unit (PACU; in total minutes), and the primary dichotomous outcome was if a preoperative nerve block was administered (yes or no; probability of yes was modeled). All modeling was performed in parallel with R v4.2.1 (Vienna, Austria), SAS version 9.4 (Cary, NC, USA), and STATA version 17.0 (College Station, TX, USA).

Table 2. Statistical components to include in a Bayesian data analysis plan

Analysis of Time to Readiness for Discharge

Linear regression was used to model the association between treatment and discharge with and without adjustment for covariates. In the unadjusted models, the intercept represented the mean time to readiness for discharge in the IV fentanyl treatment group and the treatment coefficient represented the differences in the mean time to readiness for discharge from PACU between the sufentanil and IV fentanyl treatment groups. The adjusted models included centered procedure length (minutes) and sex (dichotomous).

The primary prior chosen for our analysis was a vague N(mean = 0, variance = 10,000) prior giving almost equal a priori weight to a large range of plausible values (Table 3; Fig. 1 panel A). In this prior distribution, there was an 80% a priori probability that the treatment difference was between −128 and 128 minutes (the values used to compute this range are the 10^th and 90^th percentiles of the prior distribution) and 50% of the a priori treatment effect values were less than 0. This same prior was also chosen for the regression coefficients on procedure length and sex in the adjusted model. The vague priors were considered conservative and similar to traditional analytic approaches.

Table 3. Specified priors

IG = inverse gamma (shape, scale); N = normal (mean, variance).

We also considered informative optimistic and skeptical priors (Table 3; Fig. 1 panel A) reflecting two common prior beliefs held by the clinical investigators. The optimistic prior was a N (−30, 100) reflecting a greater a priori belief in a clinically meaningful treatment difference. In this setting, the treatment difference (−30 minutes) and variance were selected based on the assumptions of the a priori power calculation, which relied upon historical data and clinician input. The optimistic prior had an 80% a priori probability of a treatment difference between −42.8 and −17.2 minutes with 99.8% of the a priori treatment differences less than 0. The skeptical prior was a N (0, 100) reflecting a greater a priori belief in no treatment effect with the variation about zero chosen from the a priori power calculation. The skeptical prior had an 80% a priori probability of a treatment difference between −12.8 and 12.8 minutes. In both scenarios, vague priors of N (0, 10000) were chosen for the regression coefficients on the covariates and intercept in the adjusted model.

We also considered a N (0, 1000) prior for all the parameters in the regression model (not just the treatment coefficient), representing the same numerical values but specified on different units. We classified this as a “pseudo” vague prior in that it was intended to be vague for all parameters, but further investigation showed it was informative for the intercept term. We selected this prior because it is a common choice among early adopters of Bayesian analyses, often thought to be a vague prior without careful consideration of the units of the outcome. Upon visual inspection of the range of parameter values allowed for the intercept compared to the distribution of the outcome, this prior was not vague and quite informative. This represents a case scenario that does not perform as intended with unintended high bias in the treatment coefficient and an example of what can go wrong in naive implementations.

For all scenarios, the model error variance had an inverse-gamma (IG) prior with a shape = 0.01 and a scale = 0.01, which are common choices for a vague prior.

PROC MCMC and PROC GENMOD with a BAYES statement were used to sample the posterior in SAS. We employed the default MCMC algorithms within each software, including N-Metropolis (SAS, PROC MCMC), Gibbs for linear regression (SAS, GENMOD + BAYES), random walk Metropolis-Hastings (STATA), and No-U-Turn sampling (RStan, BRMS package).

To improve comparability across results, parameters were initialized with regression coefficients set to “0” and the model error variance set to “1.” In practice, software default settings for initialization values are sufficient. RStan and STATA ran two MCMC chains, while SAS used a single chain. Each chain was run for 10,000 iterations including 1000 iterations discarded as burn-in and all values were stored. Software default settings were used for the target acceptance rate in Metropolis-Hastings algorithms. Convergence was assessed visually using trace, autocorrelation, and histogram plots of the posterior distributions for each parameter. Results were summarized using posterior means and 95% HPD credible intervals (CrI). We also computed the posterior probability of a positive treatment difference represented by a reduction in the mean time to readiness for discharge from PACU (i.e., treatment difference < 0).

We repeated linear modeling in a classical framework in each software. Results were presented as estimates and 95% confidence intervals.

Analysis of Administration of Preoperative Nerve Block

Logistic regression was used to model the association between treatment and the odds of administering preoperative nerve block. Unadjusted and adjusted models were specified.

As above, we considered several priors on the treatment effect. The primary prior was a vague N (0, 10) with 80% prior probability of an odds ratio (OR) between 0.02 and 60. Our optimistic prior was a N (−0.5, 2) with an 80% prior probability of an OR between 0.1 and 3.7 and a 64% prior probability of an OR less than 1. Our skeptical prior was a N (0, 2) with an 80% prior probability of an OR between 0.2 and 6.0 and a 49% prior probability of an OR less than 1. We used vague N (0, 10) for the covariate regression coefficients in the adjusted models. The “pseudo” vague prior was a N (0, 1) for the treatment and covariate regression coefficients. The MCMC algorithms were as above except for SAS GENMOD + BAYES which used a Gamerman algorithm. The posteriors were summarized using means and 95% HPD credible intervals and exponentiated to ORs for a standard clinical interpretation. The posterior probability of a reduced odds of nerve block usage (i.e., OR < 1) was also computed. Analyses were repeated in the classical framework in each software. Results were presented as OR with 95% CIs.