Marginal structural models for the estimation of the risk of Diabetes Mellitus in the presence of elevated depressive symptoms and antidepressant medication use in the Women’s Health Initiative observational and clinical trial cohorts

Background We evaluate the combined effect of the presence of elevated depressive symptoms and antidepressant medication use with respect to risk of type 2 diabetes among approximately 120,000 women enrolled in the Women’s Health Initiative (WHI), and compare several different statistical models appropriate for causal inference in non-randomized settings. Methods Data were analyzed for 52,326 women in the Women’s Health Initiative Clinical Trials (CT) Cohort and 68,169 women in the Observational Study (OS) Cohort after exclusions. We included follow-up to 2005, resulting in a median duration of 7.6 years of follow up after enrollment. Results from three multivariable Cox models were compared to those from marginal structural models that included time varying measures of antidepressant medication use, presence of elevated depressive symptoms and BMI, while adjusting for potential confounders including age, ethnicity, education, minutes of recreational physical activity per week, total energy intake, hormone therapy use, family history of diabetes and smoking status. Results Our results are consistent with previous studies examining the relationship of antidepressant medication use and risk of type 2 diabetes. All models showed a significant increase in diabetes risk for those taking antidepressants. The Cox Proportional Hazards models using baseline covariates showed the lowest increase in risk , with hazard ratios of 1.19 (95 % CI 1.06 – 1.35) and 1.14 (95 % CI 1.01 – 1.30) in the OS and CT, respectively. Hazard ratios from marginal structural models comparing antidepressant users to non-users were 1.35 (95 % CI 1.21 – 1.51) and 1.27 (95 % CI 1.13 – 1.43) in the WHI OS and CT, respectively – however, differences among estimates from traditional Cox models and marginal structural models were not statistically significant in both cohorts. One explanation suggests that time-dependent confounding was not a substantial factor in these data, however other explanations exist. Unadjusted Cox Proportional Hazards models showed that women with elevated depressive symptoms had a significant increase in diabetes risk that remained after adjustment for confounders. However, this association missed the threshold for statistical significance in propensity score adjusted and marginal structural models. Conclusions Results from the multiple approaches provide further evidence of an increase in risk of type 2 diabetes for those on antidepressants.


Background
Diabetes is a chronic illness with serious health consequences, such as adult blindness, non-traumatic limb amputation, renal failure and neuropathy. Previous literature has noted considerable diabetes and depression among postmenopausal women, with a prevalence rate that is approximately 12 % for each [1,2]. Recent literature also suggests an increased risk of diabetes among those who are depressed and on antidepressant medications [3 -6 ]. It is increasingly important to further investigate whether depression or antidepressant medication is influencing this association given that approximately 11 % of American women take antidepressant medication, and use is rising [7]. While the rates of use for depression treatment has remained the same, off-label use of antidepressants has increased significantly [8]. Examples of off-label use include treatment for certain types of pain, fibromyalgia, insomnia, and general unhappiness. In this analysis, we compare four statistical approaches to evaluate the combined effect of the presence of elevated depressive symptoms and antidepressant medication use on incident type 2 diabetes using data on approximately 120,000 women in the Women's Health Initiative (WHI).
The WHI is a longitudinal study including repeated measurements for presence of elevated depressive symptoms, antidepressant use and self-reported diagnosis of type 2 diabetes. A previously reported analysis in the WHI [9] found that elevated depressive symptoms and antidepressant medication use resulted in increased type 2 diabetes risk. This analysis was based on Cox models that are subject to bias in the presence of timedependent confounding [10][11][12]. A probable mechanism that gives rise to time-dependent confounding in this context is shown in Fig. 2 (Appendix) -participants with increased BMI could be more likely to be depressed and/or on antidepressants in the future [13]. Moreover, BMI and depression/antidepressant medication use can also significantly influence future diabetes risk. Thus, an analysis based on Cox models as in the previous literature [9] to evaluate the causal relationship between presence of elevated depressive symptoms/antidepressant use (exposure) and diabetes risk (outcome) could result in biased estimates by misattributing the effects of time-dependent confounders to the exposures of interest.
In the presence of time-dependent confounding, all Cox models including those that incorporate timevarying covariates are subject to bias [10][11][12]. Marginal structural models (MSMs) overcome limitations of Cox models -MSMs use inverse probability of treatment/ exposure weighting and, given assumptions, can yield consistent and unbiased causal estimates of the effect of exposure [10][11][12]. Other methods in this setting include propensity score adjusted models in which the probability of treatment assignment conditional on baseline covariates (or propensity score) is adjusted for as an additional covariate.
The objective of this work was to estimate the combined effect of antidepressant medication use/presence of elevated depressive symptoms on type 2 diabetes in the WHI, and to compare results obtained from MSMs and propensity score adjusted Cox models to more traditional approaches such Cox models [9,14,15]. The results presented in this paper go beyond a previously reported analysis of this hypothesis in the WHI [9] by comparing results from four different statistical approaches including MSMs that are recommended for causal inference in observational studies.

