Identifying developmental trajectories of body mass index in childhood using latent class growth (mixture) modelling: associations with dietary, sedentary and physical activity behaviors: a longitudinal study

Study population and design

The ChecKid study is a dynamic cohort study of primary school children in the ages of 4 to 12 years in the city of Zwolle, the Netherlands. The objectives of ChecKid are to investigate trends in overweight and to examine health behaviors related to childhood overweight and obesity and determinants of these behaviors within families, schools and neighborhoods. ChecKid is part of an integrated approach in which quantitative and qualitative monitoring research and environmental scans support the development, implementation and evaluation of tailored community wide interventions.

This study commenced in 2006 and in 2006, 2009 and 2012 every primary school in the city of Zwolle (n = 51/43/43) was invited to participate in the study. When a school agreed to be included in the study, all children attending the school (4-12 years) and their parents were invited to participate by means of letters distributed via the schools. Respectively 80 %, 79 % and 81 % of the schools were willing to participate in 2006, 2009 and 2012. When schools did not want to participate, it was mostly because of other priorities. Participating schools were equally spread over all neighborhoods in Zwolle. Measurements took place every three years, during 3 weeks, in October and November of 2006, 2009 and 2012, thus indicating a follow up time of six years.

Eligible children included those who had an anthropometric measurement (height and weight) and whose parents filled in a self-report questionnaire about the reported behaviors in 2006, in 2009 and in 2012. Total response in 2006 was 4,072 children (49 %) of whom anthropometric measurements and questionnaires were available, in 2009 this was 3,026 (35 %) and in 2012 5,849 (61 %). The dynamic design of the study entailed that every three years the total population at every school was invited to participate in the study. Due to this design only children who were in class 1,2 and 3 in 2006 (n = 1,599; age range 4.0 to 7.0 years) can be included as a cohort. No new participants were recruited in 2009 and 2012. Exclusion criteria for participants were not being sufficient in the Dutch language, being older than 13 years of age at follow up, attending special education, and not living in the city of Zwolle.

For the estimation of the trajectories at least two measurements were included. Because we felt that it was necessary to have a “true” measurement at baseline instead of an estimation of the trajectory, we used data of all children with a baseline measurement in 2006. Besides baseline, a minimum of 1 measurement in 2009 and/or in 2012 was included. Because of high response rates in 2006 and 2012, none of the participants had a missing measurement during these years, meaning we only have missing measurements in 2009 in our dataset. Two hundred forty-three (39.6 %) of the participants had only 2 measurements. This leaves us with a sample of 613 children, who were included and completed follow-up. Mean age at baseline is 5.1 years (range 4.0 to 7.0 years), at the second wave it is 8.1 years (range 7.0 to 9.8 years) and in 2012 11.1 years (range 10.0 to 12.9 years).

All parents gave written informed consent and medical ethical approval was obtained from the Medical Ethics Committee of the VU University Medical Centre.

BMI

Height, weight and waist circumference were measured according to protocol [40] by trained students. BMI was calculated as weight in kilograms divided by height in meters squared. The children’s BMI cut-off points defined by Cole et al. were used to define thinness, healthy weight, overweight and obesity [41, 42]. Because an increase in BMI during childhood is a feature of natural growth, we calculated individual age and gender specific BMI standard deviation scores (SDS) using the data from the National Growth Study 1997 as the reference population [7].

Demographic variables and health behaviors

We obtained information on the child’s health behavior (dietary, physical activity and sedentary behavior) during weekdays and weekends separately by self-administered questionnaires completed by one of the parents. We used information on the child’s health behavior during weekdays, as we were especially interested in finding indications for interventions that could possibly be implemented or supported in a school setting. Participation in organized sports is the only variable which was assessed during weekdays and weekends combined (>1 h a week) because a lot of sports activities are pursued during the weekend. The questionnaire for parents included socio-demographic variables such as the child’s age, gender, postal code, ethnicity (assessed by country of birth of both parents) and socio-economic status (SES) (assessed by educational level parents). Parents’ self-reported body weight and height data were used to calculate their BMI. Existing validated questionnaires on health behavior were used for the design of this questionnaire [43, 44]. As the questionnaires were originally designed for secondary school students, the questionnaires were adjusted for primary school children in this study [44]. Details of the ChecKid questionnaires are described elsewhere [45].

