Study population

The Nurses’ Health Study, a prospective cohort study initiated in 1976, enrolled 121 700 female registered nurses aged 30 to 55 years.16 The Health Professionals Follow-Up Study cohort was established in 1986 and enrolled 51 529 male health professionals aged 40 to 75 years.17 In both studies, detailed information on dietary and lifestyle factors, medical history, and disease status was collected at baseline and updated every two to four years through validated questionnaires.18 The cumulative response rate exceeded 90% for both cohorts. In the current analysis, we included participants with prevalent type 2 diabetes at baseline (1980 for the Nurses’ Health Study, and 1986 for the Health Professionals Follow-Up Study, when dietary information was first collected using a validated food frequency questionnaire), as well as participants with a diagnosis of incident type 2 diabetes during follow-up to 2018. We excluded participants if they had type 1 diabetes, CVD, or cancer at baseline; reported CVD or cancer before the diagnosis of type 2 diabetes during follow-up; left more than nine blank responses on the food frequency questionnaire or reported implausible daily energy intakes (<2510 or >14 644 kJ/day for women, and <3347 or >17 573 kJ/day for men); or had incomplete information on beverage consumption or dietary data at diabetes diagnosis. After exclusions, a total of 11 399 participants of the Nurses’ Health Study and 4087 participants of the Health Professionals Follow-Up Study with type 2 diabetes were included in the current analysis. For the analysis of changes in beverage consumption from before to after the diabetes diagnosis, we further excluded participants with type 2 diabetes at baseline or those with missing data on beverage consumption assessed before the diabetes diagnosis (n=2715), which left 9252 women and 3519 men for the change analysis.

Assessment of beverage intake

Intake of beverages was assessed using the validated food frequency questionnaires administered every two to four years. Participants were asked how often, on average (never to >6 times per day), they had consumed SSBs, ASBs, fruit juice, coffee, tea, low fat milk (skimmed, 1% or 2% fat), full fat milk, or plain water of a prespecified portion size (cup, glass, can, or bottle). SSBs included caffeinated colas, caffeine-free colas, other carbonated SSBs, and non-carbonated SSBs (fruit punches, lemonades, or other fruit drinks). ASBs included low calorie cola with caffeine, low calorie caffeine-free cola, and other low calorie beverages. Fruit juices included orange, apple, grapefruit, or other fruit juices. Coffee included caffeinated and decaffeinated varieties. A validation study conducted among a subsample of the participants in the Health Professionals Follow-Up Study showed reasonable validity for the assessment of beverage intake.19 The correlation coefficients between the food frequency questionnaire and multiple diet records were 0.84 for colas, 0.73 for low calorie colas, 0.75-0.89 for fruit juices, 0.93 for coffee, 0.77 for tea, 0.88 for low fat milk, 0.67 for full fat milk, and 0.52 for plain water.19 Similar correlation coefficients were found in a validation study conducted among a subsample of participants in the Nurses’ Health Study.20 Our primary dietary factors of interest were specific types of beverage consumption assessed after the diabetes diagnosis, and changes in beverage consumption before and after the diagnosis. We assessed the pre-diabetes beverage intake from the most proximal questionnaires before diabetes was ascertained.

Ascertainment of type 2 diabetes

Participants who reported a physician diagnosis of diabetes mellitus in the biennial questionnaires were sent a validated supplementary questionnaire about diagnostic tests, symptoms, and hypoglycemic treatment. The National Diabetes Data Group and American Diabetes Association criteria were applied to ascertain a diagnosis of type 2 diabetes (see supplementary appendix).2122 We excluded from the current analysis those participants who reported a diagnosis of type 1 diabetes on the supplementary questionnaire. Studies among 62 participants in the Nurses’ Health Study and 59 participants in the Health Professionals Follow-Up Study showed high validity in the supplementary questionnaire, with 98% and 97% of questionnaire confirmed type 2 diabetes diagnoses validated by medical record review in these women and men, respectively.2324

Ascertainment of outcomes

The primary endpoint was all cause mortality. We also examined the secondary outcomes of CVD incidence and mortality. Deaths were identified from reports by the next of kin or postal authorities or from searches of the National Death Index (see supplementary appendix).25 ICD-9 (international classification of diseases, ninth revision) codes were used to classify deaths from CVD (codes 390-459), cancer (codes 140-208.32), or other causes. Incident CVD was defined as fatal and non-fatal coronary heart disease, including coronary artery bypass graft surgery and non-fatal myocardial infarction, and as fatal and non-fatal stroke (see supplementary appendix).

Assessment of covariates

In both cohorts, information on lifestyle factors and medical history was collected at baseline and in biennial questionnaires. The supplementary appendix provides details of the assessments of covariates. To assess overall diet quality, we calculated the Alternate Healthy Eating Index (AHEI) score based on intakes of 11 foods and nutrients predictive of chronic disease risk, including vegetables, fruits, whole grains, SSBs and fruit juice, nuts and legumes, red and processed meat, trans fatty acids, long chain omega 3 fatty acids, other polyunsaturated fats, sodium, and alcohol.26 In the current study, we modified the AHEI score by excluding the consumption of SSBs and fruit juices.

Statistical analysis

The Kolmogorov-Smirnov normality test was used to assess distributions of continuous variables for normality, and natural logarithm transformations of skewed variables were applied before analyses. In descriptive analyses, continuous variables were expressed as means (standard deviations) for normally distributed variables or medians (interquartile ranges) for skewed variables, and categorical variables were represented by frequency and percentage. General linear models were used to calculate mean characteristics of the study participants at the time of diabetes diagnosis, and a test for linear trend using the Wald test was performed by assigning the median value to each category of beverage consumption and modeling this variable as a continuous variable.