Women's health initiative (WHI)
The WHI enrolled 68,132 participants into clinical trials (WHI-CT) and 93,676 participants into an observational study (WHI-OS) between 1993 and 1998 [16][17][18][19]. Eligibility criteria included: postmenopausal women aged 50 to 79 years, reliable/mentally competent, and expected survival and local residency for at least three years. Medication use, presence of elevated depressive symptoms, and diabetes status were collected from participants over an average of 7.6 years of follow-up. We analyzed data with follow-up to 2005.

Study variables Incident diabetes
Diabetes status was assessed by self-report at baseline and at each annual follow-up visit, which has been found to be a reliable indicator of diagnosed diabetes [20]. Time to diabetes was calculated as the interval between study enrollment and development of diabetes as evidenced by an annual medical history update, or censorship (death or end of study participation).

Antidepressant medication use
WHI-CT participants were instructed to bring all current prescription and nonprescription medications in original containers to clinic visits at baseline and years 1, 3, 6, and 9. WHI-OS medication data were collected at baseline and year 3. The Master Drug Data Base (MDDB: Medi-Span, Indianapolis, IN) was used to categorize the medications. Based on the MDDB classification, a binary indicator for antidepressant medication use was created. Antidepressants include the following major groups: 1) Selective serotonin reuptake inhibitors; 2) Monoamine oxidase inhibitors; 3) Tricyclic antidepressants; 4) Tetracyclics; 5) Serotonin/norepinephrine reuptake inhibitors (SNRIs); 6) Aminoketones; 7) Triazolopyridines; and 8) Dibenzoxazepine. A dichotomous indicator of antidepressant medication use was then created for each measurement period [9]. Due to sample size limitations, we did not perform analyses by class of medication.

Elevated depressive symptoms
Elevated depressive symptoms were measured using the 6-item Center for Epidemiological Studies Depression Scale (CES-D) [21]. A participant was determined to have elevated depressive symptoms if their score was 5 or higher on the CES-D. Presence of elevated depressive symptoms was available at baseline and year 3 in the WHI-OS and at baseline, year 1, and close out in the WHI-CT. Only 6 % of participants were assessed at years 3, 6, and 9. Due to high levels of missing values at year 1 and later, analysis in the CT cohort adjusted only for baseline presence of elevated depressive symptoms.

BMI
BMI was available at baseline and year 3 in the WHI-OS and at baseline, years 1, 3, 6 and 9 in the WHI-CT.

