Performance of Bayesian Propensity Score Adjustment for Estimating Causal Effects in Small Clinical Trials

Airi Takagi¹, Takuhiro Yamaguchi²

¹Clinical Research Data Center, Tohoku University Hospital, Sendai, Japan.
²Department of Biostatistics, Tohoku University Graduate School of Medicine, Sendai, Japan.

DOI: 10.4236/ojs.2023.131001                

Abstract

Propensity score (PS) adjustment can control confounding effects and reduce bias when estimating treatment effects in non-randomized trials or observational studies. PS methods are becoming increasingly used to estimate causal effects, including when the sample size is small compared to the number of confounders. With numerous confounders, quasi-complete separation can easily occur in logistic regression used for estimating the PS, but this has not been addressed. We focused on a Bayesian PS method to address the limitations of quasi-complete separation faced by small trials. Bayesian methods are useful because they estimate the PS and causal effects simultaneously while considering the uncertainty of the PS by modelling it as a latent variable. In this study, we conducted simulations to evaluate the performance of Bayesian simultaneous PS estimation by considering the specification of prior distributions for model comparison. We propose a method to improve predictive performance with discrete outcomes in small trials. We found that the specification of prior distributions assigned to logistic regression coefficients was more important in the second step than in the first step, even when there was a quasi-complete separation in the first step. Assigning Cauchy (0, 2.5) to coefficients improved the predictive performance for estimating causal effects and improving the balancing properties of the confounder.

Keywords

Bayesian Estimation, Causal Inference, Propensity Score, Quasi-Complete Separation, Prior Distribution

1. Introduction

Adjustment by propensity score (PS) is a valuable method for comparing treatment effects in non-randomized trials or observational studies. Rosenbaum and Rubin proposed a PS adjustment method to provide an unbiased estimate of the average treatment effect [1] . PS is the treatment assignment probability, calculated from the measured confounding variables using logistic regression. Subjects with the same PS can be considered comparable or exchangeable. Thus, adjustment with PS provides an unbiased comparison between groups.

The conventional frequentist PS adjustment method involves two sequential steps: estimating PS and estimating causal effects. In the first step, PS is the point-estimated probability of treatment assignment. In the second step, causal effects are estimated using the PS calculated in the first step. However, this conventional two-step sequential method does not consider the impact of uncertainty in the PS on the performance of causal effect estimation [2] [3] [4] . Conventional sequential estimations that use regression adjustment with the PS under the correct specification of the causal model in the second step result in a smaller RMSE than other methods, such as stratification, IPTW, and doubly robust estimation [4] . For comparison with the Bayesian simultaneous method, we use the maximum likelihood estimation of the regression model as our conventional two-step method.

The Bayesian method, in which the first and second steps are estimated simultaneously, has been proposed [5] [6] . There has been also a two-step Bayesian estimation method in which the impact of uncertainty in the PS on the performance of causal effect estimation is evaluated [7] . However, the most useful feature of Bayesian estimation is a simultaneous estimation. The simultaneous Bayesian method accounts for the uncertainty in PS estimation by treating the PS as a latent variable during the estimation of causal effects. When using the maximum likelihood estimation of the regression model, if there are insufficient participants for the number of confounders, quasi-complete separation can easily occur in logistic regression used for estimating the PS. Although the Bayesian simultaneous PS method presents a potential solution for the issue of quasi-complete separation faced by small trials, research is lacking in this area. We focus on the Bayesian method to address these issues and compare the performance of Bayesian simultaneous estimation with that of conventional maximum likelihood estimation. In Bayesian estimation, the priors are important. The utility of priors for providing interpretable shrinkage and conducting causal sensitivity analyses has been discussed [8] . Some prior distributions for Bayesian logistic regression have been proposed as useful methods for dealing with separation in discrete outcomes [9] [10] . Simulations were conducted to evaluate the performance of the Bayesian simultaneous estimation assigned to each prior distribution. Furthermore, in Bayesian simultaneous estimation, the first step is affected by the second step through the simultaneous model fitting of PS and causal effects. The balancing properties of PS are important for unbiased comparisons between the treated and untreated groups. Therefore, we also evaluate covariate balancing properties between treatment and control groups [11] [12] [13] [14] .

