Chapter6.tex

\chapter{Data application}
\label{chap:data}
\section{Introduction}
Obesity has become one of the most burdensome public health issues that faces the United States today \cite{surgeon2001}.  This epidemic impacts children and adults alike.  A recent study estimates that about 32\% of children in the U.S. are overweight or obese \cite{ogden2008}. This alarming statistic has been a call to action for public officials and parents to change how we live and eat. However, encouraging healthy diets and physical activity through increased education and community involvement has seen limited success \cite{hardeman2000,lombard2009}.  It is generally accepted that genetic, metabolic, and environmental factors in addition to behavioral factors play a role in regulating a child's weight \cite{surgeon2001}. Yet, the exact determinants of childhood obesity are poorly understood.

As it is difficult to change the behavior of individuals, recent research has focused on other potential causes and venues for intervention. One hypothesis is that exposure to hormone-altering chemicals prenatally and during early postnatal development disrupt and reprogram the metabolic system to favor weight gain \cite{tuma2007}. A number of rodent studies suggest that prenatal exposure to endocrine disrupters such as bisphenol A (BPA) is associated with increased body weight after adjusting for exercise and amount and type of food \cite{rubin2001,akingbemi2004,hiyama2011, howdeshell1999,miyawaki2007,somm2009,wei2011,xu2011}. A small number of cross-sectional human studies suggest a positive association between BPA and obesity in children and adults \cite{carwile2011,shankar2012,trasande2012}. However, only one study has investigated human prenatal BPA exposure and they found a significant negative association between prenatal exposure and body mass index at age 9 but only in girls \cite{harley2013}.
 
 To fully understand the potential impact of prenatal chemical exposure on childhood growth development, it is necessary to observe the height and weight of children during the time between infancy and adolescence. Quality longitudinal growth studies are becoming more common, but most data analyses do not fully take advantage of the longitudinal nature of the data and only look at two cross-sectional points to determine change over time. In the past few years, researchers have noted the importance of the entire growth trajectory of individuals with a view to determine if there are distinct growth trajectory patterns \cite{pryor2011,carter2012,li2007,garden2012}. In this case, we determine whether growth trajectory patterns are associated with early factors such as prenatal BPA exposure.

\subsection{Body mass index trajectories}
Quantifying and evaluating children's growth involves measuring both height and weight. Body mass index (BMI), calculated as weight (kg)/height$^{2}$ (m$^{2}$), is used as an indicator of total body fat for most individuals \cite{roche1981}. BMI is also used to clinically classify individuals as overweight or obese. The Center for Disease Control (CDC) publishes clinical charts that include BMI-for-age reference growth curves in terms of percentiles and z-scores using cross-sectional data from National Health and Nutrition Examination Survey (NHANES). The z-scores are not intended for longitudinal analysis as the CDC performs different normalizing transformations for each age in months to resolve the skewed sex and age-specific distributions for BMI. Since the reference data is cross-sectional, the charts may not reflect typical age-related patterns of BMI change. Therefore, we concentrate on the original BMI trajectories rather than of BMI z-scores over time.

One way to detect alterations in growth patterns is to apply cluster analysis methods to childhood BMI trajectory data. Then, additional analysis can be completed to estimate the relationship between growth patterns and early life factors such as exposure to BPA as well as other important factors such as maternal BMI, maternal smoking, maternal gestational weight gain, maternal age, birth weight, and duration of breastfeeding which have been suggested to be potential predictors of increased risk of a high-rising BMI growth as well as the early and late onset of obesity \cite{pryor2011,carter2012,li2007}.

\subsection{Clustering}
The clustering methods used in many longitudinal applications including growth trajectories are typically based on a finite mixture model, which is a probabilistic model for representing groups within the overall population that occur with different frequencies but are not known \emph{a priori}. In the simplest form, the mixture density for the random outcome vector $\B y$ is a weighted sum of $K$ group densities written as
$$f(\B y|\BS\theta) = \sum^{K}_{k=1}\pi_{k}f_{k}(\B y|\BS\theta_{k})$$
where $\pi_{k}$ is the probability of belonging to the $k$th group and $\pi_{k}>0$ for all $k=1,..,K$ and $\sum^{K}_{k=1}\pi_{k}=1$. The probability densities $f_{k}$ are typically assumed to be multivariate Gaussian for continuous outcomes such as BMI. With longitudinal data, the main goal is to study the change over time so the mean within each group is modeled conditional on the observations times $\B t.$  This model is the basis for many popular methods and software that have rapidly become commonplace in the growth trajectory literature.