Dietary intake was measured by exploring frequency and amount of fruit intake, vegetable intake and frequency and amount of sugared drink consumption. Consumption of vegetables was dichotomized as eating vegetables daily (5 days a schoolweek) vs not daily (<5 days a schoolweek). Fruit consumption was dichotomized as less than 2 portions a weekday (<2 portions) vs 2 or more portions a weekday (≥2 portions). This was done according to guidelines issued by the Dutch Nutrition Centre [46, 47]. The consumption of sugared drinks (fruit juices, soft drinks and sweetened milk drinks) was dichotomized as 3 or less glasses of sugared drinks per weekday (≤3 glasses) versus more than 3 glasses of sugared drinks per weekday (>3 glasses), according to guidelines used in previous studies [45, 48].

Sedentary behavior was measured by investigating total screen time (TV viewing and using the computer) and TV viewing separately. Parents were asked to report frequency and duration of time (in categories) spent watching TV or using the computer. Averages were calculated by multiplying the number of days that the child spent on the behavior by the mid-category values of duration of the item in 5 categories: <0.5, 0.5-1, 1-2, 2-3, and >3 h a day, and dividing this by the number days per week. Screentime and TV viewing were both dichotomized as less than 2 h (<2 h) or 2 h or more (≥2 h) per weekday.

Physical activity was measured by investigating time spent on outside play and in organized sports. Parents were asked to report frequency and duration of time (in categories) spent on outside play or participation in organized sports. Averages were calculated by multiplying the number of days that the child spent on the behavior by the mid-category values of duration of the item in 5 categories: <0.5, 0.5-1, 1-2, 2-3, and >3 h a day, and dividing this by the number of days per week. Outside play was dichotomized according to the Dutch standard for exercise for children as playing outside 1 h or more (≥1 h) and less than 1 h (<1 h) per weekday [49]. Participation in organized sports was dichotomized as less than 1 h (<1 h) and 1 h or more (≥1 h) per week.

Statistical analyses

Statistical analyses were conducted using the Mplus 7.11 and the PASW 20.0 software packages. Descriptive statistics were used to investigate population characteristics and differences in behaviors over time, using Chi square and t-tests.

Estimation of BMI trajectories

Longitudinal BMI trajectories were analyzed with latent class growth mixture modeling [21, 25, 50, 51], LCGMM. LCGMM is a contemporary longitudinal technique in obesity research. The underlying aim of the technique is to capture the heterogeneity in BMI development over time in k number of subgroups (or classes), each with a distinct BMI trajectory shape. Each class has its own growth parameters (i.e., intercept and linear slope), where the intercept represents the average BMI-value at baseline and the slope is indicative of the rate at which the BMI-values change over time.

Various LCGMM models were run before we chose a final model. First, several linear LCGMM with fixed intercept and slope variance within classes were investigated. Next, quadratic slopes were added to the model allowing for curved developmental patterns. The timescores were based on the follow-up measurements and were held equal across time (0,1,2). To determine the optimal number of classes, a common forward procedure [26, 35] was performed, starting with a one class solution (i.e., assuming all individuals follow a similar BMI trajectory over time), then adding more classes one at a time to investigate whether or not the model fit improves due to the additional class(es). To explore the possibility of different trajectories by gender, these steps were initially performed separately in boys and girls. Because the trajectories were very similar in nature in both girls and boys, as also found in Pryor et al. [31] and Magee [34], the analyses were performed on the total sample. The results of these analyses will be reported in this paper.

To compare fit between subsequent models, we used several model fit indices and other criteria, as suggested by the literature [52]. First, we used two model fit indices; the Bayesian Information Criterion (BIC) and the Bootstrapped Likelihood Ratio Test (BLRT). The BIC [53, 54] is common in latent class models and considers both the likelihood of the model as well as the number of parameters in the model. A lower value is indicative of a better fitting model. The BLRT [55] provides a p-value, comparing a model with k-1 classes with a model with k classes. The BLRT has been shown to be a very consistent indicator of the optimal number of classes [56]. To further look at the optimal number of classes, we also took the posterior probabilities into account. Children were assigned to their most-likely class based on these probabilities. The probability for the class to which a child is assigned to should be considerably higher that the probabilities for the other classes to make the classes easily distinguishable from each other.

Finally, we considered the usefulness of the models in practice by examining the shape of the trajectories and the sample size of each class. We ran several quadratic models, forcing all within-class variances at zero (Nagin type). Although the BIC for the 2-class model was slightly lower (3579.191 versus 3607.081 for the linear model), the class size and the trajectory shapes were similar. The linear GMM included estimated (free) intercept variances within classes (Mplus default) and the slope variance was forced at zero (Nagin type model). We based these choices on the comparison of model fit indices described earlier. We rejected models with clinically uninterpretable classes as well as models including classes with a sample size <1 % of the total sample [26].

