Association of Cardiovascular Outcomes and Mortality With Sustained Long-Acting Insulin Only vs Long-Acting Plus Short-Acting Insulin Treatment

Key Points Question Among adults with type 2 diabetes receiving long-acting insulin with a glycated hemoglobin A1c level of 6.8% to 8.5%, is there a difference in cardiovascular events and mortality among individuals who do or do not initiate additional treatment with short-acting insulin? Findings In this cohort study of 57 278 adults, addition of short-acting insulin was associated with increased all-cause mortality compared with long-acting insulin alone and a decreased risk of acute myocardial infarction, with limited evidence of a difference in congestive heart failure. Meaning Given the lack of evidence demonstrating an increase in major cardiovascular events or cardiovascular mortality, the increased mortality with the combination of long- and short-acting insulin may be explained by noncardiovascular events or unmeasured confounding.


≥ eMethods 1 -Cohort Construction
We searched the electronic health records and administrative databases of the HP, KPCO, KPNC, KPSC health plans to identify all diabetes patients between 1/1/2000 and 12/31/2013 with a first insulin dispensing composed of only long-acting insulin between 1/1/2005 and 12/31/2013 who had a first elevated A1c between 6.8% and 8.5% 28 days or more after initiation of long-acting insulin. The algorithm used to identify diabetes patients is described below (second bullet). The elevated A1c date is referred to as the index date. Each patient who met all of the following criteria was included in the main study cohort: • age on date of first long-acting insulin dispensing and index date ≥21 and ≤89 • diabetes recognition occurred before or on first long-acting insulin dispensing where the diabetes recognition date was defined from the patient's diagnoses from inpatient, ambulatory, laboratory, and pharmacy encounters. Specifically, diabetes recognition was defined as the earlier of one inpatient diagnosis (ICD-9-CM 250.x, 357.2, 366.41, 362.01-362.07) or any combination of two of the following events occurring within a 24month period of time, using the date of the first event in the pair as the identification date: 1) A1C > 6.5% (48 mmol/mol); 2) fasting plasma glucose > 126 mg/dl (7.0 mmol/L); 3) random plasma glucose > 200 mg/dl (11.1 mmol/L); 4) an outpatient diagnosis code (same codes as inpatient); 5) any anti-hyperglycemic medication dispensed. For example, an individual with an A1C of 7.5% (57 mmol/mol) followed by an outpatient diagnosis of diabetes would be identified with diabetes on the (earlier) date of the A1C, with a laboratory result as the primary source. When the two events used for identification came from the same source (e.g., two outpatient diagnoses), they were required to occur on separate dates, but no more than 24-months apart. Note the following exception: two dispensings of metformin, thiazolidinediones, or liraglutide -with no other indication of diabetes -was not counted because these agents could be used for diabetes prevention, weight loss or to treat polycystic ovarian syndrome. Events that were identified during a pregnancy (within 270 days prior to a delivery) were excluded from consideration • the elevated A1c on the index date occurred while on long-acting insulin and no short-acting insulin was filled between the first long-acting dispensing and the index date • minimum of 12 months of health plan enrollment before first long-acting insulin dispensing and allowing for multiple gaps not exceeding 90 days combined. This criterion was also required for the 12 months before index date.
• minimum of 12 months of drug coverage before first long-acting insulin dispensing and allowing for multiple gaps not exceeding 90 days combined. This criterion was also required for the 12 months before index date.
• not pregnant on first long-acting insulin dispensing and on index date • no evidence of bariatric surgery in the 2 years before first long-acting insulin dispensing, i.e., no record of the following ICD-9 procedure and CPT- • no evidence of end stage renal disease in the 2 years before first long-acting insulin dispensing, i.e., no record of the following ICD-9 diagnosis, ICD-9 procedure, and CPT-4 codes (kidney transplant): v42.0, 996.81 ; 55.6, 55.61, 55.69 ; 50360, 50365, 50380 and most recent GFR laboratory result (if any) 15 and no record of 2 or more of the following ICD-9 diagnosis, ICD-9 procedure, and CPT-4 codes dated >90 days apart as primary or secondary diagnosis (dialysis): 585. 6 • no evidence of a stage 4 cancer diagnosis in the 2 years before first long-acting insulin dispensing, i.e., no record of the following ICD-9 diagnosis codes 197. x, 198.x, 199.x. This criterion was also required in the 2 years before index date.
− ≤ − ∈ • no evidence of hospice or palliative care in the 2 years before first long-acting insulin dispensing, i.e., no record of an hospice encounter and no record of the ICD-9 diagnosis code v66.7 and no record of the CPT code 99377 and 99378. This criterion was also required in the 2 years before index date.
• at least one A1c laboratory measurement recorded in the 2 years before first long-acting insulin dispensing. This criterion was also required in the 2 years before index date.
• insulins dispensed on first long-acting insulin dispensing do not include animal, inhaled, or short-acting insulins • diabetes of type 2 defined by the following ratio being strictly lower than 50%: the number of ICD-9 diagnosis codes 250.x1 and 250.x3 (type 1) in the 2 years before first long-acting insulin dispensing divided by the sum of this number and the number of ICD-9 codes 250.x0 and 250.x2 (type 2) in the 2 years before first long-acting insulin dispensing. If this ratio is not defined (i.e., denominator is 0), the diabetes type is unknown and the patient excluded from the study cohort. This criterion was also required in the 2 years before index date.
In addition to these criteria above, KPCO patients living outside the Denver/Boulder area were excluded due to incomplete data capture.

