Leveraging a Bayesian Approach in a Comparative Effectiveness Trial of Major Adverse Cardiovascular Events

Introduction

There have been several pivotal studies of sodium–glucose cotransporter-2 inhibitors (SGLT2i) and dipeptidyl peptidase 4 inhibitors (DPP4i) that evaluated major adverse cardiovascular events and heart failure hospitalization.^1–12 Several studies have shown benefit in Major Adverse Cardiovascular Events and Heart Failure Hospitalization (MACE+HF) among those with preexisting cardiovascular disease (CVD) (secondary prevention), but results were more varied for those without a history of cardiovascular disease (primary prevention). Multiple cohort evaluations which utilized Cox proportional hazards analysis failed to find statistically significant effects and others found protective associations.^5–12 Further, in a population without a history of heart failure, the CANVAS trial failed to identify a strong protective effect of SGLT2i (hazard ratio (HR): 0.87, 95% confidence interval (CI): 0.72, 1.06), but it did identify a strong reduction in MACE outcomes for those with a history of heart failure (HR: 0.61, 95% CI: 0.46, 0.80).^5,8 Additionally, previous studies have estimated different effect estimates for SGLT2i for each type of cardiovascular event. A meta-analysis of SGLT2i trials on the composite of myocardial infarction, stroke, and cardiovascular death estimated an HR of 1.00 with a 95% CI of (0.87, 1.16) for a fixed effects model on patients with multiple risk factors. When focusing on heart failure (HF) hospitalizations and cardiovascular death, that same meta-analysis estimated an HR of 0.84 with a 95% CI of (0.69, 1.01) for a fixed effects model on patients with multiple risk factors.⁵

We previously reported that among a cohort of Veterans without a history of cardiovascular disease (primary prevention), SGLT2i for diabetes treatment was associated with reduced risk of heart failure hospitalization when compared with DPP4i. When performing a Cox proportional hazard analysis using propensity score weighting to balance covariates at baseline, the HR for SGLT2i versus dipeptidyl peptidase 4 inhibitors (DPP4i) for time to HF hospitalization only was estimated to be 0.73 with 95% CI of (0.54,0.97). When considering the composite outcome of MACE+HF, the HR for SGLT2i versus DPP4i was of 0.87 (95% CI: 0.74, 1.03).²

The inclusion of the null effect in the confidence interval for our previous study was met with skepticism as this result conflicted with previous studies that found a statistically significant beneficial association between SGLT2i and cardiovascular outcomes, but were restricted to a different population with cardiovascular disease.^7,9 Additionally, we were concerned that some readers misinterpreted the Cox proportional hazard results for the composite outcome as a null effect of SGLT2i on cardiovascular outcomes. The failure to reject the null hypothesis does not necessarily provide strong support for the null hypothesis. The oversimplification that the null value in the 95% confidence interval demonstrates proof of the null persists even among well-trained medical professionals.

To address this skepticism for null effects and potential misinterpretations, we apply a Bayesian approach to further investigate the association of SGLT2i with MACE+HF and its components in a propensity score matched primary prevention cohort. The Bayesian approach allowed us to further investigate the differences in effect estimates seen for the individual components of the MACE +HF outcome (HF hospitalization, CVD hospitalization, and CVD death) using probability statements as opposed to hazard ratios and 95% confidence intervals, which can be misinterpreted. It also allowed us to incorporate prior knowledge and skepticism regarding a null effect using informative prior distributions.

Bayesian statistics assume that model parameters have an unknown distribution rather than one true value and use Bayes’ Rule to incorporate prior knowledge about this distribution. Inference is made using the posterior distribution of the parameters, which is proportional to the product of a prior distribution and a likelihood distribution.^13,14 Prior information can be included via the choice of the prior distribution. If we are confident about the range of potential values for a parameter of interest, a highly informative prior with a small variance can be used. The posterior distribution will not be heavily influenced by new data, and the posterior distribution will resemble the prior distribution. On the other hand, in scenarios where the investigators are unsure about the values a parameter will take, the use of a less informative or “uninformative prior” with a larger variance and flatter shape can be chosen. The posterior is influenced to a greater degree by the new data and will resemble the likelihood more than the prior distribution.