The Proc Traj add-on package for SAS \cite{jones2001}, used in the BMI literature \cite{pryor2011,carter2012}, fits group-based trajectory models. These models are finite mixture models such that for a continuous outcome measure, the group distributions are assumed to be multivariate Gaussian with a polynomial mean function and the repeated measures within an individual are assumed independent conditional on the group membership. Thus, the model assumes that the outcome data $\B y_{i} = (y_{i1},...,y_{im_{i}})$ observed at times $\B t_{i} = (t_{i1},t_{i2},...,t_{im_{i}})$ for subject $i$ who belongs to group $k$ is generated according to
$$\B y_{i} = \beta_{0k}+\beta_{1k} \B t_{i} + \beta_{2k}\B t_{i}^{2}+ \beta_{3k} \B t_{i}^{3} + \BS\epsilon_{i}$$
where $\BS\epsilon_{i}\sim N(0,\sigma^{2}I)$ for a cubic mean model. Additionally, the group probabilities can be modeled using a generalized logit function to allow time-stable covariates determine group membership (setting $\BS\gamma_{K}=0$ for identifiability),
$$\pi_{k}=\frac{\exp(\BS\gamma_{k}^{T}\B w_{i})}{\sum^{K}_{l=1}\exp(\BS\gamma_{l}^{T}\B w_{i})}$$
where $\B w_{i}$ is a design vector based on time-stable covariates. 

Another software program, Mplus \cite{muthen2010}, commonly used in the literature \cite{li2007,garden2012}, fits a generalization of this mixture model, termed a growth mixture model, which allows random effects in the mean structure to account for within-group variation that is ignored in the Proc Traj model specification. Assuming a distribution for the cluster-specific slope and intercept coefficients attempts to model within-individual dependence that is inherent in repeated measures. These random effects can be incorporated into the distribution of the errors. Thus, the observed data for subject $i$ is assumed to be generated by
$$\B y_{i} = \beta_{0k}+\beta_{1k} \B t_{i} + \beta_{2k}\B t_{i}^{2}+ \beta_{3k} \B t_{i}^{3} +\BS \epsilon_{i}$$
where $\BS\epsilon_{i}\sim N(0,\BS\Sigma_{ik})$ with $\BS\Sigma_{ik}=\BS\Lambda_{i}\BS\Psi_{k}\BS\Lambda_{i}^{T}+\BS\Theta_{k}$, where $\BS\Theta_{k}$ is the covariance matrix of the random errors, $\BS\Psi_{K}$ is the covariance matrix for the random effects, and $$\BS\Lambda_{i}=\left(\begin{array}{cccc}1&t_{i1}&t_{i1}^{2}&t_{i1}^{3}\\ 1&\vdots&\vdots&\vdots\\ 1&t_{im_{i}}&t_{im_{i}}^{2}&t_{im_{i}}^{3} \end{array} \right)$$
This model is equivalent to the Proc Traj model when $\BS\Psi_{k}=\B 0$ and $\BS\Theta_{k}=\sigma^{2}\B I$ for all $k=1,...,K$. It is important to note that the random effects indirectly specify a complex, non-stationary covariance structure. This is problematic since misspecifying the covariance structure in finite multivariate Gaussian mixture models can cause bias in the mean estimates and result in an incorrect number of chosen groups as seen in Chapter \ref{chap:misspecify}. Special attention should be paid to modeling the dependence inherent in repeated measures of an outcome over time. I recommend fitting a model using a well-known stationary correlation structure such as exchangeable or exponential correlation before fitting a complex non-stationary model via random effects. 

These model-based methods are typically understood to group individuals such that members of the same group share a similar pattern of change over time \cite{garden2012}. This suggests that two individuals with similar patterns of change but different vertical levels would be grouped together. However, model parameters and subsequently group membership are estimated by maximizing a likelihood function, which is based on Euclidean distance between observed data and group means if normality is assumed (see Chapter \ref{chap:motivate}). This process results in groups being primarily determined by the level despite subtle differences in the shape especially if the level and shape of a trajectory are independent or weakly dependent. Therefore, these methods do not directly group trajectories on shape and resulting groups may include growth trajectories of the same level but different shapes over time. This translates into estimated group memberships and means not accurately representing shape groups present in the data.

Due to this misunderstanding, researchers present the results of a mixture model analysis by discussing the shape of the estimated mean growth pattern for each group \cite{pryor2011,carter2012}. However, the mean only represents the average shape of all trajectories in the group; perhaps no one individual follows the path of the mean especially when the groups are not homogeneous in terms of shape. Additionally, the estimated relationships between risk factors and the resulting groups are deceitful when the group descriptions are inaccurate and misunderstood. Care needs to be taken when making conclusions based on these methods especially in describing the shape of the trajectories within each group.  

