Comparative Effectiveness of SGLT2 Inhibitors, GLP-1 Receptor Agonists, DPP-4 Inhibitors, and Sulfonylureas on Risk of Kidney Outcomes: Emulation of a Target Trial Using Health Care Databases

Overall Study Design and Specification of a Target Trial

We followed the framework proposed by Hernán and Robins (23) to emulate a four-arm RCT (the target trial) of the comparative effectiveness of SGLT2i, GLP-1, DPP-4, or sulfonylureas on risk of kidney outcomes defined as a composite end point of eGFR decline >50%, ESKD, or all-cause mortality. First, the target trial protocol that would address the research question was specified (Supplementary Table 1). We then emulated the target trial by following the specified protocol to design the cohort, define exposure definitions and outcomes, select confounders, and specify the analytic approach (24,25). The protocol of the target trial, as well as key components of the target trial that was emulated, are specified in Supplementary Table 1.

Cohort Design

Participants were users of the VA Health Care System with a record of use of SGLT2i, GLP-1, DPP-4, or sulfonylureas after 1 October 2016 (N = 605,345), where date of first prescription was designated time 0 (T₀). Participants with no prescription record of the four antihyperglycemics within 1 year before T₀ met the inclusion criteria (n = 269,348). Participants were not included in the cohort if they had a history of type 1 diabetes, eGFR <15 mL/min/1.73 m², dialysis, or kidney transplant (selecting n = 243,289). Participants were further excluded if there was no recorded measurement of outpatient eGFR, hemoglobin A_1c (HbA_1c), height, weight, blood pressure, and LDL within the year before T₀, yielding an analytic cohort n = 216,558 (Supplementary Fig. 1). Participants were followed until the occurrence of an outcome or administrative end of follow-up (30 November 2019), where participant status was updated every 90 days.

Data Sources

We used VA databases, which have provided high-dimensional electronic health records since 1 October 1999 (26–29). The VA operates an integrated health care system consisting of >1,400 health care facilities, including 143 hospitals and 1,241 outpatient clinics. All veterans who are enrolled for receipt of care have access to medical benefits, including inpatient and outpatient services; preventive, primary, and specialty care; medications; mental health care; home, geriatric, and extended care; and medical equipment and prosthetics. The Corporate Data Warehouse (CDW) Outpatient Encounters domain and Inpatient Encounters domain were used to collect ICD-10 diagnosis codes, Current Procedural Terminology (CPT) codes, ICD-10 procedure codes, and type of health care system from the patients’ outpatient and inpatient encounters (30). Pharmacy data and provider information were obtained from the CDW Outpatient Pharmacy domain. The Managerial Cost Accounting data set and CDW Laboratory Results domain were used to collect information about laboratory data (31). Vital measurements were collected from the CDW Vital Signs domain, and demographic information was collected from the CDW Patient domain and VA Vital Status database (32).

Outcomes

The primary study outcome was time until the occurrence of the composite outcome of eGFR decline >50%, ESKD, or all-cause mortality. Time of ESKD was identified by the first record of dialysis, kidney transplant, or eGFR <15 mL/min/1.73 m² in inpatient or outpatient laboratory data sets, while date of death was obtained from the VA Vital Status database. eGFR values were estimated on the basis of the four-variable CKD Epidemiology Collaboration creatinine equation (33).

Exposure

Prescriptions of SGLT2i, GLP-1, DPP-4, or sulfonylureas were identified from outpatient pharmacy records; medications are listed in Supplementary Table 2. We defined two causal contrasts of interest: the intention-to-treat effect and the per-protocol effect. The intention-to-treat effect is the effect of being assigned to a treatment strategy at baseline, where treatment was defined by having a record of prescription of the assigned medication class at T₀. The per-protocol effect is the effect of continued use of the medication class assigned at T₀ throughout follow-up. Participants were considered to have discontinued use of the medication (nonadherent to the protocol) if a prescription was not refilled 90 days after the end of supply.

Covariates

