Statistical analyses

Longitudinal analyses were conducted to identify associations between beverage intakes, NARs, MAR, or energy intake and height. A linear mixed model was fit for each dietary variable under study, with the dietary variable as a predictor and height as the outcome. Because height is approximately normally distributed, traditional Gaussian linear mixed-effects models were used, with random slopes and intercepts used to account for the correlation between repeated height measurements on a single subject. Linear mixed models allow for participants to be included in the model if the subject has full data available (dietary intake information and height outcome information) for ≥1 time point. The models were fit by applying maximum likelihood estimation. All of the modeling results were obtained by using PROC MIXED in SAS (version 9.4; SAS Institute, Inc.).

We were interested in the associations between dietary intakes and height throughout childhood and adolescence. To improve statistical power, for each particular dietary variable, we included all participants and all available time points in a single model examining the effect of that dietary intake on height. In order to accommodate changing growth patterns starting in early adolescence and the differences in growth patterns for female and male participants, a change point was incorporated into the model at age 15.5 y for males and 13.5 y for females (17). The change point allowed us to model 4 different slopes for the association between height and time: males before 15.5 y, males after 15.5 y, females before 13.5 y, and females after 13.5 y.

A baseline model was developed that allowed for 4 slopes for the association between height and time for the 4 different groups described above by using indicators to represent slopes for 3 groups and defining a reference intercept and slope for the fourth group. This modeling framework allowed for different slope and intercept variables for the 4 relevant partitions of the data. On the basis of a visual inspection of scatterplots, age 13.5 y was a clear change point for female subjects, whereas the deceleration in growth appeared to be more gradual for males. Therefore, 4 different change points were considered for male participants: ages 13.5, 15, 15.5, and 16 y. After comparison of the log-likelihoods for the baseline models with the 4 change-point options for male participants (and the change point for female participants held constant at age 13.5 y), the change point at age 15.5 y for males was shown to have substantially better fit than the other 3 change-point options for the male participants.

In the initial set of models, a covariate representing intake of a beverage, NAR, MAR, or energy intake was added to the baseline model to identify if the beverage, NAR, MAR, or energy intake was associated with a difference in height over time. In a second set of models, each covariate representing a dietary intake and the interaction between that intake and the indicator for sex was used to determine if such an association differed between females and males. In order to decide between the 2 sets of models, we examined the P value for the interaction between the intake covariate and the indicator for sex. In addition, for each intake variable under study, we compared the Bayesian Information Criterion (BIC) of the model with and without the interaction to determine which exhibited better penalized fit. Models were then adjusted for potential confounders, including MAR, energy intake, and baseline socioeconomic status (SES; i.e., low, middle, and high based on household income and mother's educational attainment). The effect of beverage, the selected NAR, MAR, or energy intake on height for a particular model was summarized, whereas the fixed effects estimates from the baseline covariates were omitted for clarity.