The goal of this research is to group individuals on the basis of the shape of their growth trajectory, so it is appropriate to compare the new methods proposed in Chapter \ref{chap:methods} with those used in practice on a real data set. The vertically shifted mixture model is based on the same foundation as those described above. The main difference is that rather than fitting a model to the original data, the input for the model is the transformed or vertically shifted data vector. Each individual's mean BMI is subtracted from all of their measured BMI values. In effect, each individual is normalized to have mean zero and the level is removed without eliminating the variability at any time point. By removing the mean prior to fitting a multivariate Gaussian mixture model, we allow the clustering method to focus directly on the shape rather than the level. Thus, the resulting groups can be honestly summarized by describing the shape of the mean trajectory of the group and estimated associations with risk factors can be interpreted accurately in terms of shape groups. I illustrate the difference between the methods with longitudinal BMI data on a sample of children living in the Salinas Valley, CA.


\section{Data}
The Center for the Health Assessment of Mothers and Children of Salinas (CHAMACOS) Study is a longitudinal birth cohort study designed to assess the health effects of pesticides and other environmental exposures on the growth and development of low-income children living in the agricultural Salinas Valley, CA \cite{eskenazi2004,eskenazi2005}. Of 601 pregnant women enrolled in the study in 1999-2000, a total of 527 mostly Latino mother-child singleton pairs were followed through a live-birth delivery and 327 pairs continued to be followed through the 9-year interview. Baseline maternal characteristics were measured at the start of the study and maternal urine and blood samples were taken twice during pregnancy and then again shortly after delivery to measure levels of pesticide and chemical exposure. Child height and weight were measured without jackets and shoes by trained staff at interviews that occurred at birth and when the child was approximately 1, 2, 3 $\frac{1}{2}$, 5, 7, and 9 years of age. BMI was calculated as weight (kg) divided by height squared ($m^{2}$) for ages 2 and over. The exact age of the child was also recorded due to variability in the interview times. For this thesis, we limit our analysis to 303 children who have at least four recorded BMI measures. Details of the study are published elsewhere \cite{eskenazi2003}. All study activities were approved by the Committee for the Protection of Human Subjects at the University of California, Berkeley. 

For illustrative purposes, we limit the discussion to two baseline risk factors, maternal pre-pregnancy BMI and BPA exposure. One of the strongest known predictors of a child's BMI is maternal pre-pregnancy BMI. At the start of the study, the maternal height was measured and used together with self-reported pre-pregnancy weight to calculate the BMI. During the first and second half of pregnancy, maternal BPA exposure was measured via urine samples. See \textcite{harley2013} for details of how the concentrations were measured.

Characterizing BMI growth trajectories of children and distinguishing distinct patterns are of great importance to public health officials. We use this data set to illustrate the differences between the clustering methods and how they work in practice with real data.

\section{Cluster methods}
I use the standard and novel methods described in this thesis and compare the clustering groups and the inference on baseline factors that may be related to cluster membership. Specifically, we fit a multivariate Gaussian mixture with conditional independence like Proc Traj (Model 1), a multivariate Gaussian mixture with exponential correlation (Model 2), and the vertically shifted mixture model proposed in Chapter \ref{chap:methods} (Model 3). There is some variability in the observation times, but the growth is rather gradual so the independence assumption may be the best approximation for the covariance of the transformed data in Model 3.

All of these model-based methods require a mean structure. I use a B-spline basis to model the mean function. Based on the visual inspection of the data, I determine that quadratic basis functions with one internal knot should be sufficient to model the complexity given the limited number of time points.  The internal knot is place at the median age of 5 years. 

To compare the methods in terms of baseline factor relationships, I fit the model with two different factors: maternal pre-pregnancy BMI, which is known to impact the child's BMI, and the log base two transformed maternal BPA exposure during pregnancy, which is hypothesized to impact growth trajectories. The relationships are estimated separately from each other but simultaneously with the other model parameters.

To estimate parameters for all of the models, I use maximum likelihood estimation via the EM algorithm. The algorithm is run five times for each model with random starts and the final estimate are from the fit with the highest likelihood. For each model, the number of groups is chosen by fitting the models for $K=2,3,4,5,6$ and choosing the value of $K$ that minimizes the BIC. To make inferences about the parameters, I estimate the robust sandwich standard error for the parameter estimators \cite{white1982}. Odds ratios for the baseline factors are calculated by exponentiating the coefficients and are presented with exponentiated confidence intervals.