A major benefit of the Bayesian approach is that it can provide an intuitive summary of the results, as the posterior distribution allows one to make probability statements about the parameters. For example, the Bayesian approach allows us to calculate 95% credible intervals by calculating the 0.025 and 0.975 percentile of the posterior distribution for a parameter of interest. We can then say there is a 95% probability that the true parameter value lies within the interval. Conversely, the interpretation of a frequentist confidence interval is not a probability statement.^13,14

We apply Bayesian survival modeling to a large observational propensity score matched primary prevention cohort of United States Veterans to explore the association of SGLT2i with MACE+HF and its components. We incorporate prior beliefs and skepticism for a null association via an informative prior. We use the posterior distribution to estimate the probability of a meaningful protective association between SGLT2i and MACE+HF events and its components to provide an easily interpretable summary measure.

Methods

Data Set for the Bayesian Analysis

The Bayesian analysis is performed on a propensity score matched cohort obtained from a cohort of Veterans Health Administration (VHA) patients with diabetes mellitus (DM) who had a first prescription for a hypoglycemic medication between January 1, 2001, and December 31, 2016. Additional cohort data were obtained through December 31, 2019. The results and details of this cohort have been previously published.² The study cohort includes Veterans aged 18 years or older with DM who were using metformin, sulfonylurea, or insulin alone, or in combination and later fill one of the following newer drug classes: SGLT2i and DPP4i. We further subset to patients without a history of cardiovascular diseases. The institutional review board of the VHA approved this study with a waiver of informed consent. In compliance with the Declaration of Helsinki, the research involved no more than minimal risk as it was a Phase 4 study of real prescribing patterns. The data accessed was limited to variables regularly stored in medical records and no individual participant was identifiable as the results were reported in aggregate. Patient data was kept on a secure VHA server, and the analysis was performed on anonymized data on secure VHA servers. The research could not have been carried out without the waiver as it was a retrospective study of a large number of Veterans and current contact information would likely be unavailable for a meaningful proportion of Veterans in the study.

The primary outcome was the time to the MACE + HF event. The outcome date was the hospital admission date for acute myocardial infarction, ischemic or hemorrhagic stroke, acute heart failure, or cardiovascular death date. The primary discharge diagnosis or underlying cause of death identified each event. We also evaluated each component separately, not censoring on the other components, except for cardiovascular death. The start of follow-up was the medication fill date of SGLT2i or DPP4i. The follow-up continued until an outcome or censoring event where censoring criteria included: death, study end date (December 31, 2019); non persistence (90 days without medication); crossover/addition of diabetes drug in a different class (eg, SGLT2i user who starts DPP4i) and the 181st day of no VHA contact (inpatient, outpatient, or pharmacy use).²

Propensity Score Matching

Propensity score matching, which is commonly used in large EHR-based studies, was used here as a data-preprocessing method to improve model performance and achieve covariate balance at baseline between the SGLT2i and DPP4i groups.^15–17 The propensity score (PS) is the conditional probability of assignment to a treatment given a set of covariates.¹⁸ The propensity score is often calculated using a logistic regression model: , where is the propensity score, Z is the exposure with 1 being exposed and 0 being unexposed, and X is a set of observed baseline characteristics. Once PS are calculated, they are used in varied manners as a data reduction technique and to achieve balance within a population. PS strategies including matching, stratification, inverse probability of treatment weighting (IPTW), and as a measure for covariate adjustment.^15,19,20 In this analysis, we implemented propensity score matching using the nearest-neighbor approach and a caliper equal to 0.05, as it was simple to implement and adequately balanced the groups. To perform nearest neighbor matching, treated subjects are ordered randomly and are matched one at a time to the control subject with the closest propensity score. Once a pair is matched, they are removed from the matching pool. Matching continues until all treated subjects are matched or no additional acceptable matches can be made. Remaining observations are dropped if the difference in propensity scores is over the caliper.^21,22 As a comparator to the Bayesian time-to-event analysis, we ran a Cox proportional hazards regression analysis on this newly created propensity score matched cohort.

Bayesian Time-to-Event Model

Like the Cox proportional hazards model, the proposed Bayesian approach also uses a hazard model with direct covariate adjustment. We assumed the following hazard function for episode at time :

where