In Section 2, we describe the conventional frequentist approach and Bayesian simultaneous estimation and outline the methods and prior distributions applied to the data. In Section 3, we present a simulation comparing Bayesian simultaneous and conventional frequentist methods. Data are generated from optional distributions and models in Part 1, and conditions similar to reports of postoperative statin trials comparing the development of acute kidney injury (AKI) between treatment groups in Part 2 [15] . In Section 4, the results of these simulations are reported. In Section 5, we summarize our findings and discuss how the predictive performance of the Bayesian simultaneous method can be improved in small trials.

2. Methodology

2.1. Propensity Score Adjustment

PS is defined as the probability of treatment assignment given measured covariates if the treatment assignment is strongly ignorable and there are no unmeasured confounders. In other words, the PS is estimated under the condition that the treatment assignment cannot be a result of confounding and all confounders have been observed. When estimating PS, a multidimensional vector of measured confounders can be reduced to a one-dimensional vector by using logistic regression. Comparing treated and untreated subjects with the same PS provided an unbiased estimate of the treatment effect. This is similar to the random assignment of the treatment status in clinical trials.

The PS adjustment for a binary outcome is described below with treatment variable X, outcome variable Y, a vector of p covariates plus intercept C (1, C₁, …, C_p), and regression coefficients vector γ:

$\mathrm{logit}\left[\mathrm{Pr}\left(X=1|C\right)\right]=\gamma C$ (1)

and

$\mathrm{logit}\left[\mathrm{Pr}\left(Y=1|X,C\right)\right]=\alpha X+\beta g\left(C,\gamma \right).$ (2)

First, the PS is estimated using a logistic regression model that regresses the treatment assignment on observed confounders, as described in Equation (1). Subsequently, the treatment effect is regressed on the treatment assignment and PS, as described in Equation (2).

Equation (1) defines the association between X given C, and Equation (2) defines the association between Y and X given C. g(C, γ) in Equation (2) is the PS. In the simultaneous Bayesian method, posterior estimates of the parameters calculated from Bayesian logistic regression were obtained using the models given in Equations (1) and (2). The PS is modeled as a latent variable and the marginal posterior for the treatment effect is integrated over this variable. In this study, we sample from posterior distributions using the Metropolis-Hastings MCMC algorithm. Updating γ from its conditional posterior distribution corresponds to the simultaneous updating of PS, α, and β in the outcome model.

2.2. Separation and Prior Information

2.2.1. Quasi-Complete Separation

Quasi-complete separation occurs when one or more covariates almost perfectly predict the outcome variable in models for binary outcomes.

Consider logistic regression for PS estimation with the dependent variable X_i and a vector of covariates C_i with coefficient vector γ:

$\mathrm{Pr}\left(X_i=1|C_i,\gamma \right)=\frac{1}{1+\exp \left(-C_i\gamma \right)}.$

The log-likelihood for this model is:

$\ln L\left(\gamma |X\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}\left\{X_i\ln \left[\frac{1}{1+\exp \left(-C_i\gamma \right)}\right]+\left(1-X_i\right)\ln \left[1-\frac{1}{1+\exp \left(-C_i\gamma \right)}\right]\right\}}.$

In quasi-complete separation, the maximum likelihood estimate γ is large. Separation is related to little or no overlap in covariate distributions between the treatment groups.

2.2.2. Firth Bias Reduction

A penalized likelihood-based Firth logistic regression method may provide an effective solution to separation [16] . Firth’s approach adds a bias term to the log-likelihood function as follows:

$\ln L{\left(\gamma |X\right)}^{*}=\ln L\left(\gamma |X\right)+\frac{1}{2}\ln \left|I\left(\gamma \right)\right|.$

This penalized likelihood eliminates bias even in the presence of complete separation.

2.2.3. Cauchy Distribution