For each participant, we calculated person years of follow-up from the date of diabetes diagnosis to the date of occurrence of study outcomes, last return of a valid follow-up questionnaire, or end of follow-up (30 June 2018 for the Nurses’ Health Study, and 30 January 2018 for the Health Professionals Follow-Up Study), whichever came first. Because changes in diet after a diagnosis of cancer could distort the associations of interest, for the CVD incidence analyses we stopped updating dietary variables after participants reported a diagnosis of cancer. For mortality analyses, dietary intake was not updated after a diagnosis of cancer or CVD. Time varying Cox proportional hazards models, conditioned on age and follow-up cycle, were applied to estimate hazard ratios and 95% confidence intervals for the associations of each beverage intake with all cause mortality, CVD incidence, and CVD mortality. Changes in beverage intake from before to after the diabetes diagnosis were defined as the absolute difference in beverage consumption (time varying post-diabetes beverage intake minus pre-diabetes beverage intake). The time varying covariates assessed during follow-up were considered in the multivariable models. Missing data for beverage consumptions and covariates during follow-up were replaced by the most recent valid assessments. In the multivariable model, we adjusted for age (years), duration of diabetes (years), sex (men or women), white ethnicity (yes or no), physical activity (<3.0, 3.0-8.9, 9.0-17.9, 18.0-26.9, ≥27.0 metabolic equivalents of task-hours/week), smoking status (never, former, current 1-14 cigarettes/day, current ≥15 cigarettes/day), alcohol consumption (0, 0.1-4.9, 5.0-14.9, ≥15.0 g/day), menopausal status and post-menopausal hormone use (pre-menopause, post-menopause (never, former, or current hormone use), or missing; Nurses’ Health Study only), family history of type 2 diabetes (yes or no) or myocardial infarction (yes or no), intake of total energy, and modified AHEI score (all in fourths). To further reduce the impact of confounding by existing comorbidities, disease management, and weight change, we further included history of hypertension (yes or no) or hypercholesterolemia (yes or no), use of antihypertensive (yes or no) or lipid lowering drug (yes or no), aspirin use (yes or no), diabetes drug use (oral drug only, insulin use, or others), and change in body mass index (BMI) before to after the diabetes diagnosis in the fully adjusted model.27 We mutually adjusted for different types of beverage intakes in the analysis of specific types of beverages. To obtain overall estimates for men and women and to increase statistical power, we pooled the hazard ratios from each model from the two cohorts with the use of an inverse variance weighted meta-analysis by the random effects model, which accounted for between study heterogeneity.28 CVD incidence and mortality were also examined according to the per serving intake of beverages. In the analysis of changes in beverage consumption from before to after the diabetes diagnosis, we further adjusted for beverage intake before the diagnosis in the multivariable model.

In the current study, we tested the proportional hazards assumption by using a likelihood ratio test comparing models with and without multiplicative interaction terms between beverage consumptions and calendar year, and we did not find evidence of violation of the assumption. Tests for trend were performed by assigning a median value to each beverage consumption category as a continuous variable. To examine the dose-response relationships between beverage intake and the outcomes, we used restricted cubic spline regression with three knots. Tests for non-linearity were based on the likelihood ratio test comparing two models: one with only the linear term and the other with the linear and the cubic spline terms.

We estimated the association of substituting a serving of one beverage for another by including both as continuous variables in the same multivariable model. Differences in their β coefficients were used to calculate the hazard ratios for the substitution effects, and their variances and covariance matrix were used to derive the 95% confidence intervals for the point estimate.

Several sensitivity analyses were conducted to test the robustness of our findings. First, we restricted our analyses to adults with incident type 2 diabetes by excluding those with prevalent diabetes at baseline. Second, we excluded deaths that occurred within four years after the diabetes diagnosis to examine whether the results were impacted by reverse causation bias. Third, a four year and eight year lag were placed between the assessment of beverage intake and outcome incidence, respectively. In these analyses, beverage intake was used to predict disease occurring four years or eight years later. Fourth, given that weight change can be an intermediate outcome, in our final model we adjusted for BMI before the diabetes diagnosis, instead of change in BMI before to after the diagnosis to examine the robustness of our observed associations. Fifth, we examined potential confounding from measures of socioeconomic status by adding partner’s education and self-rated socioeconomic status to the final model. Sixth, we used beverage intake assessed before the diabetes diagnosis instead of the cumulative average after diagnosis to evaluate whether changes in consumption pattern immediately after the diagnosis might impact the associations of interest. Seventh, as it is likely that participants might quit drinking unhealthy drinks immediately after the diabetes diagnosis, we skipped the first food frequency questionnaire after the diagnosis and used the rest to calculate cumulative averages and re-examine the associations. Eighth, we also conducted a sensitivity analysis excluding current and former smokers to further reduce confounding by smoking status. Ninth, we performed an analysis restricted to adults with asymptomatic type 2 diabetes to assess the impact of diabetes screening on associations of interest. Tenth, we controlled for the number of diabetes related symptoms as a measure of disease severity. Lastly, to reduce potential confounding by glucose control, we further adjusted for the self-reported levels of glycated hemoglobin HbA_1c in a subset of the participants (n=5192).

All statistical analyses were performed with SAS software, version 9.4 (SAS Institute, Cary, NC). Two sided P<0.05 was considered statistically significant.