eMethods 2 -Data Structure and Notation
All analyses in this report are based on analytic datasets constructed with the MSMstructure SAS macro 1 to coarsen daily EHR data using the 90-day unit of time, i.e., time-dependent variables are updated every 90 days in the resulting analytic datasets. More specifically, for each of the five failure time outcomes considered (eTable 1), a separate analytic dataset is constructed by collecting the realizations of the random variables described below for all patients in the main or CVD study cohort.
Follow-up time (expressed in 90-day units) is denoted by t and, by convention, the first 90 days of follow-up are denoted by t = 0. The time when the patient's follow-up ends is denoted by T and is defined as the earliest of the time to failure denoted by T or the time to a right-censoring event denoted by C. When a patient is right-censored, i.e., C < T, the type of right-censoring event experienced by the patient is recorded and denoted by Γ with possible values 1-5 to represent the administrative end of study, disenrollment from the health plan, start of a pregnancy, initiation of a non standard insulin (i.e., inhaled or animal insulin), or death, respectively. The indicator that the end of follow-up is due to the occurrence of a failure event is denoted by ∆ = I(T ≤ C), i.e., ∆ = 1 implies that T = T and ∆ = 0 implies that T = C. Treatment with insulin therapy at time t is represented by the categorical variable A1(t) with four possible levels: long-acting insulin only (encoded by 0), both long-acting and short-acting insulin (encoded by 1), no insulin therapy (encoded by 2), short-acting insulin only (encoded by 3). The indicator of the patient's right-censored status at time t is denoted by A2(t). We thus have A2(t) = 0 for t = 0, . . . , T − 1 when T ≥ 1 and A2(T ) = 1 − ∆. The exposure variable denoted by A(t) is defined by A(t) = (A1(1), A2(t)). At each time point t = 0, . . . , T, covariates such as A1c measurements (eTables 2-3) are denoted by a component Lj(t) of the random vector L(t) and defined from measurements that occur before the exposure at time t, A(t), or are otherwise assumed not to be affected by the exposures at time t or thereafter, (A(t), A(t + 1), . . .). If no such measurements were collected, each variable Lj(t) is defined by convention using last observed value carried forward at t > 0. If no baseline measurements were collected for a continuous variable in L(0), the variable is defined by convention as the median of the baseline values from patients with observed measurements at t = 0. For categorical variables in L(0), a separate level is defined to encode missing baseline measurements. For each time-independent or time-dependent covariate Lj with at least one missing measurement (at baseline or at t > 0), an indicator of missing covariate measurement at time t is created and included as a distinct variable (e.g., to encode intensity of clinical monitoring) in the random vector L(t) for all time points t. In addition, the vector of covariates L(t) at time t include an outcome measurement denoted by Y(t), i.e., Y(t) L(t) for t = 0, . . . , T. For each time point t = 1, . . . , T + 1, the outcome is the indicator of past failure, i.e., Y(t) = I(T t 1) and Y(0) = 0 by convention. By definition, the outcome is thus 0 for t = 0, . . . , T , not observed at t = T + 1 if ∆ = 0 and, 1 at t = T + 1 if ∆ = 1.
In short, the observed data in each analytic dataset are realizations of n copies Oi of the random process O = T , ∆, (1 ∆)Γ, L (T ), Ā (T ), ∆Y(T + 1) where n = 57, 278 in each of the four analytic datasets to evaluate AMI, CHF, CVA, and all-cause mortality, and n = 39, 279 in the analytic dataset to evaluate CVD mortality. In the analyses of each dataset, we assumed 2 that the random variables Oi are independent and identically distributed and we denote their common distribution with P0. ≈ ≈ · · To simplify expressions below, we use the overbar notation ¯ to denote the history of a variable from baseline to time t (e.g., Ā (t) = (A(0), . . . , A(t))) and, by convention, L(t) and A(t) are nil when t < 0. Lower case notation is used to represent a possible realization of a random variable.