In Bayesian linear regression, the prior distributions of the regression coefficients are important. Non-informative priors are often applied when there is no prior information. In extreme conditions such as separation, Lasso regression is equivalent to using the Laplace prior in Bayesian interpretation [17] . The Laplace distribution is centered around zero and has long tails. When applying such distributions to the regression coefficients, the posterior estimates shrink towards zero. Cauchy priors have features similar to those of the Laplace prior, whereas the normal prior is more diffuse around zero.

Gelman recommends the Cauchy distribution with a center of 0 and a scale of 2.5 as a default prior for regression coefficients in logistic regression because it allows inferences to be drawn even in the presence of separation [9] . Previous research has found that logistic regression coefficients are typically between −5 and 5. In logistic regression, a change of 5 corresponds to moving a probability from 0.01 to 0.5, or from 0.5 to 0.99 [9] . Here, we used a Cauchy distribution with a center of 0 and a scale of 2.5 (Cauchy (0, 2.5)).

2.2.4. Inclusion of Data Information

There are several approaches to the inclusion of data information in prior distributions for Bayesian analysis. Sullivan and Greenland calculated prior information, such as the median or variance from confidential intervals [18] [19] . In this study, we use a normal prior distribution that has information arising from true models; the mean is the true parameter value and the variance is 6.25, because the coefficients of logistic regression are almost always between −5 and 5, as described above.

3. Application

Simulations were conducted to evaluate the performance of the Bayesian simultaneous method. The simulation study comprised two parts. In “Part 1”, we generated data from arbitrary distributions to confirm the basic performance of each model and prior distribution. In “Part 2”, we generated data using a similar condition to that reported in the postoperative statin trial [15] . Each model and prior distribution were applied to the generated data. Simulations were performed using SAS, version 9.4 (SAS Institute, Cary, NC, USA).

3.1. Data Generation

In Part 1, we generated replicated datasets containing N = 160 observations for small trials, which led to quasi-complete separation, and N = 1000 observations for large trials, which is sufficient to eliminate quasi-complete separation. All the simulated datasets contain eight covariates C_i (C₁, …, C₈), four of which are continuous variables, and the remaining are discrete. Continuous variables were generated from a normal distribution with a mean of 0 and variance of 1, and discrete variables were assumed to follow a Bernoulli distribution with p = 0.3. X_i and Y_i are generated from Bernoulli distributions, and the relationship between them is described below.

$\mathrm{Pr}\left(X_i=1|C_i,\gamma \right)=\frac{\exp \left({\gamma }_0+{\gamma }_1{C}_{i1}+\cdots +{\gamma }_8{C}_{i8}\right)}{1+\exp \left({\gamma }_0+{\gamma }_1{C}_{i1}+\cdots +{\gamma }_8{C}_{i8}\right)}$ (3)

and

$\mathrm{Pr}\left(Y_i=1|X_i,C_i\right)=\frac{\exp \left({\alpha }_0+{\alpha }_1{X}_i+{\beta }_1{C}_{i1}+\cdots +{\beta }_8{C}_{i8}\right)}{1+\exp \left({\alpha }_0+{\alpha }_1{X}_i+{\beta }_1{C}_{i1}+\cdots +{\beta }_8{C}_{i8}\right)}.$ (4)

We used logistic regression in both steps and regarded Equations (3) and (4) as the true models used to generate data under the given simulation design. Table 1 summarizes the characteristics of each simulation and the true values of the model coefficients. The simulation scenario and regression coefficients in Part 1 follow in Table 1. We assigned each value to the coefficient γ_n in Equation (3). Each scenario had different coefficient values (α₀, α₁, β_n) in Equation (4). We set γ₈ = −4.4 as the coefficient of confounder 8 to cause quasi-complete separation within 0 < g(C, γ) < 1 and 0 < Pr(Y_i = 1|X_i, C_i) < 1. The allocation rate was set to 1:1. Under this condition, quasi-complete separation occurred in 16% of replicated datasets. This is not severe and can occur in clinical trials with small sample sizes. Using the conditions for small trials, we generated 8000 datasets

Table 1. Catalog of simulation designs including regression coefficient, outcome event rate per treatment group, number of subjects, and separation rate in replicated datasets. Each row corresponds to the true values of the regression coefficient used to generate the given simulation design. X = 1 and X = 0 are set for each treatment group.