We selected a set of predefined covariates, including predictors of antihyperglycemic treatment assignment on the basis of clinical knowledge and informed by prior RCTs. Predefined covariates were evaluated in the year before T₀ and included demographics, such as age, race (White, Black, and other), and sex; clinical measurements, such as HbA_1c, eGFR, albuminuria, blood pressure, LDL, and BMI; history of diseases, including congestive heart failure, cardiovascular diseases, cancer, alcoholism, hypoglycemia, diabetic ketoacidosis, acute kidney injury, bladder and urinary tract infections, venous thromboembolism, pancreatitis, and bone fracture; and history of medication use, including metformin, insulin, thiazolidinediones, other diabetes medications (α-glucosidase inhibitors and amylin analogs), statins, ACE inhibitors or angiotensin receptor blockers, β-blockers, diuretics, and calcium channel blockers. We also included other covariates, such as smoking status (never, former, current), number of outpatient visits, hospitalization length of stay, type of health care system where the antihyperglycemic was prescribed at T₀ (categorized as hospital system or outpatient clinic), and the calendar year of T₀. Medications, excluding insulin, were categorized into two levels on the basis of whether they were used for >90 days within 1 year before T₀. Insulin was categorized into a four-level variable on the basis of insulin use within 1 year before T₀: 1) never used; 2) used for a total of >90 days but not used within 90 days before T₀; 3) used within 90 days before T₀, but total duration of use was <90 days; and 4) total duration of use was >90 days and was used within 90 days before T₀. All laboratory measurements, besides albuminuria, were treated as continuous variables. Albuminuria was categorized as yes or no on the basis of an albumin-to-creatinine ratio >30 mg/g. Predefined continuous covariates, including age, eGFR, and HbA_1c, were controlled for as restricted cubic spline functions, where knots were set at the 5th, 30th, 60th, and 95th percentiles. Covariates were time updated for per-protocol analyses.

To more closely emulate the random treatment assignment in the target trial, we additionally leveraged the high-dimensional information available in VA electronic databases by using a high-dimensional variable selection algorithm to reduce potential imbalances in patient characteristics that influence treatment and are not captured by the predefined covariates (25,34). Variables within the year before T₀ were selected from outpatient records, inpatient records, pharmacy data, and laboratory data. Data were then organized into seven dimensions: outpatient ICD-10 diagnostic codes, outpatient CPT codes, inpatient ICD-10 diagnostic codes, inpatient CPT codes, and inpatient ICD-10 procedure codes for surgeries, pharmacy records, and laboratory results. For each dimension, the top 300 items (e.g., ICD codes, procedures, medications, laboratory tests) that occurred most frequently within cohort participants were selected. Abnormal laboratory values were categorized as above or below reference ranges. As a proxy of disease severity, for each participant, we evaluated the number of occurrences for each of the top 300 ∗ 7 dimensions = 2,100 items and then categorized each item’s occurrence into a set of three binary values: ever occurred, sometimes occurred (occurred at a frequency higher than that of 50% of cohort participants), and frequently occurred (occurred at a frequency higher than that of 75% of cohort participants), which generated, in total, 300 ∗ 7 ∗ 3 = 6,300 variables from high-dimensional data. In consideration of the positivity assumption, univariate relative risk of the association between these items and treatment assignment in all six treatment pairs were computed only for variables that occurred in >50 patients in all treatment arms (35). To ensure efficiency of probability estimation, we then selected the 500 high-dimensional variables with the strongest association with differences in treatment assignment (this was done independently for each treatment pair [SGLT2i vs. GLP-1, SGLT2i vs. DPP-4, SGLT2i vs. sulfonylureas, GLP-1 vs. DPP-4, GLP-1 vs. sulfonylureas, and DPP-4 vs. sulfonylureas] for a total of six comparisons). Among these selected variables, final selection of covariates was based on predictors of differences in treatment assignment between all pairwise comparisons of the four antihyperglycemics (which appeared in all six sets of pairwise comparisons). The final set of high-dimensional variables was used along with predefined variables to predict the probability of receiving each treatment. For estimation of the per-protocol effect, we used the same predefined predictors of treatment assignment as predictors for nonadherence and additionally applied a high-dimensional covariate selection algorithm at each time point. High-dimensional variables were selected on the basis of their relationship with treatment adherence (36). Both predefined and high-dimensional variables were time updated as appropriate and were used to predict the probability of following the protocol at each time point. Supplementary Fig. 2 provides an overview of the high-dimensional variable selection algorithm.

Statistical Analyses

For the overall cohort and across treatment arms, participant characteristics are presented as means and SDs for continuous variables and frequencies and proportions for categorical variables. Unadjusted survival probabilities across the four treatment arms were estimated by the Kaplan-Meier method.