Other covariates
Other variables available at baseline for inclusion in multivariable models include: age, race/ethnicity (American Indian/Alaskan Native; Asian or Pacific Islander; Black/ African American; Hispanic/Latino; White; Other), education (<=high school; high school or GED; > = high school, but less than 4 years of college; 4 or more years of college), minutes of recreational physical activity per week, total energy intake, hormone therapy use (never, former, current), family history of diabetes (no, yes, don't know) and smoking status (never smoked, past smoker, current smoker).

Statistical analysis Analysis datasets
WHI OS: included 68,169 women after exclusions for self-reported diagnosis of diabetes at baseline or missing data on one or more of the following: baseline diabetes status, race/ethnicity, presence of elevated depressive symptoms, antidepressant medication use, BMI or women for whom the Year 3 visit occurred more than 3.5 years post enrollment (Fig. 1). 3624 (5.3 %) women reported diagnosis of type 2 diabetes during follow-up. Antidepressant use, BMI, presence of elevated depressive symptoms was available at baseline and year 3 and included in analyses.
WHI CT: included 52,326 women after the following exclusions -women with a self-reported diagnosis of diabetes at baseline or missing data on one or more of the following: diabetes status at baseline, race/ethnicity, baseline presence of elevated depressive symptoms, antidepressant medication use, BMI (Fig. 1). 4171 (8.0 %) women reported diagnosis of type 2 diabetes during follow-up. Antidepressant use and BMI were recorded at baseline, year 1 and year 3. Models adjusted for presence of elevated depressive symptoms at baseline because presence of elevated depressive symptoms was only measured on a small percentage (6.26 %) of women after year 1.

Statistical models
We compared results from four approaches for estimating the association between the presence of elevated depressive symptoms/antidepressant medication use with incident type 2 diabetes, in the WHI-OS and WHI-CT datasets separately. Results reported by Ma and colleagues [9] were based on multivariable Cox models as in Approach 1 and Approach 2 described below. The results from these models were compared to propensity score adjusted models (Approach 3) and MSMs (Approach 4) [22][23][24][25]. The interaction of elevated depressive symptoms and antidepressant medication use was investigated in all models, but all interactions were insignificant (p > 0.26 in WHI-OS, p > 0.06 in WHI-CT).
Approach 1 is a Cox proportional hazards (PH) model including the following covariates at baseline: elevated depressive symptoms, antidepressant use, BMI, age, ethnicity, education, minutes of recreational physical activity per week, total energy intake, hormone therapy use, family history of diabetes and smoking status.
Approach 2 is a Cox model with time-varying values of elevated depressive symptoms (WHI-OS only), antidepressant medication use and BMI, and baseline values of age, ethnicity, education, minutes of recreational physical activity per week, total energy intake, hormone therapy use, family history of diabetes and smoking status.
Approach 3 is similar to Approach 2, with additional adjustment for propensity for taking antidepressants at baseline. The propensity score was calculated to predict antidepressant use at baseline (outcome) from a logistic model. Predictors included in the model were baseline measures of BMI, age, ethnicity, education, minutes of recreational physical activity per week, total energy intake, hormone therapy use, family history of diabetes and smoking status.
Approach 4 is a MSM for a time to event outcome [10]. Two separate logistic regression models were fit to calculate the probability of being on antidepressants and the probability of censoring (not observing the outcome by that point in time). The models incorporated timevarying data on the presence of elevated depressive symptoms and BMI while adjusting for baseline values of age, ethnicity, education, minutes of recreational physical activity per week, total energy intake, hormone therapy use, family history of diabetes and smoking status. These models were used to determine the stabilized Inverse Probability of Treatment (IPTW) weights. In the WHI-CT analysis, models were also adjusted for clinical site. To relax the linearity assumption, we added a quadratic function of time to the model. Model details and associated SAS code are included in the Appendix. Table 1 presents baseline characteristics of our study population, by antidepressant medication use. Baseline characteristics by presence of elevated depressive symptoms are shown in Table 4 (Appendix). In both cohorts, the mean age was approximately 63 years old (63.6 (OS); 62.8 (CT)). Women were primarily White (86.

OS cohort Elevated depressive symptoms
An unadjusted model including only baseline values of antidepressant medication use and presence of elevated depressive symptoms showed a significant increase in diabetes risk for those with elevated depressive symptoms (HR 1.34; 95 % CI: 1.23-1.45). Approaches 1 and 2 that adjust for other confounders also resulted in a statistically significant association between the presence of elevated depressive symptoms and increased risk of diabeteshowever, this association did not remain statistically significant in propensity score adjusted models and MSMs.

Antidepressant medication use
All models showed a consistent, statistically significant increase in diabetes risk for those exposed to antidepressant medications vs. those who were not ( ]. There were no significant differences in the HRs estimated for the presence of elevated depressive symptoms between Approaches 2 and 3. The results from MSMs (Approach 4) were almost identical to Approaches 2 and 3 -the HR (CI) for antidepressant medication use was 1.35 (95 % CI 1.21 -1.51). In this application, the MSM approach yielded similar results to the traditional extended Cox model.

CT cohort Elevated depressive symptoms
As in the WHI-OS, the unadjusted model showed a significant increase in diabetes risk for those with elevated depressive symptoms. While this association was statistically significant in Approaches 1 and 2, it did not remain so in the propensity score adjusted model (Approach 3) and the MSM (Approach 4).

Antidepressant medication use
As in the WHI OS, all models in the WHI-CT showed a consistent, statistically significant increase in diabetes risk for those exposed to antidepressant medications vs. those who were not. Table 3 presents an estimate of variation in antidepressant medication use and presence of elevated depressive symptoms at baseline and year 3, by cohort. In the WHI-OS, 4.9 % of women were using antidepressant medication at both time-points, whereas 88.5 % never used them. 2.2 % of women who were using antidepressant medication at baseline had stopped by year 3, and 4.5 % who were not using antidepressant medication at baseline had started by year 3. 6.4 % experienced elevated depressive symptoms at both baseline and year 3, while 76.3 % never did. 8 % of women who experienced elevated depressive symptoms at baseline did not report experiencing those symptoms at year 3, and 9.3 % without elevated depressive symptoms at baseline did experience them by year 3. Similar patterns were observed in the WHI-CT. Table 5 (Appendix) presents the HRs and 95 % CIs for all covariates in the MSMs, in the WHI-OS and CT cohorts. Table 6 (Appendix) presents the distributions of the IPTW weights, including the estimated probability of having one's own observed treatment history and censoring history at follow-up time points. The probability of remaining uncensored was close to 1 for both cohorts at each follow-up time point given both the baseline and time-varying covariates. There was variation in the probability of having one's own observed treatment history, but the mean and median were close to 1 at 36 month follow-up in the WHI-OS and 12 and 36 month follow-up in the WHI-CT.

Discussion
Previous research has shown that the prevalence of elevated depressive symptoms and diabetes is high in postmenopausal women [1,26]. Ma and colleagues [9] found an increased risk of type 2 diabetes among women in the WHI cohorts who reported elevated depressive symptoms  Of the four different modeling approaches considered in this paper, MSMs are the gold standard for use in observational studies in which the presence of time-dependent confounding (Fig. 2, Appendix) is a possibility. For observational study settings in the absence of time-dependent confounding, propensity score adjusted models can be used to adjust for bias in exposure assignment. Propensity score adjusted models correct for confounding by indicationthis bias is present in studies in which individuals who are prescribed or take a medication are inherently different in their risk profile with respect to outcome when compared to those who do not take the drug. In the absence of confounding by indication and time-dependent confounding, simpler Cox models with or without timevarying covariates are appropriate.    We did not observe differences in hazard ratio estimates between the four modeling approaches that were considered. One explanation for the concordance of estimates from these different approaches suggests that time-dependent confounding by BMI was not a substantial factor in these datathus, BMI measured over the course of observation may not be strongly affected by exposure (i.e., as an intermediate) and/or did not exert a strong influence on our exposure and outcome of interest (i.e., as a confounder). However, other potential explanations for this observed concordance of results include limited longitudinal measurements of the key exposure variables, measurement error of the confounder and/or incorrect model specification of the dose response relationship with respect to the effects of the confounder on both exposure and outcome. However, this study was not designed to pinpoint the specific factor that caused the concordance of estimates from the various models.
A limitation of this work is that we had limited longitudinal follow up. The WHI-OS had two repeated measures, at baseline and at year 3. The WHI-CT had more repeated measures available, but we were only able to utilize three time points (baseline, year 1, year 3) due to high levels of missing data at later time points. In addition, in the WHI-CT, because presence of elevated depressive symptoms was measured on only a small percentage of participants after year 1, our models could not incorporate this factor as a time varying exposure. Due to limitations of the available data, the analysis did not account for antidepressant dose and adherence. While a bidirectional association between depression and diabetes risk is biologically plausible, our study was not designed to tease out the direction of association. Lastly, due to cost considerations, diabetes status was ascertained through self-reported questionnairesthis could result in modest levels of outcome misclassification. Statistical models that account for the error-prone self-reported outcomes would be useful in this context.

Conclusions
Our analyses provide further evidence that a significant increase in diabetes risk is observed for those on antidepressant medications in the WHI. In addition, our results comparing modeling approaches demonstrates that in some settings, results from more complex methods such as MSMs may not differ substantially from traditional methods of analysishowever, we recommend that these methods be explored to establish the validity of initial findings.

SAS Code to Fit the Marginal Structural Cox Proportional Hazards Model
Here we provide details on how to set up the data, along with SAS code to perform the analyses in Approach 4 using Marginal Structural Models for a time to event outcome. The first part of the analysis uses the proc logistic procedure to calculate the probability of not being on antidepressants (Models 1 and 2). Data in the file 'iptw' contains 1 record per subject for each time point included in the analysis. For example, in the WHI-OS dataset, we have measurements at baseline and year 3. Each subject has two records, one for baseline and one for year 3. Next we use the proc logistic procedure to calculate the probability of not being censored (Models 3 and 4). In the data file 'iptw2' the data has been expanded so that a subject has 1 record per month until her time to diabetes or time to censoring. The dataset 'main' merges Models 1-4 together to calculate the IPTW weights. This again is an expanded dataset so a subject has 1 record per month until her time to diabetes or time to censoring. The dataset 'main_w' contains the inverse probability of treatment weights and is the dataset used to run the marginal structural model analysis. Last, we use the proc genmod procedure to fit the final pooled logistic regression model to obtain the estimates of our association of interest.
Antidep is a dichotomous 0/1 indicator of the participant being on or off antidepressants. Model 1 includes that dichotomous indicator as the outcome, a timedependent intercept and baseline values of the following covariates as predictors: presence of elevated depressive symptoms (base_depression), BMI (base_bmi), Age (age), Ethnicity (ethnic), Education (educ4), Minutes of physical activity per week (tminwk), Total energy intake (base_energy_cat), Family history of diabetes (diabrel), Hormone therapy use (hormstat) and smoking status (smoking).    medication use (antidepressant_tv), presence of elevated depressive symptoms (depression) and BMI (bmi). Models 3 and 4 also include variables for 'month' and 'month 2 ' to relax the linearity assumption. We merge Models 1-4 together and in the following data step use the predicted values from those models to compute our IPTW estimates, denoted by stabw and nstabw for the stabilized and non-stabilized weights, respectively.
In our last step, we use the proc genmod procedure to fit the weighted pooled logistic model to obtain estimates of our association of interest from our Marginal Structural Cox PH Model. The outcome here is diabetes, which is a dichotomous indicator of whether or not the participant developed diabetes during that month. The patient ID variable and the independent working correlation matrix (subject = id/type = ind) must be specified. We weighted the model used the stabilized weights by using the 'scwgt stabw' statement in the procedure. The 'estimate' statement asks the procedure to report the odds ratio for our main associations of interest, in addition to the coefficients in the model.