and sampled 1000 datasets that led to quasi-complete separation (small (100%)) and another 1000 that did not (small (0%)). Using the conditions for large trials, we generated 1000 datasets (large).

In Part 2, data were generated using conditions similar to those reported in a postoperative statin trial comparing AKI development between patients who used postoperative statins and those who did not [15] . In this trial, the sample size was insufficient, given the number of confounding variables, and quasi-complete separation occurred in 13% of the replicated datasets. From 10,000 generated datasets, we sampled 1000 datasets with quasi-complete separation (100%) and another 1000 without (0%). The replicated datasets contained 324 observations and 19 covariates (244 patients taking postoperative statins and 80 patients not taking postoperative statins). The outcome event rate was 21%. The association between outcome events and treatment variables, given confounders, and the distribution of confounders for patients taking postoperative statins and those who developed AKI are described in the original report. From 10,000 generated datasets, we sampled 1000 datasets that led to quasi-complete separation (100%) and another 1000 that did not (0%).

3.2. Prior Distributions and Models Applied to Replicated Datasets

In the Bayesian simultaneous approach, we assigned a normal distribution with a mean of 0 and variance of 10,000, a normal distribution with mean γ_n using information from the true parameter values of regression coefficients and variance of 6.25, and a Cauchy distribution with a center of 0 and a scale of 2.5. Table 2 lists the models applied to the data. We sampled from the posterior density using the Metropolis-Hastings algorithm and used posterior mean estimates obtained from MCMC chains of length 10,000 with the first 2000 discarded as burn-in. For comparison, we used the maximum likelihood estimation of the regression model (ML) and Firth penalized likelihood approach (Firth) as the sequential method.

Table 2. Catalog of simulation designs including the number of subjects, propensity score distributions between each group, separation rate of 1000 replicated datasets, regression coefficients, and the association between treatment and confounders (ML is a maximum likelihood estimation model).

3.3. Analysis

We evaluated the numerical performance of estimating the causal effects of 1000 replicated datasets through bias, MSE, and 95% coverage. Causal effects were obtained from estimators of the regression coefficient α₁, which is the log odds ratio of the treatment effect. We referred to 95% confidence intervals or Bayesian credible intervals for calculating coverage.

We checked the balancing properties of the PS by using standardized differences to investigate whether the resulting covariate balance was optimized [11] . We conducted PS adjustment using regression with each confounder and PS as independent variables. We assessed the degree of imbalance in the distributions of each confounder between the treatment groups after PS adjustment using score d.

For continuous variables, balancing score d is defined as:

$d=\frac{\left|{\bar{C}}_{X=1}-{\bar{C}}_{X=0}\right|}{\sqrt{\frac{S{D}_{X=1}^2+S{D}_{X=0}^2}{2}}}\times 100\left(\%\right).$ (5)

where ${\bar{C}}_{X=1}$ and ${\bar{C}}_{X=0}$ denote the sample mean of the confounder in each treatment group and $S{D}_{X=1}^2$ and $S{D}_{X=0}^2$ denote the sample variance of the confounder in each treatment group. For discrete variables, the score d is defined as:

$d=\frac{\left|{\hat{p}}_{X=1}-{\hat{p}}_{X=0}\right|}{\sqrt{\frac{{\hat{p}}_{X=1}\left(1-{\hat{p}}_{X=1}\right)+{\hat{p}}_{X=0}\left(1-{\hat{p}}_{X=0}\right)}{2}}}\times 100\left(\%\right).$ (6)

where ${\hat{p}}_{X=1}$ and ${\hat{p}}_{X=0}$ denote the prevalence of the discrete variables in each treatment group.

4. Results

4.1. Impact of Quasi-Complete Separation and Model Comparison

