Chapter1.tex

\chapter{Introduction}\label{chap:intro}

Longitudinal studies play a prominent role in health, social, and behavioral sciences as well as in the biological sciences, economics, and marketing. By following subjects over time, temporal changes in an outcome of interest can be directly observed and studied. This dissertation studies the task of clustering longitudinal data. Statisticians, machine learning communities, and applied researchers have investigated this problem for other highly structured multivariate data sets such as time series and functional data, and it is important to distinguish between these different types of data. While all of these data structures are ordered by a variable such as time, they have different characteristics that call for different methods. Time series are sequences of data points measured uniformly over time from random processes on only a few units. Examples include monthly unemployment figures, yearly global temperatures, and hourly stock prices over time \cite{shumway2011}. Functional data are assumed to provide information about smooth curves or functions that have been measured on a densely sampled grid. The data include many curves, one for each observational unit, and may be observed with substantial measurement error. Data such as hip angles over walking cycles and height during adolescence that have underlying smooth curves are categorized as functional data \cite{ramsay2002,ramsay2005}. 

In this dissertation, I focus on longitudinal data which includes sporadic repeated measurements of the same outcome observed on many subjects at sparse, potentially irregular times over a long period of time. Long-term ongoing longitudinal studies include the National Longitudinal Survey of Youth, which has attempted to interview individuals annually and then biannually for the past 30 years, and the Framingham Heart Study, which has interviewed three generations of subjects biannually over the past 60 years. These long-term studies have observed each individual as many as 30 times. The majority of longitudinal studies only observe individuals 5 to 15 times during the follow-up period due to limited resources. For the purposes of this dissertation, I often refer to a longitudinal time-ordered vector of outcomes as a trajectory as a reminder that the data vectors are ordered and represent a path of development over time. 

An important question concerns the existence of distinct trajectory patterns. Cluster analysis seeks to separate objects (individuals, subjects, patients, observational units) into homogeneous clusters. There are many ways to cluster multivariate data and there are a multitude of methods and algorithms available to use. Most methods can be categorized into one of two approaches: nonparametric and model-based methods. The first approach makes no assumptions about how the data were generated and produces a sequence of clustering results indexed by the number of clusters $k=2,3,...$ and the choice of dissimilarity measure. The later approach assumes data vectors are generated from a finite mixture of distributions. The clustering results include the conditional probability of group membership and estimates of the parameters for the assumed data generation distribution \cite{mclachlan2000,fraley1998}. Both of these approaches are widely used and the pros and cons of the approaches have been discussed in depth \cite{magidson2002, everitt1981}.

The bulk of the available clustering methods are intended for use on data vectors with exchangeable, independent elements and are not appropriate to be directly applied to repeated measures with inherent dependence \cite{everitt2009}. Recent work to cluster longitudinal data has focused on extending and adjusting the standard clustering methods. Dissimilarity measures have been designed to be sensitive to permutations in the time order by incorporating the relationships between neighboring observations \cite{chouakria2007} or by calculating dissimilarity between coefficients of a functional projection rather than the raw vectors \cite{serban2005, tarpey2003, abraham2003, tarpey2007,hitchcock2007}. In the model-based approach, one can incorporate a functional basis in the mean structure \cite{nagin1999,gaffney1999} and explicitly model the temporal correlation structure \cite{muthen1999,fraley1999,mcnicholas2010} to incorporate the longitudinal nature into the model. 

\section{Longitudinal data}
In classical univariate statistics, each individual or subject gives rise to a single measurement, termed the outcome. In multivariate statistics, that one measurement is replaced by a vector of multiple outcome measurements. In longitudinal studies, a vector of measurements is observed for each subject, but the vector represents one outcome measured repeatedly at a sequence of observation times. Therefore, the data have characteristics of both multivariate and time series data.

Longitudinal data differ from multivariate data in that the time ordering imparts dependence among measurements that is not present in a typical multivariate data set. They differ from classical time series data because the data consists of a large number of independent trajectories that are sparsely and potentially irregularly sampled over time, rather than a few random processes that are uniformly sampled over time. These properties call for methods specific to longitudinal data \cite{diggle2002}. Due to the sparsity, it is necessary to borrow strength between individuals by modeling the mean outcome as a function of explanatory variables. When there are distinct patterns in the data, averaging can be done within the groups.