eMethods 3 -Causal Estimands and Inverse Probability Estimator
The causal estimands of interest are defined by two exposure regimens over the first 16 quarters of follow-up: 1) a static regimen denoted by g0 * and defined by continuous exposure to only long-acting insulin without occurrence of right-censoring events, and 2) a stochastic exposure regimen denoted by g1 * and defined by continuous exposure to long-acting and short-acting insulin without occurrence of right-censoring events where patients are given a 4-quarter grace period 3 after the index date to intensify insulin therapy through the addition of short-acting insulin. Formally, both interventions can be defined as stochastic interventions, i.e., the collection of conditional probabilities [4, Section 6], as follows.

Continuous insulin therapy with long-acting insulin and treatment intensification in the first 4 quarters of follow-up through the addition of short-acting insulin and continuous exposure to both long-and short-acting insulin thereafter:
• no right-censoring events • patient is free to initiate short-acting insulin at any time before the last quarter of the grace period with probabilities that approximate patterns of treatment intensification observed in the cohort where the probabilities g0 are defined by the distribution of the observed data process P0. These probabilities were estimated using the sample mean and were equal to 0.97 at quarters 0, 1, and 2. They represent the proportions of patients in the real data that remained treated with only long-acting insulin at quarter t among patients who were previously treated with only long-acting insulin and who were either treated with only longacting insulin or both long-acting or short-acting insulin at quarter t. With this intervention definition, 8.7% of patients would on average initiate short-acting insulin therapy at quarter 0, 1, or 2 in the arm of an ideal randomized experiment if patients perfectly complied with the assigned intervention g1 * . This can be contrasted with the rate of treatment intensification with short-acting insulin in the observed data used, for example, to examine AMI: 8.3% (4,729) of the n = 57, 278 patients initiated short-acting insulin therapy in quarter 0,1, or 2.
The following two working, 5 logistic, marginal structural models (MSMs) for the discrete-time counterfactual hazards defined by the prior two stochastic interventions, P(Yg * x (t + 1) = 1 Yg * x (t) = 0) with x = 0, 1, were considered: • a simple working MSM whose parameterization mimics a common modeling practice that assumes constant hazard ratios over time (i.e., a model based on the proportionality assumption): • a saturated working MSM whose parameterization permits hazard ratios to change over time: for t = 0, . . . , 15 where, for each working MSM, the collection of its coefficients is denoted by β. A standard 4, 6, 7 bounded and stabilized IPW estimator approach was implemented to fit each working MSM through a weighted logistic regression where each person-time-observation (A(t), Y(t + 1)) for t = 0, . . . , 15 was duplicated for each regimen x = 0, 1 and assigned the following inverse probability weight: with the following choice of stabilizing factor λ(x) = ∏ t Pn Fx(j) = 1 | Fx (j − 1) = 1) where each factor Pn denotes a sample mean and Fx(j) is defined as the indicator that the patient followed the intervention g * x at time j, i.e., Fx(j) = I g * x (A(j) | L (j), Ā (j − 1) > 0). The resulting IPW estimator of the working MSM coefficient β is denoted by βn and define the various effect measures reported below.