Figure 1 shows the distributions of the causal effects. With a large sample size (N = 1000), the causal effect estimation variability is smaller than that with a small sample size (N = 160). The numerically summarized estimates of α₁ are presented in Table 3. Throughout the analysis, estimates of α₁ under Scenarios 1 and 4 present a large bias and MSE compared with those under Scenarios 2 and 3. Point estimates of α₁ using conventional frequentist estimation had the lowest MSE, but the coverage was slightly higher than 0.95. When using the Bayesian method with small sample sizes, Cauchy 2 shows superior performance relative to the other models, even though TrueCoef uses information from true coefficients. However, shrinking towards zero seems to lead to bias. With a large sample size, the distributions of α₁ estimates were less variable, regardless of the model and method of estimation used. In contrast, with small sample sizes, estimates of α₁ show large distributions and MSE regardless of the presence of quasi-complete separation. Bias depends on each scenario.

Figure 1. Distributions of estimate of causal effect (α₁) calculated from traditional sequential (ML, Firth) and Bayesian simultaneous method (Normal, Cauchy 2) under Scenario 2 and 3. Small-sized datasets contain 160 observations and Large-sized datasets contain 1000. Small (100%) shows all of replicated datasets lead to separation and Small (0%) shows none of them lead to separation. Horizontal dotted lines are the true parameter values.

Table 3. Numerical performance comparison of estimates of causal effect (α₁) to assess the impact of prior distribution assigned in the first step. Causal effect (α₁) is estimated from conventional sequential (maximum likelihood of regression, ML and Firth) and Bayesian simultaneous methods (Normal, Cauchy 1) with different specifications of prior distributions in the first step and the same specifications in the second step. Small datasets contain 160 observations and large datasets contain 1000 observations. “Small (100%)” shows that all of replicated datasets lead to quasi-complete separation and “Small (0%)” shows that none do.

Table 3 compares the impact of quasi-complete separation in the first step on each model. The performance of α₁ estimation is similar to that of the conventional PS estimation (ML) and Firth penalized method (Firth) despite quasi-complete separation in the first step. In the Bayesian simultaneous method, regardless of the scale parameter or shape of the prior distributions in the first step, there is little impact on the posterior estimates of α₁. The performance of models with the same specification of prior distributions in the second step (Cauchy 1, Normal) was similar, even if they had different specifications in the first step. Table 4 illustrates the performance of the models assigned different specifications for prior distributions in the second step (Cauchy 1, Cauchy 2, TrueCoef). A reduction in the scale parameter of the prior distributions in the second step led to an overall reduction in MSE and an improvement in 95% coverage. Estimates of α₁ are only slightly influenced by the specification of prior distributions in the first step but have a large impact from the second step.

4.2. Assessment of Covariate Balance of Models

Figure 2 shows the boxplot of the score d of confounder 8, which leads to quasi-complete separation, to assess the balancing properties of PS adjustment. Regardless of the prediction model used, PS adjustment significantly reduced the score d. However, the score d from the Bayesian simultaneous methods was larger than that from the conventional frequentist method. Additionally, the score d calculated from the Bayesian simultaneous methods varied by scenario, whereas those from the frequentist methods were constant.

The score d from Cauchy 1 when assigning Cauchy (0, 2.5) in the first step was superior to when using Normal (0, 10,000). Additionally, comparisons between Cauchy 1 and Cauchy 2 appear to show that the reduction in the scale parameter of prior distributions assigned in the second step is associated with an improvement in balancing properties.

4.3. Model Comparison in Postoperative Statin Use Trial

The results of the simulation using the conditions reported in the postoperative statin trial are described in Table 5. The true value of the odd ratio of causal

Table 4. Numerical performance comparison of estimates of causal effect (α₁) to assess the impact of prior distribution assigned in the second step. Causal effect (α₁) is estimated from conventional sequential (maximum likelihood, ML) and Bayesian simultaneous methods (Cauchy 1, Cauchy 2, TrueCoef) with different specifications of prior distributions in the second step and nearly identical specifications in the first step. Small datasets contain 160 observations and large datasets contain 1000 observations. “Small (100%)” shows that all of replicated datasets lead to quasi-complete separation and “Small (0%)” shows that none do.