Associations between health related behaviors and BMI trajectories

The first step of the analyses provides us with a classification of the children into distinct subgroups, based on their BMI trajectory. This classification is coded as a categorical variable with k number of categories, or classes. To study the possible associations between the classes and behaviors at baseline and the possible associations between the classes and behaviors at last measurement, several logistic regression models were conducted, in which the dichotomized behavior variable is the outcome measurement and class membership the independent variable. First crude analyses were performed. Second, adjusted analyses were carried out, controlling for potential confounding effects of gender, SES and ethnicity. Separate logistic regression models were performed for associations with behaviors at baseline in 2006 and for associations with behaviors in 2012.

Although our primary aim was to examine associations between behaviors at baseline and the trajectories to explore the predictive value of the behaviors, we also examined associations between behaviors at last measurement in 2012 and the trajectories. We wanted to explore the cross-sectional association between current behavior and BMI trajectory as parallel outcome. As behavior is known to track throughout life [18–20] it is modifiable, and as children become more autonomous in their behavior with age, we therefore think it is valuable to explore associations between health behavior at last measurement and the trajectories.

Missing data

Missing value percentages in the behavior variables were 13.9 % (TV viewing), 14.2 % (vegetable consumption), 14.5 % (outside play), 15.2 % (participation in organized sports), 15.5 % (screentime), 20.6 % (fruit consumption) and 27.7 % (consumption of sugared drinks) in 2006, and 0.8 % (consumption of vegetables), 1.1 % (consumption of sugared drinks), 3.3 % (fruit consumption), 4.4 % (TV viewing), 5.4 % (outside play), 5.4 % (participation in organized sports), and 12.7 % (screentime) in 2012. Because overweight children and children with a low SES background more often had a missing value on these variables, missing data were not completely at random. The authors used the multiple imputation procedure in PASW 20.0 software to impute the missing data. The variables: BMI of both parents, ethnicity, gender and age were included in the imputation models, and data were imputed ten times using the ‘Predictive Mean Matching ‘(PMM) procedure. For comparison reasons, data were also imputed five, twenty and fifty times using the same method, resulting in a total of 40 imputed datasets. All analyses were performed on the original data and on the imputed datasets. No major differences were found in the direction of the found associations, nor in the slopes or p-values.

Number of BMI trajectories: model selection

Sensitivity analysis

Because the two class and three class model both had equivalent fit, and because it was difficult to determine a clear normative, stable trajectory we decided to explore the three class model in an extra sensitivity analysis. We found a stable high BMI SDS trajectory (n = 212), a decreasing BMI SDS trajectory (n = 346), and a steeply increasing BMI SDS trajectory (n = 55). The three-class model did not yield a clear normative trajectory either, and when we compare the steeply increasing and increasing trajectories in the 3 class model using descriptive analyses of characteristics of the study population and the health behaviors, we find that they are nearly the same, see Additional file 1: Table S1 and Table S2 and Additional file 2: Figure S1. This indicates lower posterior probabilities and unclear class distinctiveness.

After taking into account the substantive interpretation of the two models and objectives of the analyses, a final decision for the model with two classes was made, showing a decreasing BMI SDS trajectory (n = 416) and an increasing BMI SDS trajectory (n = 197). Figure 1 shows mean trajectories of the final model. In our data we could find no children who maintained BMI SDS scores of 0 in 2006 and in 2012, thus supporting our finding lacking a normative group.

Demographics and health related behaviors in the BMI trajectories

The increasing BMI SDS trajectory had a statistically significantly higher percentage of participants of non-western ethnicity and lower SES than the decreasing BMI SDS trajectory. In the increasing BMI SDS trajectory maternal overweight prevalence was significantly higher than in the decreasing trajectory, both at baseline in 2006 and at follow-up. Paternal overweight prevalences did not differ significantly at baseline, but at follow-up the prevalence of overweight fathers was significantly higher in the increasing BMI SDS trajectory.

Associations between BMI trajectories and health related behaviors at baseline

Associations between BMI trajectories and health related behaviors in 2012

After adjustment, the children in the increasing BMI SDS trajectory scored higher odds in 2012 for drinking more than the recommended amount of sugared drinks (OR, 1.2; 95 % CI 0.8 to 1.8), participating in organized sports less than recommended (OR, 1.3; 95 % CI 0.8 to 2.0), and more than recommended use of screentime (OR, 1.3; 95 % CI 0.9 to 1.8) than those in the decreasing BMI SDS trajectory, though these results were also not significant. The OR for maternal (OR, 1.8; 95 % CI 1.2 to 2.6) as well as for paternal overweight status (OR, 1.6; 95 % CI 1.1 to 2.4) was significantly higher for children in the increasing BMI SDS trajectory compared to the children in the decreasing BMI SDS trajectory.