The cluster groups are visualized by plotting individual BMI trajectories colored according to the group assignment made by maximizing the posterior probability. The group means are represented in the adjoining panel. For the vertically shifted model, the group means are shifted such that the mean BMI at the age 2 equals the average BMI of the individuals in the particular group at the age 2 interview. 

\section{Results}
Figure \ref{fig:6-1} shows the clustering results for the three mixture models with maternal pre-pregnancy BMI as a baseline factor that impacts group membership. The chosen number of groups for each model is $K = 5$, $K=4$, and $K=5$ for Model 1, 2, and 3, respectively. Model 1 as compared to Model 2 requires more groups to model the variability within the BMI trajectories due to the limited correlation structure. Neither Model 1 nor Model 2 group individuals based on the shape of their trajectories, but it is hard to see this without comparing the plot of groups trajectories to those from Model 3. The plot from the last model brings attention to individuals with high BMI levels at age two with relatively stable trajectories over time and individuals with moderate BMI levels whose BMI drastically increased over the follow-up period. Both of these types of individuals are categorized with trajectories with very different shapes in the other two models.

It is interesting to note that the mean shapes do not drastically change between the three models even though group membership changes (right side of Figure \ref{fig:6-1}). Thus, the number of membership changes is relatively small in this case so as to not to change the shape of the means for this data set. Even though it is tempting to compare the starting levels of the mean curves, it is important to realize that the level was removed to focus on the shape in Model 3 and the artificial levels for presentation are based on the average BMI at age 2 within the groups. Thus in Model 3, group 1 includes more individuals with lower BMI at age 2 in comparison to the other models. 

\begin{figure}[ht]
\centering
\includegraphics[width=6.5in]{Chp6BMI}
\caption{Clustered BMI trajectories colored according to the group assignment made by maximizing the posterior probability and group mean functions for three mixture models (Model 1: independence, Model 2: exponential, Model 3: vertically shifted) with maternal pre-pregnancy BMI as the baseline factor.}
\label{fig:6-1}
\end{figure}

Table \ref{tab:6-1} presents the estimated odds ratios comparing group membership odds for a one unit increase in maternal pre-pregnancy BMI. The reference group has the highest numerical label with the most stable shape over time. There seems to be a consistent monotone trend in the odds ratios for pre-pregnancy BMI as the derivative of the trajectory shape increases, but the inference depends on the model used. We see that it highly predicts group memberships in Model 1 as the individuals are divided mainly by BMI level. However, generalizing the correlation structure in Model 2 increases the standard error enough to decrease the strength of the evidence. While the point estimates are similar for Model 3, the standard errors are slightly smaller but not enough to provide strong evidence that maternal pre-pregnancy BMI also impacts the growth pattern over time. 

% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Thu Mar 14 16:06:36 2013
\begin{table}[ht]
\centering
\begin{tabular}{lccc}
  \thickhline
 & Odds Ratio & 95\% CI& $P>|z|$ \\ 
  \hline
Model 1: Group  1 & 1.32 & (1.2, 1.45) & 0 \\ 
  Model 1: Group  2 & 1.24 & (1.14, 1.35) & 0 \\ 
  Model 1: Group  3 & 1.15 & (1.04, 1.27) & 0.006 \\ 
  Model 1: Group  4 & 1.11 & (1.02, 1.21) & 0.017 \\ 
  Model 1: Group  5  (ref.) & - & - & - \\ 
  Model 2: Group  1 & 1.3 & (0.85, 1.98) & 0.232 \\ 
  Model 2: Group  2 & 1.23 & (0.88, 1.73) & 0.23 \\ 
  Model 2: Group  3 & 1.11 & (0.88, 1.39) & 0.387 \\ 
  Model 2: Group  4  (ref.) & - & - & - \\ 
  Model 3: Group  1 & 1.28 & (0.92, 1.79) & 0.136 \\ 
  Model 3: Group  2 & 1.24 & (0.93, 1.67) & 0.143 \\ 
  Model 3: Group  3 & 1.18 & (0.89, 1.55) & 0.247 \\ 
  Model 3: Group  4 & 1.04 & (0.85, 1.28) & 0.7 \\ 
  Model 3: Group  5  (ref.) & - & - & - \\ 
   \thickhline
\end{tabular}
\caption{Odds ratio estimates and confidence intervals for a one unit increase in maternal pre-pregnancy BMI comparing each group to the reference group for the three mixture models (Model 1: independence, Model 2: exponential, Model 3: vertically shifted). }
\label{tab:6-1}
\end{table}