Figure 2. Distributions of balancing score d of confounder 8 which is cause of separation. Each regression coefficient was adjusted PS from traditional sequential (ML, Firth) and Bayesian simultaneous method (Normal, Cauchy 2, TrueCoef) and calculated balancing properties. Small-sized datasets contain 160 observations and Large-sized datasets contain 1000. Small (100%) shows all of replicated datasets lead to quasi-complete separation and Small (0%) shows none do.

Table 5. The consequence of estimates of causal effect (α₁) from conventional sequential (maximum likelihood, ML) and Bayesian simultaneous (Normal, Cauchy 2) methods in a postoperative statin trial setting. Replicated datasets contain 324 observations; “All (100%)” shows all of replicated datasets lead to quasi-complete separation and “None (0%)” shows none of them lead to quasi-complete separation.

effects (0.3) was obtained from the trial report. Quasi-complete separation occurred in 12% of the replicated data. Assigning Cauchy 2 yielded a similar performance to that of conventional ML.

5. Discussion and Conclusions

In this study, we conducted a simulation assuming a small trial with quasi-complete separation caused by one specified confounder and compared the performance of each model.

With a small sample size, despite the presence of quasi-complete separation and using Firth’s bias-eliminated estimation in the first step, the performance of the causal effect estimation is only slightly affected. Dimensionality reduction and sequential estimation in frequentist PS methods circumvent the adverse impact of a high-dimensional vector of confounders. Additionally, regression adjustment by PS exhibits superior performance if the outcome model is true [4] . In the Bayesian simultaneous method, the shrinkage estimation of regression coefficients in the first step has little impact on the performance of the causal effect estimation because calculations are fitted to the model of the second step [5] [6] . Hence, prior distributions applied to the coefficients in the second step are important, and the variability of α₁ is influenced by the specification of prior distributions in this step. In contrast, with large sample sizes, the impact of giving prior specifications is small. An insufficient sample size with respect to the number of confounders produced large variability in the PS estimates and led to a remarkable increase in the MSE of α₁. Therefore, an adequate sample size should be used.

Using Cauchy 2 drastically reduces the variance of causal effect estimates compared with other models, but the value of bias depends on the scenario. An adequate scale of prior distributions and shrinkage estimation in the second step seems to reduce the variability, and assigning Cauchy (0, 2.5) may be valuable regarding variance and coverage. However, further research that considers the trade-off between bias and variability is required. Using Cauchy 2 is superior to using TrueCoef because the normal distribution tends to be more diffuse around the mean [10] [17] . The conventional frequentist method (ML) exhibits a small bias and low MSE, but the coverage of α₁ is observed slightly higher than 0.95, without considering residuals.

The balancing property of confounder 8, which led to quasi-complete separation, significantly improved in all models. When using the Bayesian method, Cauchy 2 is superior because the shrinkage estimation of the coefficients caused by Cauchy (0, 2.5) improves the balancing properties. Score d calculated from posterior estimates varies depending on the scenario because simultaneous fitting for the model in the first and second steps provides flexibility in the estimation of γ_n. However, the first step is affected by the second step because it must be adjusted for outcome model fitting such that the posterior distribution of γ_n involves additional information from the second step. It has been argued that variables included as measured confounders should not depend on the outcome model [12] .

In small trials, the variability of causal effect estimation is large, without considering quasi-complete separation; therefore, it is necessary to select adequate prediction models and sample sizes related to the number of confounders. The Bayesian simultaneous PS method, assigning Cauchy (0, 2.5) to coefficients as the prior distribution, is superior in terms of predictive performance for estimating causal effects and improving the balancing properties of the confounder, when leading to quasi-complete separation.

Acknowledgements