is a vector of the values of the baseline covariates for episode , and is a vector of the corresponding coefficients.²³ Also, is the baseline hazard function evaluated at and is estimated using monotone splines, also known as M-splines. M-splines are non-increasing or non-decreasing functions and are monotone piecewise polynomials that can be expressed in the following way:

where denotes the th basis term for a degree M-spline function evaluated at a vector of knot locations, , and is the lth M-spline coefficient.²⁴ M-splines capture a non-linear but monotonic relationship, meaning that the hazard is either non-increasing or non-decreasing. They are flexible and allow for the approximation of complex functions by piecing together simpler polynomial functions in a smooth and continuous manner.^23,24

We assumed and let be an indicator of whether the episode is assigned to the treatment group, which in our case is SGLT2i versus DPP4i. Like in the Cox model, is the cause-specific hazard ratio (HR) for the treatment effect. A HR less than one suggests a protective association (lower association with the event in the next time period compared to a control subject). The posterior distribution for the HR was found by exponentiating posterior draws for . Bayesian modeling was performed using R Statistical Software²⁵ (version 4.4.1) and the stan_surv() function in the rstanarm package.^23,26 These packages fit the Bayesian model using an implementation of the Hamiltonian Monte Carlo method known as the No-U-Turn Sampler (NUTS).²⁷

Specifying the Prior Distribution

We assumed fixed values for and . We fixed to be a vector with knots at the minimum timepoint, maximum timepoint and 3 internal knots at equal spaced percentiles, 25%, 50% and 75%. We fixed . For all parameters besides , we used the default suggestions in the stan_surv() function in the rstanarm package in R.^23,26 A normal prior with mean 0 and a standard deviation of 20 is used for the intercept, . For the M-spline coefficients, , a Dirichlet prior with hyperparameter of a vector with all elements equal to 1 is used.

Priors Based on Prior Studies

For , we can incorporate information from other studies by using a normal prior with the mean centered around the log of an estimate of HR in the literature and pick a standard deviation based on the corresponding confidence interval for HR. We performed a literature search and focused on a meta-analysis of cardiovascular outcome trials.⁵ This meta-analysis combined data from three cardiovascular outcome trials, finding a hazard ratio of 1.00 (0.87,1.16) in patients without cardiovascular disease and a hazard ratio of 0.86 (0.80,0.93) in patients with cardiovascular disease (CVD).⁵ Based on the hazard ratio 1.00 (0.87,1.16), we created an informative prior with normal distribution, a mean of , and a standard deviation (sd) of 0.220, referred to as the No CVD prior. The sd was first approximated based on the confidence interval (0.87, 1.16), and then tripled to allow for more variability in posterior draws.⁵ The sd based on (0.87,1.16) was approximated by taking the mean of and and multiplying by 3 to get sd . We include the No CVD prior because it is a summary value from a systematic review of the effect of SGLT2i in studies done on patients without cardiovascular disease, which was our population of interest.

We also created a prior based on 0.86 (0.80,0.93) and refer to it as the CVD prior. Using a similar approach, we selected a normal distribution with mean of and sd=. Although the subjects in our population do not have a history of cardiovascular disease, we included this prior because it incorporates information from published studies that suggest a protective effect of SGLT2 for the MACE+HF outcome.

Uninformative Normal Prior

The last prior considered in our Bayesian survival model is a normal distribution with a large variance (sd = 25) centered at 0. We called this prior the Uninformative prior as its shape is relatively flat compared to the other priors and will have a smaller influence on the posterior distribution of and HR. In addition to fitting a Bayesian survival model for the MACE+HF outcome using the Uninformative prior, we evaluated separate Bayesian models using the Uninformative prior for the cause-specific hazard for each individual component of the MACE+HF outcome: heart failure hospitalization (HF), cardiovascular hospitalization (CVD Hospitalization), and cardiovascular death (CVD Death).

Probability of a Meaningful Protective Association

The time-to-event Bayesian model allows estimation of the probability of a meaningful protective association by using the posterior distribution of HR = . We consider an HR less than 0.90 to be meaningful protective association based on the meta-analysis with estimate of HR = 0.86 for the cardiovascular cohort.⁵ The estimated probability that HR < 0.90 can be found by taking the number of posterior draws where and dividing by the number of posterior draws.