In this dissertation, I let $Y_{ij}$ represent an outcome random variable and $\B x_{ij}$ be a design vector of length $p$ based on explanatory variables observed at time $t_{ij}$ for individual $i=1,...,n$. I let $\B w_{i}$ be a design vector of length $q$ based on baseline variables observed at or before $t_{i1}$. The set of repeated outcomes for individual $i$ together in a vector is denoted as $\B Y_{i} = (Y_{i1},...,Y_{im_{i}})$ with mean $E(\B Y_{i}) = \BS\mu_{i}$ and $m_{i}\times m_{i}$ covariance matrix $Cov(\B Y_{i}) = \BS\Sigma_{i}$, where the element in the $j$th row and $l$th column of $\BS\Sigma_{i}$ is the covariance between $Y_{ij}$ and $Y_{il}$. The covariance matrix can be written as $\BS\Sigma_{i}=\B V_{i}^{1/2}\B R_{i}\B V_{i}^{1/2}$ where $\B R_{i}$ is a $m_{i}\times m_{i}$ correlation matrix of $\B Y_{i}$ and $\B V_{i}$ is a $m_{i}\times m_{i}$ diagonal matrix with variances as the non-zero elements. In this dissertation, individuals are assumed independent of each other; therefore, $\B Y_{i}$ and $\B Y_{i'}$ are independent for $i\not = i'$. In general, I use capital letters to represent random variables or matrices, bold font for vectors and matrices, and small letters for specific observations. Hence, I denote $\B y_{i}$ as the observed outcome vector for individual $i$. 

\section{Cluster analysis}
One of our most natural abilities is grouping and sorting objects into categories. We simplify our world by recognizing similarities, grouping similar objects, and naming the groups. For example, when children first learn colors, different shades and hues of blue are grouped together as blue. Similarly, golden retrievers, poodles, and terriers are all dogs. These group labels allow us to summarize a large amount of information in a way that can be easily understood. Identifying and labeling groups of similar items is also important when summarizing a large data set so that information can be retrieved more efficiently by providing a concise description of patterns of similarities and differences in the data. 

Cluster analysis is one tool to explore groups within a data set and it has found its way into various scientific disciplines. From mathematics and statistics to biology and genetics to economics and market research, each field uses its own terminology to describe the grouping process. It is called numerical taxonomy in biology, unsupervised pattern recognition in the machine learning literature, and data segmentation in market research, but the problem is same: to find groups of similar objects.
 
The problem is well defined, but the solution is not. There is no agreed upon criteria that determines that one grouping is better than another. In general, most statisticians agree that clusters should be formed by maximizing the similarity within groups and maximizing the dissimilarity between groups. Some believe that there are naturally occurring groups that need to be discovered, but others remark that cluster analysis can always find clusters even if there is no true underlying group structure. However, \textcite{bonner1964} suggests that the user is the ultimate judge of the meaning and value of the clusters. Clusters should not be judged on whether they are true or false but rather in their usefulness. No method works well on every data set to find meaningful groups. For this reason, there is a wealth of clustering algorithms in the literature. 

In this dissertation, the general goal of clustering is to maximize similarity within groups understanding that the analysis is highly dependent on how the user determines individuals to be similar. In this way, clusters are meaningful to the user while maintaining the basic mathematical philosophy of clustering.  I overview two approaches to clustering as they provide a foundation for the methods discussed in this dissertation. 

\section{Nonparametric clustering methods}\label{sec:nonpar}
The first clustering approach focuses on explicitly defining similarity between subjects rather than making assumptions about how the data were generated. The three key ingredients to these methods are the dissimilarity measure, the clustering algorithm, and the number of clusters.

\subsection{Dissimilarity measures}
There are many terms used for a measure to quantify the distinctiveness of a pair of subjects. Metric, distance, dissimilarity, and similarity are all related concepts. A metric is any function, $d$, that satisfies the following three mathematical properties,
\begin{align*}
(i)& \quad d(\B x,\B y)\geq 0 \text{ with equality if and only if }\B x =\B y\\
(ii)& \quad d(\B x,\B y) = d(\B y,\B x)\\
(iii)&\quad d(\B x,\B z) \leq d(\B x,\B y) + d(\B y,\B z)
\end{align*}
for vectors $\B x,\B y,\B z \in\mathbb{R}^{m}$. Some people use the term distance function to describe a function that satisfies the first two conditions but not necessarily the third condition, known as the triangle inequality. Therefore, all metrics are also distance functions by this definition.

The most commonly used metric on Euclidean space is the Minkowski distance of order $p$ defined as
$$d(\B x,\B y) = \|\B x - \B y\|_{p} = \left(\sum^{m}_{j=1}|x_{j}-y_{j}|^{p}\right)^{1/p}$$
where $\B x = (x_{1},...,x_{m})$ and $\B y = (y_{1},...,y_{m})$. The metric equals the Euclidean distance when $p=2$,
$$||\B x - \B y||_{2}=\sqrt{\sum_{j=1}(x_{j}-y_{j})^{2}},$$
and the Manhattan distance when $p=1$,
$$||\B x - \B y||_{1}=\sum_{j=1}|x_{j}-y_{j}|.$$
As $p$ increases to infinity, the distance converges to the Chebyshev distance,
$$||\B x - \B y||_{\infty}=\max_{j}|x_{j}-y_{j}|.$$

Other distance measures are based on variants of the correlation coefficient. These do not satisfy the triangle inequality and are not considered metrics. One popular distance measure based on the Pearson correlation coefficient is defined as
$$d(\B x,\B y) =1 - \frac{(\B x - \bar{x})^{T}(\B y - \bar{y})}{\|\B x - \bar{x}\|_{2}\|\B y - \bar{y}\|_{2}}$$
where $\bar{x}$ represents the average of the elements in $\B x$ and similarly for $\bar{y}$. The codomain of this distance function is $[0,2]$. Perfect negative linear correlation results in the largest distance, perfect positive linear correlation results in the smallest, and no linear correlation is in the middle. The Pearson correlation coefficient can be replaced with its uncentered version to calculate the cosine-angle distance, defined as
$$d(\B x,\B y) =1 - \frac{\B x ^{T}\B y }{\|\B x\|_{2}\|\B y \|_{2}}.$$
When vectors are exactly the same, the angle between the vectors is 0 and the distance equals 0. On the other hand, the maximum value of 2 is attained when two vectors are 180 degrees apart. Both of these distance measures are invariant to scaling such that $d(\B x, a\cdot \B x) = 0$ for $a>0$ and the Pearson correlation distance is additionally invariant to shifts such that $d(\B x, a\cdot \B x +b) = 0$ for $a>0$ and $a,b\in \mathbb{R}$. 

It is also important to note that the distance functions introduced thus far require vectors of equal length and are invariant to permutation of the order of the elements in the vector.

More generally, a similarity function satisfies the following three conditions,
\begin{align*}
(i)&\quad S(\B x,\B y)\geq 0 \\
(ii)&\quad S(\B x,\B y) = S(\B y,\B x)\\
(iii)& \quad S(\B x,\B y)\text{ increases in a monotone fashion as }\B x\text{ and }\B y\text{ are more similar}
\end{align*}
where objects $\B x$ and $\B y$ need not be vectors in Euclidean space. A dissimilarity function is defined in the same way except that the function increases as the objects are more dissimilar. This general definition allows the function to be flexible enough to take into account how one deems two objects to be similar or dissimilar. In many cases the term distance and dissimilarity is used interchangeably if the objects are numeric vectors. The choice of the dissimilarity measure needs to be made after careful study as it can have a dramatic effect on the final clustering results. 

\subsection{Algorithms}
Clustering algorithms aim to optimize a criterion based on the chosen dissimilarity measure. Most clustering algorithms can be classified into two group according to their general search strategies: hierarchical and partitioning algorithms. Hierarchical methods involve constructing a tree of clusters in which the root node is a cluster containing all objects and the leaves are clusters each containing one object. The tree can be constructed in a divisive manner (i.e., top down by recursively dividing groups) or agglomerative manner (i.e., bottom up by recursively combining groups). In order to combine groups, there are different ways of measuring the distance between groups of objects: single-linkage, complete-linkage, average-linkage. This approach should be used when there is a priori scientific knowledge of a hierarchical structure in the data.

Partition methods aim to map objects into $k\geq 2$ disjoint groups by maximizing a criterion without any larger hierarchical structure. Two popular methods include K-means \cite{macqueen1967,hartigan1979} and partitioning around medoids (PAM) \cite{kaufman1990}. 

K-means is one of the simplest unsupervised learning algorithms used to solve the well-known clustering problem. Given a set of vectors $(\B y_{1},\B y_{2},...,\B y_{n})$ where $\B y_{i}\in\mathbb{R}^{m}$ for all $i=1,...,n$, the K-means clustering algorithm aims to partition the $n$ vectors into $K$ sets $\{C_{1},...,C_{K}\}$ so as to minimize the sum of squared Euclidean distance to the assigned cluster centroids denoted as
$$\min_{\{C_{k}\}_{k=1}^{K}}\sum_{k=1}^{K}\sum_{i: \B y_{i}\in C_{k}}\|\B y_{i} -\BS \mu_{k}\|_{2}^{2}$$
where $\BS\mu_{1},...,\BS\mu_{K}$ are the centroids or mean vectors of the $K$ clusters. To find the sets that minimize the criterion, the algorithm starts by strategically or randomly choosing $K$ data vectors as the centroids. The general K-means algorithm proceeds by alternating between two steps:
\begin{itemize}
\item Assignment step: Assign each observation to the cluster with the closest centroid. The $k$th set equals
$$ C_{k} = \{\B y_{i}: \|\B y_{i} - \BS\mu_{j}\|_{2}^{2}\geq \|\B y_{i} - \BS\mu_{k}\|_{2}^{2} \;\;\forall\;\; 1\leq j\leq k,\;\; i=1,...,n\}$$
such that every $\B y_{i}$ is in one and only one set.
\item Update step: Calculate the centroid of each new set,
$$\BS\mu_{k} = \frac{1}{|C_{k}|}\sum_{i: \B y_{i}\in\B C_{k}} \B y_{i}.$$
\end{itemize}
The algorithm continues to iterate until the sets no longer change. Depending on the initial partition, the algorithm is expected to converge to local optima; therefore, completing multiple random starts offers the best way to obtain a global optimum. The K-means algorithm works best when data clusters are about equal in size and shape. The algorithm attempts to find spherical clusters since the K-means algorithm is based on squared Euclidean distance. If the groups are not dense spherical densities with many spurious outliers, the centroid may not be representative of any of the cluster members.

The partitioning around medoids (PAM) algorithm attempts to overcome some of the issues with K-means. The algorithm operates on a user-provided dissimilarity matrix rather than squared Euclidean distance. It is robust to outliers since the medoid or middle vector for each group is selected from the observed data vectors rather than based on mean calculations. To find $K$ sets of vectors to minimize the sum of dissimilarity of the vectors with their respective medoids, $K$ individuals are first randomly chosen as medoids denoted as $\BS \nu_{1},...,\BS \nu_{K}$. The PAM algorithm alternates between the following two steps:
 \begin{itemize}
\item Build step: Assign each observational unit to the closest medoid with respect to the given dissimilarity matrix.
\item Swap step: For each $k=1,...,K$, swap the medoid, $\BS\nu_{k}$, with each non-medoid observation and compute the sum of the dissimilarities of the vectors with their closest medoid. Select the configuration with the smallest sum of dissimilarities.
\end{itemize}
The algorithm continues until there is no change in the medoids. Although it has a longer run time, this building and swapping procedure can return a smaller sum of dissimilarity in contrast to a simple assignment and update procedure in the K-means algorithm. 

\subsection{Choosing the number of components}
A major challenge in clustering is to determine the optimal number of clusters. For the partition methods, a maximum number of clusters is chosen $K_{max}<n$ and the algorithm is run for each value $K=2,3,....,K_{max}$. Then, some take a direct approach to choose the optimal $K$  by optimizing functions of within and between cluster dissimilarity \cite{mulligan1985} or the average silhouette \cite{kaufman1990}. Others take a testing approach by comparing a cluster summary with its expectation under an appropriate null distribution using the Gap statistic \cite{tibshirani2001} or the CLEST approach \cite{dudoit2002}. 

The average silhouette has the advantage of working well with any cluster routine and dissimilarity measure. For each vector $i$, the silhouette $s(i)$ is defined as
$$s(i)=\frac{b(i)-a(i)}{\max\{a(i),b(i)\}}$$
where $a(i)$ is the average dissimilarity of the $i$th vector to all other vectors in its cluster, $A$, $d(i,C)$ is average dissimilarity of the $i$th vector to all vectors in cluster $C$ and $b(i)=\min_{C\not= A} d(i,C)$. The overall average silhouette width $\bar{s}$ is the average of $s(i)$ over all observational units in the data set. The chosen clustering algorithm is run and the overall average silhouette is calculated for each possible value of $K$. Then the optimal $K$ is chosen to maximize the average silhouette while minimizing the within-group dissimilarities in comparison to the between-group dissimilarities. 

\section{Model-based clustering methods}
For a chosen number of clusters, partitioning and hierarchical algorithms divide vectors into non-overlapping, non-empty sets. Each subject is in one and only one cluster. This type of clustering is referred to as hard clustering in contrast to soft or fuzzy clustering where each subject can be in multiple groups to varying degrees. This distinction largely impacts subjects that are near the edge between two groups in terms of the chosen dissimilarity measure. With hard clustering, those subjects are forced in one group. With soft clustering, those subjects on the edge contribute to multiple clusters to varying degrees and the clustering results include association levels that indicate the uncertainty in group membership. Probabilistic model-based methods allow for soft clustering with the possibility of formal inference \cite{fraley2002}. Assuming a probability distribution for the data provides a framework in which to estimate the probability of group membership as well as to estimate parameters and make inferences on the relationship between covariates and group membership. The main model used for clustering is the finite mixture model.

\subsection{Finite mixture model}
In a general finite mixture, the density for a random variable $\B y$ takes the form
$$f(\B y|\BS\theta)=\pi_{1}f_{1}(\B y|\BS\theta_1)+\cdots+\pi_{K}f_{K}(\B y|\BS\theta_K)$$
where $f_k$ and $\BS\theta_k$ are the density and parameters of the $k$th component and $\pi_{k}$ is the probability an observation belongs to the $k$th component ($\pi_{k}>0$; $\sum^{K}_{k}\pi_{k}=1$). The full parameter vector $\BS\theta$ includes the prior probabilities $\pi_k$ and the component parameters $\BS\theta_k$. In many situations, the component densities are assume multivariate Gaussian parameterized by a mean vector $\BS\mu_{k}$ and covariance matrix $\BS\Sigma_{k}$,
$$f_k(\B y| \BS\mu_k,\BS\Sigma_k) = \frac{\exp\left(-\frac{1}{2}(\B y - \BS\mu_k)^T\BS\Sigma_k^{-1}(\B y-\BS\mu_k)\right)}{\sqrt{\det(2\pi \BS\Sigma_k)}}.$$
The Gaussian distribution can be reparameterized to accommodate longitudinal data and to reduce the number of parameters. If the mean outcome values are thought to depend on explanatory variables, the mean vectors are replaced with $\BS\mu_k = \B x \BS \beta_k$ such that $\B x$ is a design matrix based on those variables that impact the mean. The covariance matrix can be simplified by assuming a structure such as independence ($\BS\Sigma_{k}=\sigma^{2}_{k}\B I$), compound symmetry ($\BS\Sigma_{k}=\sigma^{2}_{k}(\rho_{k}\B1\B1^{T}+(1-\rho_{k})\B I)$), or exponential covariance  ($[\B \Sigma_{k}]_{jl} = \sigma_{k}^{2}\exp(-| t_{ij}-t_{il}| / r_{k})$). On the other hand, the covariance structure can be parameterized through the eigenvalue decomposition of the form
$$\BS\Sigma_k=\lambda_k \B D_k\B A_k\B D_k^T$$
where $\B D_k$ is the orthogonal matrix of eigenvectors, $\B A_k$ is a diagonal matrix whose elements are proportional to the eigenvalues, and $\lambda_k$ is a proportional constant \cite{banfield1993}. The mixture model can also be extended to include more complexity through further parameterizations. If baseline factors are thought to influence group membership, the prior group probabilities can be parameterized using a generalized logit model
$$\pi_k(\B w,\BS\gamma)= \frac{\exp(\B w^T \BS \gamma_k) }{\sum_{j=1}^K \exp(\B w^T \BS \gamma_j)}$$
where $\BS\gamma_{K}=\B 0$ and $\B w$ is a design vector based on baseline factors. 

\subsection{Expectation-maximization algorithm}\label{sec:em}
Under the assumption that $\B y = \{\B y_{1},...,\B y_{n}\}$ are independent realizations from the mixture distribution $f(\B y | \BS\theta)$ defined above, the log-likelihood function for the full parameter vector $\BS \theta$ is given by
$$\log L(\BS\theta)=\sum^{n}_{i=1}\log f(\B y_{i}|\BS \theta).$$
The maximum likelihood estimate of $\BS\theta$ is obtained by finding an appropriate root of the score equation, $\partial \log L(\BS\theta)/\partial \BS\theta=\B 0.$ Solutions of this equation corresponding to local maxima can be found iteratively through the expectation-maximization (EM) algorithm \cite{dempster1977}, which was originally presented in the context of an incomplete data problem. In the clustering case, $\B y_{i}$ is assumed to have stemmed from one of the components and the label denoting its originating component is missing. Let $z_{ik}$ equal 1 if $\B y_{i}$ was generated from component $k$ and 0 otherwise. The component label vectors $\B z_{i}=(z_{i1},...,z_{iK})$ for $i=1,..,n$ are assumed to be realizations of random vectors $\B Z_{i}=(Z_{i1},...,Z_{iK})$ that follow a categorical distribution with probabilities $\pi_{1},...,\pi_{K}$. The complete data log-likelihood function, based on these component labels and the observed data $\B y$, is given by
$$\log L_{c}(\BS \theta) =  \sum^{K}_{k=1}\sum^{n}_{i=1}z_{ik}\left[\log \pi_{k}+\log f_{k}(\B y_{i}|\BS\theta_k)\right].$$

The EM algorithm treats the $z_{ik}$ as missing data and iteratively imputes the missing values and then  estimates the parameters. The expectation step (E-step) involves taking the conditional expectation of the complete log-likelihood given the data and the current value of the parameter estimates. Since the complete log-likelihood is linear in the unobserved data $z_{ik}$, this step involves calculating the conditional expectation of $Z_{ik}$ given the observed data. On the $t$th iteration, the expectation is
\begin{align*}
E_{\BS\theta^{(t-1)}}( Z_{ik}|\B y) &= P_{\BS\theta^{(t-1)}}( Z_{ik}|\B y)\\
&=\alpha_{ik}^{(t)}
\end{align*}
using the estimates of the parameters from the previous iteration. The quantity $\alpha_{ik}^{(t)}$ is the posterior probability that the $i$th individual belongs to component $k$, written as
$$\alpha_{ik}^{(t)}=\pi_{k}^{(t-1)}f_{k}(\B y_{i}|\BS\theta^{(t-1)}_{k})/\sum_{j=1}^{K}\pi_{j}^{(t-1)}f_{j}(\B y_{i}|\BS \theta_{j}^{(t-1)})$$
for $i=1,...,n$ and $k=1,...,K$.

In the maximization step (M-step), the parameter estimates for the prior probabilities and component parameters are updated by maximizing the conditional expectation of the complete-data log-likelihood from the E-step. On the $t$th iteration, the updated estimates of the prior probabilities $\pi_{k}^{(t)}$ are calculated independently of the estimates of the component parameters $\BS\theta_{k}^{(t)}$. The prior probability estimates equal 
$$\pi_{k}^{(t)}=n^{-1}\sum^{n}_{i=1}\alpha_{ik}^{(t)}.$$
If the generalized logit function is used for the prior probabilities to include dependence on baseline covariates such that $\pi_{k}=\pi_{k}(\B w_{i},\BS\gamma) = \frac{\exp(\B w_{i}^T \BS \gamma_k) }{\sum_{j=1}^K \exp(\B w_{i}^T \BS \gamma_j)}$, then the updated value of $\BS\gamma^{(t)}=(\BS\gamma_{1}^{(t)},...,\BS\gamma_{K-1}^{(t)})$ is an appropriate root of 
$$\sum^{K}_{k=1}\sum^{n}_{i=1}\alpha_{ik}^{(t)}\partial\log \pi_{k}(\B w_{i},\BS\gamma)/\partial \BS\gamma = \B 0$$
calculated using a numerical optimization routine. The component parameters are updated by finding an appropriate root of 
$$\sum^{n}_{i=1}\alpha_{ik}^{(t)}\partial\log f_{k}(\B y_{i}|\BS\theta_k) /\partial \BS\theta_{k} =\B 0$$
at the $t$th iteration for each component $k=1,...,K$. Closed-form solutions to this equation exist for many Gaussian mixture models.  

The E-step and M-step are alternated repeatedly until convergence. \Textcite{dempster1977} showed the incomplete likelihood increases monotonically at each iteration of this algorithm. Hence, the EM algorithm guarantees convergence to a local maximum; global convergence may be attained by running the algorithm multiple times randomly assigning individuals to initial components and using the estimates associated with the one with the highest log-likelihood. For more details about the EM algorithm, see \textcite{mclachlan1997}. 

When the algorithm has converged, the values of the parameters and posterior probabilities from the last iteration are taken as the final estimates. The component estimates characterize the mean and variances within the clusters and the posterior probabilities provide the strength of group association for the soft clustering. This clustering can be translated into a hard clustering by choosing the group label that maximizes the posterior probability, $\arg\max_k \alpha_{ik}$, for each subject. Then, one measure of uncertainty in the group labels is $(1-\max_k\alpha_{ik}).$ A variant of the EM algorithm called the classification EM \cite{celeux1992}, in which the posterior probabilities are converted to group indicators before the M-step, is equivalent to the K-means algorithm for a Gaussian mixture with $\Sigma_{k} =\sigma \B I$. For further discussion about other estimation approaches besides maximum likelihood via the EM algorithm, see \textcite{mclachlan2000}.

\subsection{Estimation issues}
Although an estimation tool exists, there are potential issues of parameter identifiability with mixture models. \Textcite{fruhwirth2006} distinguished between three types of nonidentifiability: invariance to  relabeling of components, potential overfitting, and nonidentifiability due to the family of component distribution and the covariate design matrix. The first two issues are resolved through constraints such as $\BS\theta_{k}\not=\BS\theta_{k'}$ for all $k,k'=1,...,K$, $k\not=k'$. The last concern is often solved by assuming Gaussian components since finite mixtures of multivariate Gaussians are identifiable \cite{teicher1963,yakowitz1968}. However, \textcite{hennig2000} suggested that the introduction of a regression structure to a Gaussian mixture requires a full rank design matrix as well as a rich covariate domain for regression parameters to be identifiable. On the other hand, prior probabilities parameters from a generalized logit are identifiable by setting the parameters of one component to zero such as $\BS\gamma_{K}=\BS0$ \cite{jiang1999}.

Besides identifiability, there are other known issues with finite mixture models. McLachlan and Peel \cite{mclachlan2000} noted that the sample size must be quite large for asymptotic theory to accurately describe the finite sampling properties of the estimator. Also, when component variances are allowed to vary between components, the mixture likelihood function is unbounded, and each observation gives rise to a singularity on the boundary of the parameter space \cite{day1969,kiefer1956}. However, Kiefer \cite{kiefer1978} outlined theory that guarantees that there exists a particular local maximizer of the mixture likelihood function that is consistent, efficient, and asymptotically normal if the mixture is not overfit. To avoid issues of singularities and spurious local modes in the EM algorithm, Hathaway \cite{hathaway1985} considered constrained maximum likelihood estimation for multivariate Gaussian mixtures based on the following constraint on the smallest eigenvalue of the matrix $\BS\Sigma_{h}\BS\Sigma_{j}^{-1}$, denoted as $\lambda_{\text{min}}(\BS\Sigma_{h}\BS\Sigma_{j}^{-1})$,
$$\min_{1\leq h\not=j\leq K}\lambda_{\text{min}}(\BS\Sigma_{h}\BS\Sigma_{j}^{-1})\geq c > 0$$ for some positive constant $c\in[0,1]$ to ensure a global maximizer. 

\subsection{Choosing the number of components}
Just as with partitioning methods, it is difficult to choose the number of mixture components and the problem has not been completely resolved \cite{mclachlan2000}. A natural approach is to test the hypothesis that $K=K_{0}$ against the alternative that $K=K_{1}$ with a likelihood ratio test. However, the regularity conditions do not hold for the test statistic to have its usual asymptotic null distribution. This problem has been considered by many authors (see \textcite{mclachlan2000} for a summary) who have suggested modified test statistics \cite{wolfe1971}, approximate null distributions \cite{lo2001}, and resampling methods to empirically estimate the null distribution \cite{mclachlan1987}.

Another approach is to use information theoretic methods to select a model. A set of candidate models are compared by measuring the information lost when the model is used to approximate the true reality. The model with the lowest information loss and lowest information criterion is chosen as the best model. Two popular information criterion are the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). 

The information loss can be described in term of the Kullback-Leibler (KL) information \cite{kullback1951} of the true distribution with respect to the fitted model. If $g(\B y)$ is the true density, the KL information of $g(\B y)$ with respect to the estimated model $f(\B y|\hat{\BS \theta})$ is
$$\int \log[g(\B y)/f(\B y|\hat{\BS\theta})]g(\B y)d\B y=\int \log[g(\B y)]g(\B y)d\B y-\int \log[f(\B y|\hat{\BS\theta})]g(\B y)d\B y.$$
Since the first term does not involve the model, only the second term is relevant. When comparing models, this integral can be estimated using the empirical distribution, which places equal mass $1/n$ at each observation $\B y_{i}$,
\begin{align*}
\int \log[f(\B y|\hat{\BS\theta})]g(\B y)d\B y &\approx \frac{1}{n} \sum^{n}_{i=1}\log f(\B y_{i}|\hat{\BS\theta})\\
&=\frac{1}{n} \log L(\hat{\BS\theta}).
\end{align*}
This is a biased estimate of the expected log density. Akaike \cite{akaike1973,akaike1974} showed that the bias is asymptotically equal to $d$, where $d$ is equal to the number of parameters in the model. Thus, he proposed choosing a model that minimizes
$$AIC = -2\log L(\hat{\BS\theta}) + 2d.$$
The AIC criterion is often used to choose the number of components in a mixture model. However, many authors have observed that it is inconsistent \cite{koehler1988} and tends to overestimate the number of components in a mixture \cite{soromenho1994,celeux1996}.

The BIC was originally derived in a Bayesian framework as an approximation to the integrated likelihood  \cite{schwarz1978}, which is used in Bayes factors to compare two models. For further information about Bayes factors, see \textcite{kass1995}. The BIC is calculated as
$$BIC= -2 \log L(\hat{\BS\theta}) + d\log n $$
where $d$ is the number of parameters in the model and $n$ is the sample size. The criterion can be justified in a frequentist framework, but the regularity conditions breakdown for mixture models. Nevertheless, \textcite{fraley1998} noted considerable theoretical and practical evidence to support the use for clustering. Unlike the AIC, the BIC has been shown not to asymptotically underestimate the true number of components \cite{leroux1992}. In the context of nonparametric density estimation,  \textcite{roeder1997} showed that the density estimate chosen with the BIC is consistent. However, if the component densities assumptions do not hold, \textcite{biernacki2000} found that it tends to fit too many components. 

To overcome the shortcomings of the BIC for the mixture model, \textcite{biernacki2000} proposed the ICL-BIC based on the integrated complete likelihood for a mixture with $K$ components,
$$ICL-BIC = -2\log L(\hat{\BS\theta}) + d\log n - 2\sum^{K}_{k=1}\sum^{n}_{i=1}\hat{z}_{ik}\log \alpha_{ik}$$
where $\alpha_{ik}$ is the posterior probability of individual $i$ belonging to component $k$, $\hat{z}_{ik}$ is the maximum a posteriori estimate of group membership, $d$ is the total number of parameters, and $n$ is the sample size. The last term is an entropy term based on the strength of the group classification. If the components of the mixture are well separated, then the entropy will be close to zero, but when the components are largely overlapping, the entropy is large and acts as a penalty for having poorly separated clusters. While this criteria has an additional penalty, it is not used as frequently as BIC in practice.

When comparing information criterion between competing models, large differences correspond to strong evidence for one model over the other. \Textcite{kass1995} provide standard conventions for calibrating differences in BIC which can guide the final choice of models. When the differences are small, the smallest number of components that fit the data is chosen for the sake of parsimony.

\section{Thesis outline}
This dissertation is organized in the following manner. In Chapter \ref{chap:misspecify}, I develop a series of simulations to study the impact of misspecifying the covariance structure in finite mixture models for longitudinal data and present finite-sample and asymptotic bias results. In Chapter \ref{chap:motivate}, I motivate the rest of the dissertation by discussing the general clustering methodology for longitudinal data and illustrate how these methods do not cluster based on the shape of change over time despite the natural tendency to interpret the results using shape terminology. Chapter \ref{chap:methods} presents three proposed methods to cluster based solely on the shape and compares them to standard methods through a simulation study in Chapter \ref{chap:sim}. Chapter \ref{chap:data} applies these methods to longitudinal growth data of children. In Chapter \ref{chap:concl}, I overview the results of the dissertation, discuss the contributions, and present areas of future work.