I wish to thank Department of Biostatistics of Tohoku University Graduate School of Medicine for their support and encouragement.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Rosenbaum, P.R. and Rubin, D.B. (1983) The Central Role of the Propensity Score in Observational Studies for Causal Effects. Biometrika, 70, 41-55. https://doi.org/10.1093/biomet/70.1.41
[2]	Austin, P.C. and Mamdani, M.M. (2005) A Comparison of Propensity Score Methods: A Case-Study Estimating the Effectiveness of Post-AMI Statin Use. Statistics in Medicine, 25, 2084-2106. https://doi.org/10.1002/sim.2328
[3]	Gelman, A. and Hill, J. (2007) Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, New York. https://doi.org/10.1017/CBO9780511790942
[4]	Lunceford, J.K. and Davidian, M. (2004) Stratification and Weighting via the Propensity Score in Estimation of Causal Treatment Effects: A Comparative Study. Statistics in Medicine, 23, 2937-2960. https://doi.org/10.1002/sim.1903
[5]	McCandless, L.C., Gustafson, P. and Austin, P.C. (2009) Bayesian Propensity Score Analysis for Observational Data. Statistics in Medicine, 28, 94-112. https://doi.org/10.1002/sim.3460
[6]	Zigler, C.M., et al. (2013) Model Feedback in Bayesian Propensity Score Estimation. Biometrics, 69, 263-273. https://doi.org/10.1111/j.1541-0420.2012.01830.x
[7]	Liao, S.X. and Zigler, C.M. (2020) Uncertainty in the Design Stage of Two-Stage Bayesian Propensity Score Analysis. Statistics in Medicine, 39, 2265-2290. https://doi.org/10.1002/sim.8486
[8]	Oganisian, A. and Roy, J.A. (2021) A Practical Introduction to Bayesian Estimation of Causal Effects: Parametric and Nonparametric Approaches. Statistics in Medicine, 40, 518-551. https://doi.org/10.1002/sim.8761
[9]	Gelman, A., Jakulin, A., Pittau, M.G. and Su, Y.S. (2008) A Weakly Informative Default Prior Distribution for Logistic and Other Regression Models. Annals of Applied Statistics, 2, 1360-1383. https://doi.org/10.1214/08-AOAS191
[10]	Greenland, S. and Mansournia, M.A. (2015) Penalization, Bias Reduction, and Default Priors in Logistic and Related Categorical and Survival Regressions. Statistics in Medicine, 34, 3133-3143. https://doi.org/10.1002/sim.6537
[11]	Austin, P.C. (2009) Balance Diagnostics for Comparing the Distribution of Baseline Covariates between Treatment Groups in Propensity-Score Matched Samples. Statistics in Medicine, 28, 3083-3107. https://doi.org/10.1002/sim.3697
[12]	McCandless, L.C., Gustafson, P., Austin, P.C. and Levy, A.R. (2009) Covariate Balance in a Bayesian Propensity Score Analysis of Beta Blocker Therapy in Heart Failure Patients. Epidemiologic Perspectives & Innovations, 6, Article No. 5. https://doi.org/10.1186/1742-5573-6-5
[13]	Rosenbaum, P.R. and Rubin, D.B. (1985) Constructing a Control Group Using Multivariate Matched Sampling Methods That Incorporate the Propensity Score. Journal of the American Statistical Association, 39, 33-38. https://doi.org/10.1080/00031305.1985.10479383
[14]	Rosenbaum, P.R. (1989) Optimal Matching for Observational Studies. Journal of the American Statistical Association, 84, 1024-1032. https://doi.org/10.1080/01621459.1989.10478868
[15]	Billings 4th., F.T., Pretorius, M., Siew, E.D., Yu, C. and Brown, N.J. (2010) Early Postoperative Statin Therapy Is Associated with a Lower Incidence of Acute Kidney Injury Following Cardiac Surgery. Journal of Cardiothoracic Vascular Anesthesia, 24, 913-920. https://doi.org/10.1053/j.jvca.2010.03.024
[16]	Firth, D. (1993) Bias Reduction of Maximum Likelihood Estimate. Biometrika, 80, 27-38. https://doi.org/10.1093/biomet/80.1.27
[17]	Tibshirani, R. (1996) Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society: Series B, 58, 267-288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x
[18]	Greenland, S. (2006) Bayesian Perspectives for Epidemiological Research: I. Foundations and Basic Methods. International Journal of Epidemiology, 35, 765-775. https://doi.org/10.1093/ije/dyi312
[19]	Sullivan, S.G. and Greenland, S. (2013) Bayesian Regression in SAS Software. International Journal of Epidemiology, 42, 308-317. https://doi.org/10.1093/ije/dys213