Results

Analytic Cohort

In the retrospective cohort, 129,834 episodes were included with 23,107 episodes for SGLT2i medication and 106,727 assigned DPP4i (reference). After fitting a logistic regression propensity score model using the covariates given in Table 1 and performing one-to-one matching using a caliper value of 0.05, we matched 21,821 episodes of SGLT2i to 21,821 episodes of DPP4i. A total of 1286 SGLT2i episodes were dropped due to the difference in propensity scores being above the caliper. The distributions of the propensity scores for SGLT2i and DPP4i groups are given in Supplemental Figure 1. After matching, all standardized mean differences of observed baseline covariates were below the threshold of 0.2 and all but one were at or below 0.05, indicating good balance between the groups. The study cohort prior to and after matching is summarized in Table 1.

Table 1 Summary of Patient Characteristics Prior to and After Propensity Score Matching

Primary Outcome: MACE+HF

In the matched cohort, there were 188 MACE+HF events for SGLT2i and 277 MACE+HF for DPP4i. The distributions of the three priors for (No CVD, CVD, and Uninformative) used in the Bayesian time-to-event analysis are shown in Supplemental Figure 2. For each prior for , Figure 1 demonstrates the posterior distribution for the cause-specific hazard ratio, , for SGLT2i. The prior distributions influence both the location and shape of the corresponding posterior distributions. The estimated hazard ratios and corresponding 95% credible intervals differed by prior: No CVD: 0.82 (0.68,0.96), CVD: 0.81 (0.71,0.94), and Uninformative: 0.79 (0.65,0.94) but all suggest a protective association with the MACE+HF outcome (Table 2). The estimates of the HR and 95% credible intervals for the Bayesian analysis using the Uninformative prior were similar to the results found performing the Cox proportional hazards regression with an HR of 0.78 and 95% CI of (0.65, 0.95). The estimated values for the probability of a meaningful protective association also support a protective association and the probability of a meaningful protective association for the Bayesian model with the No CVD, CVD, and Uninformative priors were 88%, 92%, and 93%, respectively (Table 2). Focusing on the Bayesian time-to-event model with the Uninformative prior, the estimated survival function demonstrated greater survival probability for the SGLT2i for the MACE+HF outcome (Figure 2A).

Table 2 Event Rates and Cause-Specific Hazard Ratios for Major Adverse Cardiovascular Events Plus Heart Failure

Figure 1 Posterior Distributions for the Cause-specific Hazard Ratio. The posterior distributions for the cause-specific hazard ratio for sodium–glucose cotransporter-2 inhibitors (SGLT2i) versus dipeptidyl peptidase 4 inhibitors (DPP4i) for the major adverse cardiovascular events plus heart failure (MACE+HF) outcome from the three different Bayesian models (No CVD: yellow, CVD: blue, Uninformative: black). The vertical dashed lines (No CVD: yellow, CVD: blue, Uninformative: black) denote the posterior means, and the vertical grey line denotes a null association.

Figure 2 Survival Plots. Probability of survival for major adverse cardiovascular events plus heart failure ((A): MACE+HF) and its components ((B): HF, (C): CVD Hospitalization, and (D): CVD Death) among sodium–glucose cotransporter-2 inhibitors (SGLT2i) vs dipeptidyl peptidase 4 inhibitor (DPP4i) users without cardiovascular disease from the cause-specific Bayesian time-to-event analyses using the Uninformative prior. The blue line and blue shading provide the survival function and 95% pointwise prediction interval for SGLT2i users, respectively, and the black line and grey shading provides the survival function and 95% pointwise prediction interval for DPP4i users, respectively.

Components of MACE+HF