eMethods 4 -Denominator of the Inverse Probability Weights
The denominators of the IP weights (1) used to fit each of the two working MSM described above require estimation of the conditional probabilities g0( These probabilities can be factorized based on the following 11 propensity scores (PS) for: • continuation of only long-acting insulin therapy in the first quarter of follow-up among patients not right-censored in the first quarter (PS denoted by µ1(0)): • initiation of short-acting insulin therapy and continued use of long-acting insulin in the first quarter of followup among patients not right-censored in the first quarter and who interrupted treatment with only long-acting insulin during quarter 1 (PS denoted by µ2(0)): • continuation of only long-acting insulin therapy in any given quarter t after the first quarter of follow-up among patients not right-censored in quarter t and who were previously continuously exposed to only long-acting insulin (PS denoted by µ3(t)): • switching to only long-acting insulin therapy in any given quarter t after the first quarter of follow-up among patients not right-censored in quarter t and who previously initiated short-acting insulin within the first 4 quarters of follow-up and who remained continuously exposed to both long-acting and short-acting insulin thereafter through quarter t − 1 (PS denoted by µ4(t)): • continuation of both long-acting and short-acting insulin therapy in any given quarter t after the first quarter of follow-up among patients not right-censored in quarter t and who previously initiated short-acting insulin within the first 4 quarters of follow-up and who remained continuously exposed to both long-acting and short-acting insulin thereafter through quarter t 1 and who were not exposed to only long-acting insulin at quarter t (PS denoted by µ5(t)): • exposure to both long-acting and short-acting insulin in any given quarter t after the first quarter of follow-up among patients not right-censored in quarter t and who were previously continuously exposed to only longacting insulin except at quarter t (PS denoted by µ6(t)): • right-censoring due to administrative end of study at any given quarter t (PS denoted by µ7(t)): • right-censoring due to disenrollment from the health plan at any given quarter t among patients who were not right-censored at time t due to administrative end of study (PS denoted by µ8(t)): where I A2(t) = 1, Γ = j is the indicator that the patient experienced the j th type of right-censoring event • right-censoring due to start of pregnancy at any given quarter t among patients who were not right-censored at time t due to administrative end of study or health plan disenrollment (PS denoted by µ9(t)): where L (0) denotes the indicator that the patient is female and I A2(t) = 1, Γ 1, . . . , k = 0 is • right-censoring due to initiation of a non-standard insulins (animal or inhaled) at any given quarter t among patients who were not right-censored at time t due to administrative end of study, health plan disenrollment, or start of pregnancy (PS denoted by µ10(t)): • right-censoring due to death at any given quarter t among patients who were not right-censored at time t due to administrative end of study, health plan disenrollment, start of pregnancy, or initiation of non-standard insulins (PS denoted by µ11(t)): We note that the probability of a patient not experiencing a right-censoring event at quarter t given past covariates and exposures can then be derived from the last 5 PS as follows: For the AMI, CHF, CVA, and CVD mortality outcomes, we constructed the denominators of the IP weights (1) for all outcomes contributing to the MSM fits (i.e., when ∏ t g * x (A(j) | L (j), Ā (j − 1)) > 0) using the formulas below and estimates of the 11 PS above for t = 0, . . . , 15: • for replicates of person-time observations indexed by x = 0: where the product terms indexed by j are nil when t = 0.
• for replicates of person-time observations indexed by x = 1: where the product terms indexed by j are nil when t = 0.
For the all-cause mortality outcome, the same formulas for the denominators of the IP weights were used except that the factors involving the PS for death (i.e., (1 µ11(j))) were ignored (i.e., replaced by 1's) because death is then the failure outcome of interest (i.e., there is no right-censoring due to death). Each of the first three approaches considered for estimating the denominators of the IP weights using the formulas just described consists in fitting a separate logistic model for each of the the 11 PS µj(t) defined above. The three approaches only differ by the set of covariates that define each of the main terms included in each logistic model. We describe these sets in the next section.

eMethods 5 -Standard Propensity Score Estimation with Three Covariate Adjustment Sets
In the first approach implemented to estimate the denominators of the IP weights, the main terms included in a given PS logistic model were those associated with covariates presumed to impact both failure and the PS outcome as indicated in eTables 4-5. For instance, in the analyses of CHF, the PS logistic model for continuation of only long-acting insulin therapy in the first quarter of follow-up (µ1(0)) included main terms for all covariates in these tables where a value 1 is found in both the treatment and CHF columns. For the time-dependent covariates selected based on this rationale, only main terms for their current values L(t) were included in the PS logistic models, i.e., no main terms for other summary measures of the covariate histories were considered (e.g., latest change in value L(t) L(t 1) or a lagged value L(t 1)). All PS logistic models fitted with pooled data over time (i.e., µj(t) for j = 3, . . . , 11) also included main terms for time t (expressed in 90-day intervals). PS logistic models for right-censoring events (i.e., µj(t) for j = 7, . . . , 11) included two main terms indicating whether the patient followed one of the two interventions (g0 * or g1 * ) through quarter t 1. For the PS logistic models for administrative end of study (µ7(t)) and for the initiation of non-standard insulins (µ10(t)), only main terms for time t and the two indicators of prior exposures being consistent with interventions g 0 * and g 1 * were included in the model. For the PS logistic model for start of pregnancy (µ9(t)), only main terms for time t, age at index date, and the two indicators of prior exposures being consistent with interventions g0 * and g1 * were included in the model. All continuous variables considered by the various PS logistic models were discretized using the cutoffs given in eTable 6 and main terms for the resulting dummy variables (for the non-reference level) were included in the models. eTable 7 provides an example of the logistic model fit for µ1(0) based on the PS estimation approach 1.
The second approach implemented to estimate the denominators of the IP weights followed the same principles with the difference that the main terms included in a given PS logistic model (including for start of pregnancy (µ9(t)) and administrative end of study (µ7(t)) were those associated with covariates presumed to, at least, impact failure as indicated in eTables 4-5. However, for the PS logistic model for the initiation of non-standard insulins (µ10(t)), only main terms for t and the two indicators of prior exposures being consistent with interventions g0 * and g1 * were included in the model because <5 patients initiated non-standard insulins which limited the number of covariate that could be considered. All other modeling decisions were identical to those of the first approach described above. eTables 8-9 provide an example of the logistic model fit for µ1(0) based on the PS estimation approach 2.
The third approach implemented to estimate the denominators of the IP weights followed the same principles with the difference that the main terms included in a given PS logistic model were those associated with the covariates presumed to impact either failure or the PS outcome as indicated in eTables 4-5. The PS logistic models for treatment decisions (i.e., µj(t) for j = 1, . . . , 6) also included interaction terms between the study site indicators and index year indicators. The PS logistic models for the start of pregnancy (µ9(t)) and administrative end of study (µ7(t)) included main terms for all covariates presumed to affect failure. However, for the PS logistic model for the initiation of non-standard insulins (µ10(t)), only main terms for t and the two indicators of prior exposures being consistent with interventions g0 * and g1 * were included in the model because <5 patients initiated non-standard insulins which limited the number of covariate that could be considered. All other modeling decisions were identical to those of the first approach described above. eTables 10-12 provide an example of the logistic model fit for µ1(0) based on the PS estimation approach 3.
Thus, the three sets of variables that define the main terms included in any given PS logistic model according to the three approaches just described are nested and of increasing size.

eMethods 6 -Data-adaptive Propensity Score Estimation
In the fourth approach implemented to estimate the denominators of the IP weights, a separate super learner 8 was used to estimate each of the 10 PS µj(t) with j = 1, . . . , 9, 11 instead of a separate logistic model (as done in the first three approaches). Because <5 patients initiated non-standard insulins, the same logistic model for estimating µ10(t) (initiation of non-standard insulins) as the one used in the prior three approaches was also used in approach 4. Each super learner was constructed based on 10-fold cross-validation and the following 15 learners: • the same logistic model as in approach 1 • the same logistic model as in approach 2 • the same logistic model as in approach 3 • a logistic model constructed using the same principles described for approach 2 but with the difference that only main terms for covariates presumed to at least impact failure as indicated by a rank or 1 in eTables 4-5 • a logistic model constructed using the same principles described for approach 2 but with the difference that only main terms for covariates presumed to at least impact failure as indicated by a rank or 1 or 2 in eTables 4-5 • a logistic model constructed by including main terms for the first 30 covariates most associated with the PS outcome, i.e., the 30 covariates with the smallest p-values for the test that the Pearson's product moment correlation coefficient is equal to 0 (implemented by the screen.corRank screener in the SuperLearner R package ? ) • a logistic model constructed by including main terms for the first 20 covariates most associated with the PS outcome (selected by the screen.corRank screener) • a logistic model constructed by including main terms for the first 10 covariates most associated with the PS outcome (selected by the screen.corRank screener) • a logistic model constructed by including main terms for the first 5 covariates most associated with the PS outcome (selected by the screen.corRank screener) • a regression based on linear splines and their tensor products that considers only the first 10 covariates most associated with the PS outcome (selected by the screen.corRank screener). The regression is implemented by the SL.polyclass routine given below that implements the polyclass learner 9 based on the Bayesian Information Criterion (BIC) as the model selection criterion. To improve computing speed, this learner was favored over the SL.polymars routine that is available by default in the SuperLearner R package but that relies on cross-validation for model selection. • a regression based on linear splines and their tensor products that considers only the first 5 covariates most associated with the PS outcome (selected by the screen.corRank screener). The regression is implemented by the SL.polyclass routine given above that implements the polyclass learner 9 based on the Bayesian Information Criterion (BIC) as the model selection criterion.
• a random forest regression (implemented by the SL.randForest routine) that considers only the first 10 covariates most associated with the PS outcome (selected by the screen.corRank screener) • a random forest regression (implemented by the SL.randForest routine) that considers only the first 5 covariates most associated with the PS outcome (selected by the screen.corRank screener) • an extreme gradient boosting regression (implemented by the SL.xgboost routine) that considers only the first 10 covariates most associated with the PS outcome (selected by the screen.corRank screener) • an extreme gradient boosting regression (implemented by the SL.xgboost routine) that considers only the first 5 covariates most associated with the PS outcome (selected by the screen.corRank screener) eTable 13 provides an example of the super learner fit for µ1(0) based on the PS estimation approach 4.

eMethods 7 -Results
eTable 14 describes the distribution of reasons for end of follow-up in all primary analyses. The counts of patients following each exposure regimen at each quarter of follow-up in the AMI and CVD mortality analyses are described by the histograms in eFigure 1. These counts were similar for the CHF, CVA, and all-cause mortality analyses (data not shown). We note that it is possible for the same patient to follow both exposure regimens simultaneously during the first 3 quarters of follow-up only.