Figure \ref{fig:6-2} shows the clustering results for the three mixture models with maternal BPA exposure during pregnancy as a baseline factor that impacts group membership. The chosen number of groups for each model is $K = 5$, $K=4$, and $K=5$ for Model 1, 2, and 3, respectively. The clustering looks very similar to those from the models with pre-pregnancy BMI. The main differences occur with Model 3. The group means slightly differ in shape from those with maternal BMI. This is natural since all of the parameters include group membership are simultaneously estimated. Therefore, baseline factors have a minor influence on the group membership for a trajectory in the overlap between two groups. This change of group can be enough to change the mean curve.

\begin{figure}[ht]
\centering
\includegraphics[width=6.5in]{Chp6BPA}
\caption{Clustered BMI trajectories colored according to the group assignment made by maximizing the posterior probability and group mean functions for three mixture models (Model 1: independence, Model 2: exponential, Model 3: vertically shifted) with maternal BPA exposure during pregnancy as the baseline factor.}
\label{fig:6-2}
\end{figure}

Table \ref{tab:6-2} gives the estimated odds ratios comparing group membership odds for a one unit increase in log base 2 transformed maternal BPA exposure. The reference group is again the last group with the most stable trajectory shape. In generalizing the correlation structure (Model 2 vs. Model 1), the estimated relationship between BPA and group membership changes. Specifically, the odds ratio estimates for group 2 switches from 0.9 in Model 1 to 1.24 in Model 2. This suggests that increased BPA increases the probability of being in group 2 in comparison to group 4. This magnitude change is maintained by the vertically shifted model (Model 3), and the standard errors are smaller relative to Model 2 to provide more evidence that there is small effect of BPA exposure on the shape of the growth curve over time. 

% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Thu Mar 14 16:08:30 2013
\begin{table}[ht]
\centering 
\begin{tabular}{lccc}
  \thickhline
 & Odds Ratio & 95\% CI & $P>|z|$ \\ 
  \hline
Model 1: Group  1 & 0.84 & (0.58, 1.21) & 0.345 \\ 
  Model 1: Group  2 & 0.9 & (0.65, 1.25) & 0.542 \\ 
  Model 1: Group  3 & 1.04 & (0.75, 1.46) & 0.804 \\ 
  Model 1: Group  4 & 0.92 & (0.65, 1.29) & 0.618 \\ 
  Model 1: Group  5  (ref.) & - & - & - \\ 
  Model 2: Group  1 & 0.79 & (0.56, 1.11) & 0.181 \\ 
  Model 2: Group  2 & 1.24 & (0.84, 1.83) & 0.285 \\ 
  Model 2: Group  3 & 1.01 & (0.73, 1.39) & 0.956 \\ 
  Model 2: Group  4  (ref.) & - & - & - \\ 
  Model 3: Group  1 & 0.8 & (0.49, 1.29) & 0.356 \\ 
  Model 3: Group  2 & 1.33 & (0.95, 1.87) & 0.095 \\ 
  Model 3: Group  3 & 1.01 & (0.71, 1.43) & 0.948 \\ 
  Model 3: Group  4 & 0.88 & (0.6, 1.29) & 0.52 \\ 
  Model 3: Group  5  (ref.) & - & - & - \\ 
   \thickhline
\end{tabular}
\caption{Odds ratio estimates and confidence intervals for a one unit increase in maternal BPA exposure during pregnancy in log base 2 units comparing each group to the reference group for the three mixture models (Model 1: independence, Model 2: exponential, Model 3: vertically shifted). }
\label{tab:6-2}
\end{table}

\section{Discussion}
To explore the population in terms of their heterogeneous growth patterns over time, we used multivariate finite mixture models. Rather than averaging out all of the interesting growth patterns, mixture models allow for a finite number of relationships. Standard mixture models group these relationships based on the feature that dominate the variability. If this feature is the level of the trajectory, the estimation of the mixture model may not highlight interesting patterns over time. Therefore, we use a simple adjustment to remove the vertical level and only focus on how the children's growth changes over time.

For this data set, the mean functions do not drastically change between the standard model and the vertically shifted model, but the group memberships change. Moving individuals between groups impacts the interpretation of baseline factors and their association with the probability of being in a particular group. The well-studied association of pre-pregnancy BMI and child's BMI level \cite{whitaker1997} is supported with this data set. However, this data set does not provide strong evidence for an association with the shape of the growth pattern. This is an important distinction to make as genetic factors may impact early childhood BMI but not necessarily determine the growth trajectory over time. While this data set does not provide strong evidence to suggest that prenatal exposure to BPA is associated with the level or shape of the growth pattern, the shift in the estimates between Model 1 and Model 3 for the second group of late gainers highlights the importance of investigating the factors that impact level and those that impact shape separately.