In the matched cohort, there were 57 heart failure hospitalization (HF) events for SGLT2i and 98 HF for DPP4i (Table 3). For the Bayesian time-to-event model with the Uninformative prior and the HF outcome, the cause-specific hazard ratio estimate and 95% credible interval were 0.69 (0.49, 0.95). The HR estimate and the large probability of a meaningful protective association, 95%, suggests a protective association of SGLT2i. Evidence was less robust for a meaningful protective association of SGLT2i for the CVD hospitalization outcome. In the matched cohort, there were 91 cardiovascular hospitalization events for SGLT2i and 123 for DPP4i (Table 3). The cause-specific hazard ratio estimate was 0.86 (0.65, 1.10) and the probability of a meaningful protective association was 67%. There were 48 cardiovascular death events for SGLT2i and 81 for DPP4i (Table 3). The cause-specific hazard ratio estimate was 0.70 (0.48, 0.99) and the probability of a meaningful protective association was 93%. A comparison of the estimated cause-specific hazard ratios and probabilities of meaningful protective associations for the MACE+HF outcome and its components can be found in Figure 3. The survival functions for the components of the MACE+HF are shown in Figure 2B–D.

Table 3 Event Rates and Cause-Specific Hazard Ratios for the Components of the Major Adverse Cardiovascular Events Plus Heart Failure

Figure 3 Hazard Ratio and Probability of Meaningful Association for MACE+HF and Components. The points are the estimated hazard ratios (HR) along with 95% credible intervals (horizonal lines) and the corresponding probability of a meaningful association (HR < 0.90) for the composite outcome of major adverse cardiovascular events plus heart failure (MACE+HF) and its components (HF Hospitalization, CVD Hospitalization, and CVD Death) found using the Bayesian time-to-event models with the Uninformative prior.

Discussion

This study demonstrates that a Bayesian survival model is a useful tool for incorporating prior information into a time-to-event analysis. Using a normal prior for the coefficient of interest in the hazard function and selecting the mean and variance of the prior based on the previous effect estimates provides a straightforward way to incorporate information from previous studies. Informative priors are a systematic way to incorporate skepticism or favorability for clinically meaningful effects. The time-to-event Bayesian model provides probability statements regarding the protective association of SGLT2i using the posterior distributions of the hazard ratios. The probability statements guard against a common misinterpretation in Cox proportional hazard regression analyses, that the inclusion of 1 in the 95% confidence interval for the hazard ratio proves a null association. The Bayesian approach allows one to make probability statements about clinically meaningful effects in cardiovascular studies as opposed to simply stating whether a statistically significant effect was found. The probability of a meaningful effect is a more useful statistic for clinical decision-making than a p-value or even a confidence interval.

The results from all three Bayesian models estimated a protective association of SGLT2i versus DPP4i and the MACE+HF outcome. The strength of this association was dependent on both the observed data and prior data incorporated using the informative priors. We also note that protective effect estimates were not as strongly supported for all the components of the composite MACE+HF outcome. The probability of a protective association for CVD Hospitalization was much lower than the probability of a protective association for the composite MACE+HF outcome and the other MACE+HF components, HF hospitalization and Cardiovascular death: 67% versus 93%, 95%, and 93%, respectively. This may have contributed to the observed mixed results for the associations between SGLT2i and cardiovascular outcomes in the literature.

Limitations should be noted. Like many observational studies, this study has the possibility for unobserved confounding. Multiple steps including propensity score matching and the addition of multiple data sources (VHA, Medicare and Medicaid) address the issue of controlling for possible confounding; however, residual confounding may still be present. Since Veterans may not receive all their care at VHA facilities, misclassification of those without CVD may have occurred and outcomes may have been missed. Because we considered episodes of care in the time-to-event analysis and not individual patients, Veterans could enter the study cohort more than one time and have correlated observations. Use of a weighted propensity scores model and robust standard errors can address this limitation. Since the primary objective of this study was to explore the implementation of a Bayesian survival model, we opted to use propensity score matching instead of weighting. This approach allowed us to separate the implementation of the Bayesian survival model from the propensity score adjustment. The study population was mostly white men with a mean age above 65 years who, based on our data, did not have a history of CVD. Therefore, the results may not be generalizable to populations with lower representation in VHA. This study did not evaluate patients who initiated use of DPP4i or SGLT2i as first-line therapy. It should be noted that most patients added the DPP4i or SGLT2i onto existing combination regimens (as a third agent).

Conclusion

This study demonstrates that a Bayesian survival model can be a useful tool to incorporate prior information into a time-to-event analysis and help guard against the misinterpretations often arising from looking solely at the 95% confidence intervals of the hazard ratios. The Bayesian approach allows one to make probability statements about clinically meaningful effects.