A flowchart for the overall analytic approach for estimation of intention-to-treat and per-protocol effects is provided in Supplementary Fig. 2. Generalized propensity scores for the probability of treatment assignment were computed by a multinomial logistic regression model that predicted assignment to the four treatment groups where both the predefined variables and the high-dimensional variables were included (37,38). To take into consideration the potential effect modification of relationships with treatment assignment by past history and current use of insulin, we tested interaction terms between insulin and all predefined variables. These terms improved model fit (assessed through Bayesian information criterion) and so were retained.

To test the robustness of our emulation strategy, we applied negative and positive outcome controls to detect the presence of spurious associations and to test whether our approach would reproduce well-established associations (39). These analyses are designed to detect both suspected and unsuspected sources of spurious relationships (39). We used chronic lower respiratory diseases as a negative outcome control, where there is neither biologic plausibility nor a priori evidence suggesting the presence of a relationship (39). Changes in BMI (increase of >10% and decrease of >10% from baseline) were used as positive outcome controls, where prior evidence suggests that SGLT2i and GLP-1 are associated with weight loss and sulfonylureas are associated with weight gain (19). The associations were estimated through Cox survival models after adjusting for the probability of assignment to SGLT2i, GLP-1, and DPP-4 (the generalized propensity score). The intention-to-treat effects of medications on risk of the composite outcome of eGFR decline >50%, ESKD, or all-cause mortality was estimated through a Cox survival model after adjusting for the generalized propensity score (38). Hazard ratios (HRs) for comparisons between each treatment pair and adjusted survival curves are presented. Estimated number of events per 1,000 person-years and event rate difference were computed on the basis of the survival probability at 3 years obtained from Cox models. To evaluate intention-to-treat effects in different subgroups, subgroup analyses were conducted where subgroups were defined on the basis of baseline eGFR (≥90, <90 to ≥60, <60 to ≥45, and <45 mL/min/1.73 m²), baseline metformin use (use or non-use of metformin within 90 days before T₀), baseline cardiovascular disease (with and without cardiovascular diseases), BMI category (≤25, >25 to 30, and >30 kg/m²), and age at T₀ (≥65 and <65 years). For each subgroup analysis, high-dimensional variables were selected independently. We conducted multiple sensitivity analyses to examine the robustness of our results. We 1) estimated the treatment effect through inverse probability of treatment weighting with 0.1% truncation of stabilized weights (37,40); 2) estimated the risk of each individual component of the composite outcome, including eGFR decline >50%, ESKD, or all-cause mortality; 3) adjusted for time-varying HbA_1c with restricted cubic spline to investigate antihyperglycemic medication effects independent of glycemic control; 4) adjusted for time-varying BMI with restricted cubic spline to investigate effects independent of changes in weight; and 5) adjusted for both time-varying HbA_1c and BMI splines.

To offer more patient-centered and clinically relevant treatment effect estimates, we estimated the per-protocol effect of antihyperglycemics on the risk of the composite kidney outcome (41–43). The per-protocol effect in this study is the effect of continuous use of the assigned treatment throughout follow-up. To adjust for differences in probability of treatment assignment across treatment groups and probability of protocol adherence, we applied the approach of inverse weighting by treatment assignment and inverse weighting by the time-dependent probability of nonadherence to the protocol (43–45). Predictors of adherence (nondeviation from the protocol) included time-updated predefined variables as well time-updated high-dimensional variables, which were selected in the pseudocohort at each time point (k) among those who had adhered to the assigned treatment strategy at the prior time point (k − 1). The stabilized weight for adherence was estimated as follows:

where Z is an indicator of adherence; A is the treatment arm; V is a vector of time-independent predictors, including age, race, sex, type of health care system where the antihyperglycemic medication was initially prescribed, and year of T₀; and is the history of predefined variables and high-dimensional variables through time k − 1 (46). The time-dependent probability was estimated from logistic regressions at each time point. The time-dependent high-dimensional variables included the top 500 variables with the strongest univariate relative risk with adherence at time k, further restricted to those variables associated with differences in probability of assignment in all pairwise comparisons of antihyperglycemics (36). Weights were truncated at the 0.1th and 99.9th percentile before estimating the per-protocol effect. The per-protocol effect used pooled logistic regression, with follow-up time treated as a restricted cubic spline and knots placed at 180, 360, 540, 720, and 900 days. A 95% CI that does not cross unity was considered statistically significant. All analyses were done using SAS Enterprise Guide 7.1 software (SAS Institute, Cary, NC). The study was approved by the institutional review board of the VA St. Louis Health Care System.