Changes in health related behaviors between 2006 and 2012 in the BMI trajectories

Chi square tests were performed to explore differences in behaviors between 2006 and 2012 (Table 3), and for both trajectories a statistically significant increase was found in participation in organized sports, and in TV viewing/screentime more than 2 h daily, as well as a statistically significant decrease in the consumption of sugared drinks.

In the decreasing BMI SDS trajectory fruit consumption decreased significantly, and outside play increased significantly (see Fig. 2).

BMI trajectories in childhood

There have been a number of studies that have prospectively explored BMI trajectories in the period of (early) childhood [28–34, 56]. We found three studies on BMI trajectories and associations with early life risk factors. One study shows three trajectories of BMI (low-stable, moderate, and high-rising) in children between 5 months and 8 years of age [31]. The low-stable and moderate trajectory groups showed a similar shaped stable developmental pattern, whereas the high rising group was characterized by an increasing average BMI, which exceeded international cut-off values for obesity by age 8 years. Maternal smoking and maternal BMI were important risk factors for the high rising group. These risk factors were also found in another study which also identified three overweight trajectories: early onset overweight, late onset overweight, and never overweight in children from 2 years to 12 years old [33]. The early onset group showed early onset of overweight which persisted throughout childhood, and the late onset group had a moderately high probability of overweight at 2 years, low probability at age 4 and 6 years, but growing probability of overweight after age 8 years. The never overweight group had a low probability of developing overweight throughout childhood. In another study which identified four BMI trajectories in children from 4 to 11 years old (high-risk overweight, early onset overweight, later onset overweight and healthy weight), maternal overweight, maternal smoking and high birth weight were associated with the high-risk overweight trajectory [29]. In a follow-up study the same authors found three trajectories (healthy weight, early onset obesity, and later onset obesity) and found longitudinal associations between the early onset obesity trajectory and sleep duration [34].

In most studies concerning BMI trajectories (in childhood), besides a progressively overweight trajectory and a progressively overweight but stabilizing trajectory, a normative/stable trajectory can be found [30–32, 35]. In our study the two-class model, which had the best fit, did not comprise such a normative trajectory. Instead, most of the sample (68 %) was identified as a decreasing BMI SDS trajectory. We also explored the three-class model in an extra sensitivity analysis, but in this model no normative or stable trajectory was found either.

In a study which examined socioeconomic equalities in childhood trajectories of BMI and overweight a longitudinal tendency was found, showing two longitudinal influences: decreasing BMI in the more advantaged populations (high SES), coupled with simultaneously worsening rates of overweight in the more disadvantaged (low SES) [58]. By age 10-11 socioeconomic differences in mean BMI z-scores already present at age 4-5 had more than doubled, reflecting decreasing mean BMI among advantaged rather than increasing means among disadvantaged children. In our study the decreasing trajectory consists of significantly more children of high SES, supporting this previous finding, thus highlighting the importance of SES in the development of overweight in childhood.

Strengths and limitations of the study

Our study is strengthened by the high participation rates and the equal spread of participating schools in one city, as well as all children being measured during a three week period according to protocol. The prospective repeated measurements add to the strengths of this study. Some limitations of our study can be identified. A potential limitation is that only three measurements of BMI were used to obtain BMI developmental trajectories. The need to impute relative high percentage missing values of low SES and overweight children on the behavior variables at baseline could indicate possible bias, although our analyses on the imputed on non-imputed datasets did not yield different results. The use of parental self-report questionnaires on health behaviors can possibly lead to social desirability and difficulties recalling the actual behaviour of the child, especially if this behavior takes place outside the home [59]. This could lead to over or underestimation of the child’s behavior, possibly influencing the generalization of the results. Recommendations to avoid this type of bias would be the use of more objective measurements such as movement sensors or accelerometers or the use of a diary for measuring dietary intake.

Dichotomization of the behavior variables might be another limitation in this study, reducing the sensitivity of the data. Furthermore, use of BMI SDS scores in a longitudinal study instead of raw BMI scores is possibly a drawback because the within-subject variability in BMI values over time is much smaller for children at the end of the distribution, meaning that standardized scores are less sensitive to changes in overweight/obese and underweight children. At the follow up in 2012, children were involved up to 12.9 years, meaning that pubertal status could possibly have influenced our outcomes. Unfortunately, data on pubertal status was not available in our sample.