Chapter5.tex

\chapter{Simulations}
\label{chap:sim}
In this chapter, I compare the three proposed methods presented in the last chapter with competing clustering methods in a simulation study. The goal of this study is to examine the performance of these methods when applied to data of individuals with the same shape trajectory at a variety of levels. No method works well in all situations; therefore, the data are generated under a variety of conditions to test the methods on their robustness to model misspecification, noise, and group overlap. The ability to detect and create homogenous shape groups indicates good performance. The outcomes of interest are the data-driven number of groups, misclassification with reference to the generating shape groups, and the parameter estimates for the cluster mean shape and baseline factors. 

First, I describe the data-generating process used throughout the simulation study. Then, I discuss the specific details of the simulation implementation in terms of the variety of data conditions, methods used, and the simulation outcomes of interest. Lastly, I present the results and compare the performance of the methods.

\section{Data-generating process}
When deciding upon the data-generating process, it is important to choose a simple example that addresses the goals of the study and can generalize to other situations. This simulation study is based on a population with three trajectory shapes and three outcome levels. Since complex functions can be locally approximated with linear functions, the three shapes used in the simulation are straight lines with different slopes---positive, negative, and zero. We restrict some shapes to particular levels to induce a relationship between level and shape. The five mean functions for generating data with three shapes at different levels are  
\begin{align*}
\mu_{1}(t) &= -1 - t &&\text{ (negative slope, low level)}\\
\mu_{2}(t) &= 11 - t&&\text{ (negative slope, high level)}\\
\mu_{3}(t) &= 0&&\text{ (zero slope, middle level)}\\
\mu_{4}(t) &= -11 + t&&\text{ (positive slope, low level)}\\
\mu_{5}(t) &= 1 + t&&\text{ (positive slope, high level)}
\end{align*} 
Individuals follow these mean patterns with varying probabilities that depend on two baseline factors. The first factor $w_{1}$ impacts the shape and $w_{2}$ impacts the level. Individuals are randomly assigned values of these two binary factors with independent simulated tosses of a fair coin. 

Let $S$ be a categorical random variable that indicates the shape/slope group such that $S=1,2,3$ refers to the negative, zero, and positive slope groups, respectively. Conditional on the baseline factors, the probability of being in the $k$th shape group is
$$P(S=k |w_{1}) = \frac{\exp(\gamma_{0k}+\gamma_{1k}w_{1})}{\sum^{3}_{l=1} \exp(\gamma_{0l}+\gamma_{1l}w_{1})}$$
for $k=1,2,3$ where $\gamma_{01}=2,\gamma_{11} = -4,\gamma_{02}=1.5,\gamma_{12}=-2,\gamma_{03}=\gamma_{13} = 0$ and $w_{1}\in\{0,1\}$. Since the value of $w_{1}$ is determined by a fair coin toss, each shape group has about an equal probability, marginally. 

The second factor impacts the level. Let $L$ be a categorial random variable that indicates level group such that $L=1,2,3$ refers to the low, middle, and high group, respectively. Only those who follow the constant function have middle level outcome measures. For those in either the negative or positive slope groups, the chance of the high or low level is
\begin{align*}
P(L=k|S=1 \text{ or } S=3,w_{2}) & = \frac{\exp(\zeta_{0k}+\zeta_{1k} w_{2})}{\sum_{l\in\{1,3\} }\exp(\zeta_{0l}+\zeta_{1l}w_{2})}
\end{align*}
for $k=1,3$ and $w_{2}\in\{0,1\}$ where $\zeta_{01}=0,\zeta_{11}=0,\zeta_{03}=-3,\zeta_{13}=6$. Again, each outcome level group has about the same marginal probability.

For each individual, the chosen mean function is evaluated at five equidistant observation times $t=1,3.25,5.5,7.75,10$ that span the period 1 to 10 units and random noise is added to induce variability. The random noise is made up of two independent ingredients: individual-specific level perturbation and time-specific Gaussian measurement error. For individual $i$ ($i=1,...,n$)  following the $l$th mean function, the observed outcome at the $j$th observation time ($j=1,..,5$) is a realization of
$$y_{ij} = \lambda_{i}+\mu_{l}(t_{j})+\epsilon_{ij}\quad\text{where}\quad \epsilon_{ij}\stackrel{iid}{\sim} N(0,\sigma_{\epsilon}^{2}), \lambda_{i}\stackrel{iid}{\sim} F_{\lambda}(0,\sigma_{\lambda}^{2})$$
where $\sigma_{\epsilon}$  is the standard deviation of the measurement error, $\sigma_{\lambda}$ is the standard deviation of the level perturbation, and $F_{\lambda}$ is the probability distribution of the level perturbation. This imposes an exchangeable correlation structure within an individual's outcome vector conditional on the mean function. The distribution of the level perturbation as well as the standard deviations of the random noise are values that vary in the simulation to create different conditions on which to test the clustering methods. The values are chosen to alter the overlap of the groups and the signal to noise ratio and are discussed in the next section.

\section{Implementation}
In this section, I describe the random noise conditions, overview the clustering methods included in this study, and discuss the outcomes of interest used to compare the methods in the simulation. 

\subsection{Simulation conditions}
The first ingredient of the random noise is the level perturbation.  To test the sensitivity of the methods to distributional assumptions, I use the uniform distribution and the Gaussian distribution to generate the level perturbations. To create conditions with differing amount of overlap between groups, I let $\sigma_{\lambda}=2$ or $3$ such that the vertical shift between the mean functions with the same shape is about 4 or 6 standard deviations. The second ingredient is the measurement error, the magnitude of which influences the signal to noise ratio. I let $\sigma_{\epsilon}=0.5$ or $2$ to create two extreme conditions, strong and weak signals relative to noise. The eight possible combinations represent the conditions of the data-generating process in the simulation study. 

\subsection{Clustering methods}
For each condition, I generate a data set of $n=500$ individuals using the process described above and apply each method to cluster the data into $K=2,3,4,$ and $5$ groups. Each method produces different output making direct comparisons difficult. 

Partition methods produce group labels for each individual while model-based methods result in posterior probability estimates for each individual and group. Therefore, I translate the estimated probabilities into group labels by placing individuals into the group that has the highest posterior probability.

Besides grouping results, the group representatives are in different formats. The K-means algorithm results in a mean vector within each group, the interpretation of which depends on the input vectors. On the other hand, I model the mean with spline basis functions in the model-based clustering methods. The model output includes the estimated coefficients for the group-specific mean functions. To make the group representatives from the partition and model-based methods comparable, I transform all of the output into coefficients of a differentiated B-spline function which describes the shape of the curves. I project the data vectors onto a common spline basis, if not already in that form, differentiate the spline functions, and rearrange terms to calculate the coefficients of the spline function one order lower. 

The common spline used throughout the simulation study is a quadratic B-spline with no internal knots and boundary knots at 0.5 and 10.5. A quadratic B-spline basis without internal knots is equivalent to a Bernstein polynomial basis \cite{lorentz1953} which generalizes the standard polynomial basis; in general, I advocate using a B-spline basis that provides flexibility and is computationally stable.  As mentioned above, this function is differentiated and the coefficients are transformed in order to get the coefficients of the derivative B-spline function for final method comparisons. See Section \ref{sec:bsplines} for more details on derivatives of B-spline functions.

For every clustering method, I estimate the relationship between the two baseline factors, $w_{1}$ and $w_{2}$, and the group membership. This estimation procedure differs between the methods as it may simultaneously occur while estimating group parameters or after the clustering algorithm is complete. I now overview the clustering methods used in this study.

\subsubsection{K-means}
The K-means algorithm is a general clustering procedure for numerical vectors. $K$ groups are determined through an iterative assignment and update process to minimize the within-group sum of squares distance to the group mean vectors. Note that the objective function indirectly imposes a spherical shape on the groups in the vector space.

To choose the optimal value of $K$ in this simulation study, I use the silhouette measure, which measures the compactness of clusters, for all applications of the K-means algorithm. See Section \ref{sec:nonpar} for more details. 

To estimate the relationship between baseline factors and group membership produced from the K-means algorithm, the groups labels are taken as known and become the outcomes for a multinomial logistic regression with the observed baseline factors as explanatory variables. Let $c_{i}$ be the group label produced by the clustering algorithm and $w_{i1}, w_{i2}$ be the baseline factors for individual $i$, then the model assumes that
$$\log\frac{P(c_{i}=k)}{P(c_{i}=K)} = \gamma_{0k}+\gamma_{1k}w_{i1}+\gamma_{2k}w_{i2}$$
for $k=1,..,K-1$. Estimates for the coefficients $\BS\gamma = (\gamma_{01},\gamma_{11},\gamma_{21},...,\gamma_{0\;K-1},\gamma_{1\;K-1},\gamma_{2\;K-1})$ are calculated via maximum likelihood estimation as long as the number of individuals in each group is adequately large.

\subsubsection{K-means on original data}
The K-means algorithm applied to the 500 original data vectors results in groups based on both the level and shape. Both characteristics can play a role since this method is based on squared Euclidean distance of the original outcome measures. In order to compare the estimated group representatives to those of the following methods, the $K$ group mean vectors are projected onto the common B-spline basis and the coefficients are then transformed into derivative coefficients.  

\subsubsection{K-means on spline coefficients}
To focus on the individual trends without noise, the original data vectors are projected onto the common B-spline basis prior to clustering. This process in effect smoothes each individual trajectory into a quadratic function. The vector of basis coefficients for each individual are used as the input for the K-means algorithm. This technique originates in the functional data analysis literature \cite{ramsay2005} and has been used as a way to decrease the dimension of the data prior to clustering. 

\subsubsection{K-means on derivative spline coefficients}
This is the first proposed method of this thesis. A slight adjustment to the previous method allows the clustering to be based on the derivative function, which describes the shape, rather than the original function. Before the algorithm is applied, the spline coefficients are transformed into coefficients of the derivative function written as a B-spline function of one lower order. 

\subsubsection{K-means on difference quotients}
Another way of clustering on the basis of derivatives is to calculate difference quotients by taking the pair-wise differences within the ordered data vector and divide by the time between observations,
$$\Delta_{ij} = \frac{y_{ij+1}-y_{ij}}{t_{j+1}-t_{j}}$$
for $i=1,...,n$ and $j=1,...,4$. The vector of difference quotients, $\Delta_{i}=(\Delta_{i1},...,\Delta_{i4}),$ is used as the input for the K-means algorithm.

\subsubsection{PAM with correlation measure}
Partitioning around medoids (PAM) is a variant of the K-means algorithm. Rather than focusing on the distance to the  group mean or centroid, this method uses the middle data vector or medoid as the representative of the group, which results in a more robust procedure less sensitive to outliers. The algorithm is general enough to allow for other measures of dissimilarity beyond squared Euclidean distance. For this study, the chosen dissimilarity measure is a linear transformation of the Pearson's correlation coefficient, which measures the linear association between two vectors and is believed to detect similarity in the shape of longitudinal data \cite{chouakria2007,  eisen1998,chiou2008}. The dissimilarity between two vectors is defined as
$$d(\B y_{1},\B y_{2}) = 1-Cor(\B y_{1},\B y_{2})$$
where $Cor(\B y_{1},\B y_{2})$ is defined in Section \ref{sec:nonpar}. Perfect positive linear correlation results in a dissimilarity of 0 and perfect negative correlation results in a value of 2. If a vector is an additive shift of another vector, they are deemed to be similar. This is one desired property for a clustering method in this thesis. However, the same goes for multiplicative transformations of data vectors which change the shape of the trajectory. Therefore, this measure cannot detect magnitude differences in the slope of trend lines and cannot distinguish between a horizontal line and one with a positive or negative slope.

To choose the number of groups and estimate the relationship between baseline factors and group membership, the same procedures as the K-means algorithm are used.

\subsubsection{Finite mixture model}
In contrast to the partition methods mentioned above, using a probability model allows for uncertainty in the model parameters and group membership. The data are fit to the model using the EM algorithm in order to maximize the likelihood function (see Section \ref{sec:em} for more details). Rather than being separated into $K$ groups, mixture models allow individuals to contribute to every group to some degrees. For example, if two groups are not well separated, then an individual in the overlap of the groups could contribute equally to the estimation of both groups' parameters. 

For this simulation study, I generally assume the $K$ components of the finite mixture model are multivariate Gaussian with a component-specific mean function that can be expressed as a quadratic B-spline with boundary knots at 0.5 and 10.5. The estimated coefficients for the mean parameters are then transformed into the coefficients for the derivative spline function for the final comparison. Additionally, the prior probabilities of the components are parameterized using a generalized logit function such that the logit is equal to a linear combination of the two baseline factors. For individual $i$, 
$$\log\frac{P(c_{i}=k)}{P(c_{i}=K)} = \log\frac{\pi_{k}(w_{i1},w_{i2},\BS\gamma)}{\pi_{K}(w_{i1},w_{i2},\BS\gamma)} =  \gamma_{0k}+\gamma_{1k}w_{i1}+\gamma_{2k}w_{i2}.$$
The parameters, $\BS\gamma=(\gamma_{01},...,\gamma_{2\; K-1})$ are estimated simultaneously with the mean and covariance parameters. 
The number of components or groups is selected by minimizing the model selection criteria, BIC, over model fits with $K=2,3,4,5$ (see Chapter 1 for more details). 

\subsubsection{Finite mixture model with independence}
For the standard model fit to the original data, the covariance structure is assumed to be independent such that the variance is constant over time but differs between components. Therefore, the only covariance parameters in the model are the component-specific variances.

\subsubsection{Finite mixture model with exponential correlation}
As discussed in Chapter \ref{chap:misspecify}, assuming independence for longitudinal data when there is inherent dependence between repeated measures can result in bias and potentially misleading clustering results. Therefore, we use an exponential correlation structure that maintains the stationarity of the covariance structure but allows the correlation to decay as the time between observations increases. The number of covariance parameters increases to two per component: variance and range of dependence. 

\subsubsection{Multilayer mixture model}
This is the second proposed method of this thesis. In generalizing the standard finite mixture model to allow shape groups to have different level groups, I assume that each shape component of the mixture model is a mixture of level components assumed to be multivariate Gaussian. In order to minimize the number of parameters to be estimated, I assume independence between repeated measures. For this two-layered model, the prior probabilities for the shape components, rather than those of the level components, are parameterized to relate baseline factors to group membership. In addition to specifying the number of shape components, the number of level components within each shape component need to be set prior to estimating the parameters. It is difficult to choose both of these values in an automatic way; therefore, we fix the number of shape components to be $K=3$ since in a simulation it is known a priori. Then we use the BIC to select the number of level components within the shape components \cite{li2005}. 

\subsubsection{Vertically shifted independence mixture model}
This is the third proposed method of this thesis. Rather than trying to explicitly model the level, I remove the level prior to fitting the model by subtracting the individual-specific mean. Before the model is fit to the data, the trajectories are vertically shifted so they lie in the same range. The hope is that shape drives the group development once the level is removed. As mentioned in Chapter \ref{chap:methods}, subtracting the mean impacts the covariance structure and can result in negative correlation between repeated measures, which is hard to model with a legitimate correlation function. In this simulation, the correlation structure of the original data for this data-generating process is exchangeable and observation times are fixed. Therefore, we know that the correlation of the transformed errors is exchangeable with magnitude $\frac{-1}{4}$. The true correlation is on the boundary of the parameter space; therefore, I use the independence correlation structure to model the leftover correlation after removing the mean. The same estimation and modeling procedures as the standard finite mixture model apply.

\subsubsection{Vertically shifted exponential mixture model}
Rather than assuming independence for the transformed vector, I use the exponential correlation structure to model the error structure after removing the mean.

\subsection{Outcomes of interest}
For each method and condition, I generate a data set of $n=500$ individuals and calculate the following outcomes of interest. This is repeated so that we get 500 unique data sets under each condition on which we can apply each method and summarize the results.

\subsubsection{Number of groups}
For each replication, I use the method-specific procedure of choosing the final data-driven number of groups and record the final value of $K$ for each data set. Note that this $K$ is not chosen for the multilayer mixture as the number of shape groups is fixed to be equal to three.

\subsubsection{Misclassification rate}
To detect whether the method discovers the underlying shape structure, I fix $K=3$ and compare the cluster memberships labels to the true shape group membership using a contingency table. The cluster label columns of the contingency table are reordered such that the trace of the 3 by 3 inner matrix is maximized to match the clusters to the true shape groups. The sum of the off-diagonal elements of the newly permuted matrix represents the number of misclassified individuals. This value is divided by the sample size to get the rate. The best method results in zero misclassifications. To summarize the replications, I present the mean misclassification rate of all 500 replications.

\subsubsection{Adjusted Rand Index}
Another way to measure the similarity between the produced clusters and the true shape clusters is the Adjusted Rand Index \cite{hubert1985,mulligan1985}. This measure assesses the agreement of two partitions. In our case, the true shape groups form the reference partition and I compare it to the cluster partition produced by the method when $K=3$.
 
Given a set of $n$ objects suppose $U = \{u_{1},...,u_{K}\}$ and $V=\{v_{1},...,v_{K'}\}$ represent two different partitions of the objects such that $\cup^{K}_{i=1}u_{i}  = \cup^{K'}_{j=1}v_{j}$ and $u_{i}\cap u_{i'} = \emptyset = v_{j}\cap v_{j'}$ for $1\leq i\not=i'\leq K$ and $1\leq j\not= j'\leq K'$. The amount of overlap between the two partitions can be written in the form of a contingency table with elements $n_{ij}$ denoting the number of objects that are common to groups $u_{i}$ and $v_{j}$ and $n_{i\cdot}$ and $n_{\cdot j}$ denoting row and column sums, respectively. Let $a$ be the number of pairs of objects that are placed in the same group in $U$ and in the same cluster in $V$, $b$ be the number of pairs of objects in the same group in $U$ but not in the same cluster in $V$, $c$ be the number of pairs of objects in the same cluster in $V$ but not in the same group in $U$, and $d$ be the number of pairs of objects in different groups and different clusters in both partitions. The quantities $a$ and $d$ are interpreted as agreements and $b$ and $c$ as disagreements. The Rand Index \cite{rand1971} equals $\frac{a+d}{a+b+c+d}=\frac{a+d}{{n \choose 2}}$.  The index can yield a value between 0 and 1, where 1 indicates that the clusters are exactly the same.

One issue with the Rand Index is that the expected value of the index does not take on a constant value and thus it is hard to evaluate raw indices. \Textcite{hubert1985} suggested using the generalized hypergeometric distribution for the null model when the partitions are picked at random. This leads to the Adjusted Rand Index, which is of the form $\frac{\text{Index} - \text{Expected Index}}{\text{Max Index}-\text{Expected Index}}$. The index is bounded above by 1 and takes on the value 0 when the index equals its expected value. In the contingency table notation, the Adjusted Rand Index can be simplified to
$$\frac{\sum_{ij}{n_{ij} \choose 2} - \sum_{i} {n_{i\cdot} \choose 2} \sum_{i} {n_{\cdot j} \choose 2} / {n \choose 2}}{\frac{1}{2}\left[\sum_{i}{n_{i\cdot} \choose 2}+\sum_{j}{n_{\cdot j} \choose 2}\right] - \sum_{i} {n_{i\cdot} \choose 2} \sum_{i} {n_{\cdot j} \choose 2} / {n \choose 2}}$$
A method that performs well has an Adjusted Rand Index value close to 1. To summarize the replications, I present the mean Adjusted Rand Index for the 500 replications.

\subsubsection{Mean parameters}
As mentioned in the brief descriptions of the methods used in the simulation, the format for the group representative differs between each method. Six of the ten methods use the B-spline basis to model the mean as a function and one method explicitly uses the derivative of the B-spline function. It is not possible to map a derivative function back to the original function. Therefore, I map the clustering representatives when $K=3$ to the derivative B-spline function so that I can compare the shape of estimated representatives within each group. For reference, the true mean curves are projected onto the B-spline basis and then differentiated. 

Just as with the misclassification rate, it is necessary to match the group labels to the shape groups of the true generating distribution. This is completed by permuting the columns until the trace of the 3 by 3 inner matrix of the contingency table between the true shape labels and the clustering group labels is maximized. Once the groups are matched to the truth, the mean squared error (MSE) is estimated for each derivative coefficient using the 500 replications.

\subsubsection{Baseline factors}
Once groups are discovered, the natural question is to ask who are in the groups and how else do they differ? For each condition, method, and replication, the coefficients of the generalized logit $\{\gamma_{k0},\gamma_{k1},\gamma_{k2}\}$ for $k=1,...,K-1$ are estimated either simultaneously in the model setting or a posteriori for the partition methods when $K=3$. To summarize the estimates from all of the replications, I plot the density of the values and compare them to the true generating values. The clustering is unstable if the density plot is multimodal.

\section{Results}
Tables \ref{tab:freq1}, \ref{tab:freq2}, and \ref{tab:freq3} summarize the simulations in terms of the chosen number of groups, average misclassification rate, and average Adjusted Rand Index. It is clear from these tables that the three standard methods---the K-means algorithm applied to original vectors, K-means on spline coefficients, and the independent Gaussian mixture model---do not select three groups as the optimal number of groups (Table \ref{tab:freq1}). For almost all the conditions, the K-means algorithm on the original data obscures any shape differences in the data by vertically splitting the individuals in half: a high and a low group. If the noise in the data is smoothed out prior to clustering by applying K-means to the B-spline coefficients, five groups is favored except when high variance in the level results in overlap between groups and the algorithm selects two groups. The BIC consistently chooses five groups for the independent mixture model. It recognizes the five mean functions but no similarities in shape or level. All three of these methods also do not perform well when forced to have the true number of shape groups, $K=3$. Only about 50-60\% of the data is correctly specified in terms of the generating shape groups. Small average values for the Adjusted Rand Index confirm that groups from these methods do not match the shape groups.

If the correlation structure of the finite mixture model is generalized to exponential correlation, there is more variability in the number of groups chosen. When the time-specific measurement error is small ($\sigma_{\epsilon}=0.5$), the correct number of groups is chosen about 20\% of the time (Table \ref{tab:freq1}). Under these same circumstances, there is perfect alignment with the generating shape groups when the method is forced to created three groups. However, when the magnitude of the noise is large ($\sigma_{\epsilon}=2$), three is never chosen as the optimal number of groups and the shape groups are not discovered when $K=3$. So even though this method does not directly cluster based on shape, having a general enough correlation structure allows the method to pick up the correct groups when the number of group is known and the measurement error is small. However, it does not choose the correct number of groups in the majority of simulations. 

Of the established methods that are intended to group on shape, K-means on difference quotients picks the correct number of groups only if the magnitude of the measurement error is small relative to the strength of the signal (Table \ref{tab:freq2}). If the variability around the individual mean is large, the method chooses two groups and misclassifies about 38\% of the individuals when forced to have three groups. Using the correlation dissimilarity measure with the PAM algorithm gives slightly better results with only 24-27\% misclassification, but it never chooses three groups.

Now, I compare the three proposed methods (Table \ref{tab:freq3}). Using K-means with the derivative spline coefficients is a slight improvement over the difference quotients such that the true number of groups is chosen more frequently even under conditions of higher measurement error. However, this method did not decrease the misclassification rate or increase the Adjusted Rand Index under those same conditions. 

The multilayer mixture model is successful in correctly classifying individuals into shape groups when the level components are well separated ($\sigma_{\epsilon}=0.5, \sigma_{\lambda}=2$). The method seems to be robust to some misspecification as the results do not degrade when the level is generated by the uniform distribution in comparison with the assumed Gaussian model. However, if the distribution of the level is not distinctly bimodal, the misclassification rate is quite high.

Lastly, the method that prevails amongst the competition is the vertically shifted mixture model. In this case, either correlation assumption works well. For every condition, the method chose three groups as the optimal partition at least 99\% of the time and when forced to have $K=3$, the method discovered the shape groups with little misclassification. Only when the measurement error is large ($\sigma_{\epsilon}=2$) did the method misclassify 5\% (about 25 individuals) in terms of shape. 

% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Feb 27 13:38:21 2013
\begin{table}[ht]
\centering
\begin{tabular}{ccc|cccccc}
  \thickhline  $F_{\lambda}$&$\sigma_{\epsilon}$&$\sigma_{\lambda}$&$K=2$&$K=3$&$K=4$&$K=5$&MR&ARI\\ \hline\multicolumn{9}{c}{\textbf{K-means on Original Data}}\\ Uniform & 0.50 &   2 &   0 &   0 &   0 & 500 & 0.42 & 0.25 \\ 
  Uniform & 2.00 &   2 & 492 &   0 &   0 &   8 & 0.42 & 0.25 \\ 
  Uniform & 0.50 &   3 & 500 &   0 &   0 &   0 & 0.50 & 0.09 \\ 
  Uniform & 2.00 &   3 & 500 &   0 &   0 &   0 & 0.51 & 0.09 \\ 
  Gaussian & 0.50 &   2 &   0 &   0 &   0 & 500 & 0.38 & 0.33 \\ 
  Gaussian & 2.00 &   2 & 273 &   0 &   0 & 227 & 0.39 & 0.29 \\ 
  Gaussian & 0.50 &   3 & 412 &   0 &   0 &  88 & 0.46 & 0.16 \\ 
  Gaussian & 2.00 &   3 & 500 &   0 &   0 &   0 & 0.47 & 0.14 \\ 
   \\ \multicolumn{9}{c}{\textbf{K-means on Splines Coefficients}}\\Uniform & 0.50 &   2 &   0 &   0 &   0 & 500 & 0.46 & 0.19 \\ 
  Uniform & 2.00 &   2 &  31 &   0 &   1 & 468 & 0.47 & 0.18 \\ 
  Uniform & 0.50 &   3 & 189 &   0 &  10 & 301 & 0.48 & 0.17 \\ 
  Uniform & 2.00 &   3 & 500 &   0 &   0 &   0 & 0.48 & 0.15 \\ 
  Gaussian & 0.50 &   2 &   0 &   0 &   0 & 500 & 0.44 & 0.22 \\ 
  Gaussian & 2.00 &   2 &   0 &   0 &   0 & 500 & 0.45 & 0.20 \\ 
  Gaussian & 0.50 &   3 &   1 &   0 &   0 & 499 & 0.46 & 0.18 \\ 
  Gaussian & 2.00 &   3 & 490 &   0 &   0 &  10 & 0.47 & 0.16 \\ 
   \\ \multicolumn{9}{c}{\textbf{Independent Mixture}}\\Uniform & 0.50 &   2 &   0 &   0 &   2 & 498 & 0.44 & 0.22 \\ 
  Uniform & 2.00 &   2 &   0 &   0 &   2 & 498 & 0.40 & 0.28 \\ 
  Uniform & 0.50 &   3 &   0 &   0 &   1 & 499 & 0.49 & 0.12 \\ 
  Uniform & 2.00 &   3 &   0 &   0 &   2 & 498 & 0.49 & 0.12 \\ 
  Gaussian & 0.50 &   2 &   0 &   0 &  10 & 490 & 0.41 & 0.25 \\ 
  Gaussian & 2.00 &   2 &   0 &   0 &  24 & 476 & 0.39 & 0.29 \\ 
  Gaussian & 0.50 &   3 &   0 &   0 &   1 & 499 & 0.45 & 0.17 \\ 
  Gaussian & 2.00 &   3 &   0 &   0 &   5 & 495 & 0.45 & 0.17 \\ 
   \\ \multicolumn{9}{c}{\textbf{Exponential Mixture}}\\Uniform & 0.50 &   2 &   0 &  65 & 171 & 264 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   2 &   0 &   0 &  23 & 477 & 0.40 & 0.40 \\ 
  Uniform & 0.50 &   3 &   0 & 113 & 157 & 230 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   3 &   0 &   0 &  37 & 463 & 0.41 & 0.35 \\ 
  Gaussian & 0.50 &   2 &   0 & 114 & 136 & 250 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   2 &   0 &   0 &  34 & 466 & 0.41 & 0.38 \\ 
  Gaussian & 0.50 &   3 &   0 & 141 & 121 & 238 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   3 &   0 &   1 &  40 & 459 & 0.41 & 0.35 \\ 
   \thickhline\end{tabular}
\caption{The number of times each value of $K$ was chosen and the average misclassification rate (MR) and average Adjusted Rand Index (ARI) when $K=3$ for 500 replications of standard clustering methods applied to data generated under different conditions for the $F_{\lambda}$ and the standard deviation of $\epsilon$ and $\lambda$.}
\label{tab:freq1}

\end{table}

% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Feb 27 13:38:52 2013
\begin{table}[ht]
\centering
\begin{tabular}{ccc|cccccc}
  \thickhline $F_{\lambda}$&$\sigma_{\epsilon}$&$\sigma_{\lambda}$&$K=2$&$K=3$&$K=4$&$K=5$&MR&ARI\\  \hline\multicolumn{9}{c}{\textbf{K-means on Difference Quotients}}\\ Uniform & 0.50 &   2 &   0 & 500 &   0 &   0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   2 & 480 &  20 &   0 &   0 & 0.38 & 0.31 \\ 
  Uniform & 0.50 &   3 &   0 & 500 &   0 &   0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   3 & 480 &  20 &   0 &   0 & 0.38 & 0.31 \\ 
  Gaussian & 0.50 &   2 &   0 & 500 &   0 &   0 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   2 & 483 &  17 &   0 &   0 & 0.38 & 0.31 \\ 
  Gaussian & 0.50 &   3 &   0 & 500 &   0 &   0 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   3 & 483 &  17 &   0 &   0 & 0.38 & 0.31 \\ 
   \\ \multicolumn{9}{c}{\textbf{Correlation-based PAM}}\\Uniform & 0.50 &   2 & 380 &   0 &   0 & 120 & 0.24 & 0.49 \\ 
  Uniform & 2.00 &   2 & 500 &   0 &   0 &   0 & 0.27 & 0.46 \\ 
  Uniform & 0.50 &   3 & 380 &   0 &   0 & 120 & 0.24 & 0.49 \\ 
  Uniform & 2.00 &   3 & 500 &   0 &   0 &   0 & 0.27 & 0.46 \\ 
  Gaussian & 0.50 &   2 & 403 &   0 &   0 &  97 & 0.25 & 0.49 \\ 
  Gaussian & 2.00 &   2 & 500 &   0 &   0 &   0 & 0.27 & 0.45 \\ 
  Gaussian & 0.50 &   3 & 403 &   0 &   0 &  97 & 0.25 & 0.49 \\ 
  Gaussian & 2.00 &   3 & 500 &   0 &   0 &   0 & 0.27 & 0.45 \\ 
   \thickhline\end{tabular}
\caption{The number of times each value of $K$ was chosen and the average misclassification rate (MR) and average Adjusted Rand Index (ARI) when $K=3$ for 500 replications of clustering methods intended to group by shape applied to data generated under different conditions for the $F_{\lambda}$ and the standard deviation of $\epsilon$ and $\lambda$.}
\label{tab:freq2}

\end{table}

% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Feb 27 13:39:11 2013
\begin{table}[ht]
\centering
\begin{tabular}{ccc|cccccc}
  \thickhline$F_{\lambda}$&$\sigma_{\epsilon}$&$\sigma_{\lambda}$&$K=2$&$K=3$&$K=4$&$K=5$&MR&ARI\\ \hline\multicolumn{9}{c}{\textbf{K-means on Derivative Splines Coefficients}}\\ Uniform & 0.50 &   2 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   2 & 15 & 189 & 164 & 132 & 0.39 & 0.27 \\ 
  Uniform & 0.50 &   3 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   3 & 15 & 189 & 164 & 132 & 0.39 & 0.27 \\ 
  Gaussian & 0.50 &   2 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   2 & 18 & 191 & 171 & 120 & 0.39 & 0.27 \\ 
  Gaussian & 0.50 &   3 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   3 & 18 & 191 & 171 & 120 & 0.39 & 0.27 \\ 
   \\ \multicolumn{9}{c}{\textbf{Multilayer Mixture}}\\Uniform & 0.50 &   2 &  &  &  &  & 0.07 & 0.82 \\ 
  Uniform & 2.00 &   2 &  &  &  &  & 0.29 & 0.35 \\ 
  Uniform & 0.50 &   3 &  &  &  &  & 0.33 & 0.32 \\ 
  Uniform & 2.00 &   3 &  &  &  &  & 0.46 & 0.12 \\ 
  Gaussian & 0.50 &   2 &  &  &  &  & 0.09 & 0.75 \\ 
  Gaussian & 2.00 &   2 &  &  &  &  & 0.30 & 0.34 \\ 
  Gaussian & 0.50 &   3 &  &  &  &  & 0.31 & 0.34 \\ 
  Gaussian & 2.00 &   3 &  &  &  &  & 0.41 & 0.17 \\ 
   \\ \multicolumn{9}{c}{\textbf{Vertically Shifted Independent Mixture}}\\Uniform & 0.50 &   2 & 0 & 499 & 1 & 0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   2 & 0 & 499 & 1 & 0 & 0.05 & 0.87 \\ 
  Uniform & 0.50 &   3 & 0 & 499 & 1 & 0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   3 & 0 & 499 & 1 & 0 & 0.05 & 0.87 \\ 
  Gaussian & 0.50 &   2 & 0 & 499 & 0 & 1 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   2 & 0 & 498 & 2 & 0 & 0.05 & 0.87 \\ 
  Gaussian & 0.50 &   3 & 0 & 499 & 0 & 1 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   3 & 0 & 498 & 2 & 0 & 0.05 & 0.87 \\ 
   \\ \multicolumn{9}{c}{\textbf{Vertically Shifted Exponential Mixture}}\\Uniform & 0.50 &   2 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   2 & 0 & 500 & 0 & 0 & 0.05 & 0.87 \\ 
  Uniform & 0.50 &   3 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Uniform & 2.00 &   3 & 0 & 500 & 0 & 0 & 0.05 & 0.87 \\ 
  Gaussian & 0.50 &   2 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   2 & 0 & 499 & 1 & 0 & 0.05 & 0.87 \\ 
  Gaussian & 0.50 &   3 & 0 & 500 & 0 & 0 & 0.00 & 1.00 \\ 
  Gaussian & 2.00 &   3 & 0 & 499 & 1 & 0 & 0.05 & 0.87 \\ 
   \thickhline\end{tabular}
\caption{The number of times each value of $K$ was chosen and the average misclassification rate (MR) and average Adjusted Rand Index (ARI) when $K=3$ for 500 replications of proposed clustering methods applied to data generated under different conditions for the $F_{\lambda}$ and the standard deviation of $\epsilon$ and $\lambda$.}
\label{tab:freq3}

\end{table}
To compare the group representatives from all methods with the true mean shapes used to generate the data, we transformed the estimates when $K=3$ to coefficients of the derivative spline as previously described. There are two coefficients per group and I denote them as $\alpha_{11},\alpha_{21}$ for group 1, $\alpha_{12},\alpha_{22}$ for group 2, and $\alpha_{13},\alpha_{23}$ for group 3. Tables \ref{tab:mse1}, \ref{tab:mse2}, and \ref{tab:mse3} present the mean squared errors for the derivative spline coefficients under the simulation conditions for all clustering methods. When the clustering algorithm finds the shape groups correctly, the MSE is very close to zero for all coefficients. This occurs for the vertically shifted mixture under all conditions and both derivative-based methods and the exponential mixture model when the measurement error is small. Additionally, the coefficients in the multilayer mixture model have a small MSE only when the level groups are well separated.

 In general, the standard methods find one group with a horizontal shape on average indicated by values close to zero for $\alpha_{21}$ and $\alpha_{22}$. Large MSE values (close to or greater than 1) indicate clustering instability in that the method finds fundamentally different groups for every randomly generated data set resulting in great variation in the estimates (Table \ref{tab:mse2} and \ref{tab:mse3}). This occurs for a variety of reasons. In the case of the derivative spline coefficients, overfitting the five observations results in high variability. For the correlation-based method, the inconsistent dissimilarity between noisy horizontal vectors causes extra variability.

% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Feb 27 13:39:30 2013
\begin{table}[ht]
\centering
\begin{tabular}{ccc|cccccc}
  \thickhline $F_{\lambda}$&$\sigma_{\epsilon}$&$\sigma_{\lambda}$&$\alpha_{11}$&$\alpha_{21}$&$\alpha_{12}$&$\alpha_{22}$&$\alpha_{13}$&$\alpha_{23}$\\ \hline\multicolumn{9}{c}{\textbf{K-means on Original Data}}\\ Uniform & 0.50 &   2 & 0.71 & 0.71 & 0.09 & 0.09 & 0.72 & 0.72 \\ 
  Uniform & 2.00 &   2 & 0.78 & 0.77 & 0.06 & 0.06 & 0.80 & 0.80 \\ 
  Uniform & 0.50 &   3 & 0.87 & 0.88 & 0.01 & 0.01 & 0.88 & 0.88 \\ 
  Uniform & 2.00 &   3 & 0.88 & 0.89 & 0.01 & 0.01 & 0.89 & 0.88 \\ 
  Gaussian & 0.50 &   2 & 0.92 & 0.92 & 0.01 & 0.01 & 0.91 & 0.92 \\ 
  Gaussian & 2.00 &   2 & 0.90 & 0.91 & 0.02 & 0.02 & 0.89 & 0.90 \\ 
  Gaussian & 0.50 &   3 & 0.89 & 0.89 & 0.01 & 0.01 & 0.89 & 0.90 \\ 
  Gaussian & 2.00 &   3 & 0.89 & 0.89 & 0.01 & 0.01 & 0.88 & 0.91 \\ 
   \\ \multicolumn{9}{c}{\textbf{K-means on Splines Coefficients}}\\Uniform & 0.50 &   2 & 0.50 & 0.50 & 0.21 & 0.22 & 0.53 & 0.53 \\ 
  Uniform & 2.00 &   2 & 0.49 & 0.55 & 0.23 & 0.29 & 0.53 & 0.58 \\ 
  Uniform & 0.50 &   3 & 0.52 & 0.53 & 0.23 & 0.24 & 0.59 & 0.59 \\ 
  Uniform & 2.00 &   3 & 0.56 & 0.63 & 0.26 & 0.29 & 0.60 & 0.69 \\ 
  Gaussian & 0.50 &   2 & 0.49 & 0.49 & 0.19 & 0.20 & 0.56 & 0.56 \\ 
  Gaussian & 2.00 &   2 & 0.47 & 0.55 & 0.21 & 0.29 & 0.54 & 0.62 \\ 
  Gaussian & 0.50 &   3 & 0.68 & 0.69 & 0.14 & 0.14 & 0.69 & 0.70 \\ 
  Gaussian & 2.00 &   3 & 0.64 & 0.76 & 0.16 & 0.20 & 0.64 & 0.76 \\ 
   \\ \multicolumn{9}{c}{\textbf{Independent Mixture}}\\Uniform & 0.50 &   2 & 0.87 & 0.87 & 0.01 & 0.01 & 0.87 & 0.87 \\ 
  Uniform & 2.00 &   2 & 0.85 & 0.86 & 0.02 & 0.02 & 0.87 & 0.86 \\ 
  Uniform & 0.50 &   3 & 0.85 & 0.85 & 0.02 & 0.02 & 0.86 & 0.86 \\ 
  Uniform & 2.00 &   3 & 0.84 & 0.83 & 0.02 & 0.02 & 0.86 & 0.85 \\ 
  Gaussian & 0.50 &   2 & 0.88 & 0.88 & 0.01 & 0.01 & 0.88 & 0.88 \\ 
  Gaussian & 2.00 &   2 & 0.86 & 0.86 & 0.02 & 0.02 & 0.86 & 0.86 \\ 
  Gaussian & 0.50 &   3 & 0.91 & 0.91 & 0.01 & 0.01 & 0.90 & 0.90 \\ 
  Gaussian & 2.00 &   3 & 0.90 & 0.88 & 0.01 & 0.01 & 0.89 & 0.90 \\ 
   \\ \multicolumn{9}{c}{\textbf{Exponential Mixture}}\\Uniform & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   2 & 0.28 & 0.30 & 0.38 & 0.37 & 0.26 & 0.27 \\ 
  Uniform & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   3 & 0.12 & 0.13 & 0.39 & 0.37 & 0.10 & 0.09 \\ 
  Gaussian & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   2 & 0.26 & 0.25 & 0.41 & 0.40 & 0.24 & 0.24 \\ 
  Gaussian & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   3 & 0.12 & 0.12 & 0.39 & 0.37 & 0.11 & 0.11 \\ 
   \thickhline\end{tabular}
\caption{Mean squared error for derivative spline coefficients when $K=3$ for 500 replications of standard clustering methods applied to data generated under different conditions for the distribution of $\lambda$ and the standard deviation of $\epsilon$ and $\lambda$.}
\label{tab:mse1}

\end{table}
% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Feb 27 13:39:46 2013
\begin{table}[ht]
\centering
\begin{tabular}{ccc|cccccc}
  \thickhline $F_{\lambda}$&$\sigma_{\epsilon}$&$\sigma_{\lambda}$&$\alpha_{11}$&$\alpha_{21}$&$\alpha_{12}$&$\alpha_{22}$&$\alpha_{13}$&$\alpha_{23}$\\  \hline\multicolumn{9}{c}{\textbf{K-means on Difference Quotients}}\\ Uniform & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   2 & 0.10 & 0.09 & 0.20 & 0.17 & 0.10 & 0.10 \\ 
  Uniform & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   3 & 0.10 & 0.09 & 0.20 & 0.17 & 0.10 & 0.10 \\ 
  Gaussian & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   2 & 0.10 & 0.09 & 0.20 & 0.18 & 0.10 & 0.10 \\ 
  Gaussian & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   3 & 0.10 & 0.09 & 0.20 & 0.18 & 0.10 & 0.10 \\ 
   \\ \multicolumn{9}{c}{\textbf{PAM with Correlation Measure}}\\Uniform & 0.50 &   2 & 0.91 & 0.91 & 0.46 & 0.46 & 1.88 & 1.85 \\ 
  Uniform & 2.00 &   2 & 1.22 & 1.17 & 0.64 & 0.66 & 2.14 & 2.02 \\ 
  Uniform & 0.50 &   3 & 0.91 & 0.91 & 0.46 & 0.46 & 1.88 & 1.85 \\ 
  Uniform & 2.00 &   3 & 1.22 & 1.17 & 0.64 & 0.66 & 2.14 & 2.02 \\ 
  Gaussian & 0.50 &   2 & 0.83 & 0.87 & 0.47 & 0.46 & 1.77 & 1.77 \\ 
  Gaussian & 2.00 &   2 & 1.20 & 1.31 & 0.67 & 0.64 & 1.95 & 2.12 \\ 
  Gaussian & 0.50 &   3 & 0.83 & 0.87 & 0.47 & 0.46 & 1.77 & 1.77 \\ 
  Gaussian & 2.00 &   3 & 1.20 & 1.31 & 0.67 & 0.64 & 1.95 & 2.12 \\ 
   \thickhline\end{tabular}
\caption{Mean squared error for derivative spline coefficients when $K=3$ for 500 replications of clustering methods intended to group by shape applied to data generated under different conditions for the distribution of $\lambda$ and the standard deviation of $\epsilon$ and $\lambda$.}
\label{tab:mse2}

\end{table}
% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Feb 27 13:40:02 2013
\begin{table}[ht]
\centering
\begin{tabular}{ccc|cccccc}
  \thickhline $F_{\lambda}$&$\sigma_{\epsilon}$&$\sigma_{\lambda}$&$\alpha_{11}$&$\alpha_{21}$&$\alpha_{12}$&$\alpha_{22}$&$\alpha_{13}$&$\alpha_{23}$\\  \hline\multicolumn{9}{c}{\textbf{K-means on Derivative Splines Coefficients}}\\ Uniform & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   2 & 0.50 & 0.49 & 1.21 & 1.29 & 0.41 & 0.51 \\ 
  Uniform & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   3 & 0.50 & 0.49 & 1.21 & 1.29 & 0.41 & 0.51 \\ 
  Gaussian & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   2 & 0.46 & 0.45 & 1.26 & 1.27 & 0.48 & 0.52 \\ 
  Gaussian & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   3 & 0.46 & 0.45 & 1.26 & 1.27 & 0.48 & 0.52 \\ 
   \\ \multicolumn{9}{c}{\textbf{Multilayer Mixture}}\\Uniform & 0.50 &   2 & 0.02 & 0.02 & 0.01 & 0.01 & 0.03 & 0.03 \\ 
  Uniform & 2.00 &   2 & 0.12 & 0.12 & 0.01 & 0.01 & 0.12 & 0.11 \\ 
  Uniform & 0.50 &   3 & 0.33 & 0.33 & 0.03 & 0.03 & 0.33 & 0.32 \\ 
  Uniform & 2.00 &   3 & 0.40 & 0.40 & 0.02 & 0.02 & 0.42 & 0.40 \\ 
  Gaussian & 0.50 &   2 & 0.01 & 0.01 & 0.00 & 0.00 & 0.01 & 0.01 \\ 
  Gaussian & 2.00 &   2 & 0.13 & 0.13 & 0.01 & 0.01 & 0.13 & 0.13 \\ 
  Gaussian & 0.50 &   3 & 0.31 & 0.32 & 0.02 & 0.02 & 0.29 & 0.29 \\ 
  Gaussian & 2.00 &   3 & 0.35 & 0.34 & 0.01 & 0.01 & 0.34 & 0.35 \\ 
   \\ \multicolumn{9}{c}{\textbf{Vertically Shifted Independent Mixture}}\\Uniform & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   2 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
  Uniform & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   3 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
  Gaussian & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   2 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
  Gaussian & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   3 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
   \\ \multicolumn{9}{c}{\textbf{Vertically Shifted Exponential Mixture}}\\Uniform & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   2 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
  Uniform & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Uniform & 2.00 &   3 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
  Gaussian & 0.50 &   2 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   2 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
  Gaussian & 0.50 &   3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 
  Gaussian & 2.00 &   3 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 
   \thickhline\end{tabular}
\caption{Mean squared error for derivative spline coefficients when $K=3$ for 500 replications of proposed clustering methods applied to data generated under different conditions for the distribution of $\lambda$ and the standard deviation of $\epsilon$ and $\lambda$.}
\label{tab:mse3}

\end{table}
I have determined that some of the methods discover the shape groups under conditions of small measurement error by comparing the group labels to the true generating labels comparing the estimates of the mean shapes with known generating values. Lastly, I want to know how the choice of method impacts the estimation of the relationship between baseline factors and group membership. For simplification, we present density plots of the estimated coefficients for group 1 and 2 for the factors $w_{1}$ and $w_{2}$ for only three methods---K-means on the original data, K-means on the difference quotients, and vertically shifted exponential mixture models---after setting $K=3$. As a reminder, the first baseline factor $w_{1}$ was used to generate shape groups and the second $w_{2}$ was used to generate level groups. 

It is clear that the K-means on the original data picks up the relationship between the group membership and the second baseline factor $w_{2}$ (Figure \ref{fig:gamma1}). The density plots with dashed lines for the coefficients of $w_{2}$ are concentrated at non-zero values while the solid line densities for the coefficients of $w_{1}$ revolve around zero. For the two uniform conditions with $\sigma_{\lambda}=2$, the density plots for the coefficients of $w_{1}$ are bimodal, which indicates some variability in the groups chosen. Overall, K-means picks up the level groups since the coefficients for $w_{2}$ are significantly non-zero.

On the other hand, the K-means on the difference quotients finds a relationship with $w_{1}$ (Figure \ref{fig:gamma2}). The average coefficient values for $w_{1}$ are non-zero while the coefficients for $w_{2}$ revolve around zero. However, when $\sigma_{\epsilon}=2$, one of the density plots is bimodal indicating variability which was indicated by a higher misclassification rate.

Lastly, with the vertically shifting exponential mixture model (Figure \ref{fig:gamma3}), there is a significant relationship between group membership and $w_{1}$, the factor that was used to generate shape groups. The second factor is estimated to have no relationship with the final grouping.

\begin{landscape}
\begin{figure}[ht]
\centering
\includegraphics[height=5.5in]{Chp5Gamma1}
\caption{Density plots of estimated multinomial coefficients from analysis after running K-means on original data under  simulation conditions. Coefficients for group 3 are fixed equal to 0. Solid: $w_{1}$, dashed: $w_{2}$. Thin: group 1, thick:  group 2.}
\label{fig:gamma1}
\end{figure}

\begin{figure}[ht]
\centering
\includegraphics[height=5.5in]{Chp5Gamma2}
\caption{Density plots of estimated multinomial coefficients from analysis after running K-means on difference quotients under simulation conditions. Coefficients for group 3 are fixed equal to 0. Solid: $w_{1}$, dashed: $w_{2}$. Thin: group 1, thick:  group 2.}
\label{fig:gamma2}
\end{figure}
\begin{figure}[ht]
\centering
\includegraphics[height=5.5in]{Chp5Gamma3}
\caption{Density plots of simultaneously estimated multinomial coefficients from the vertically shifted exponential mixture model under simulation conditions. Coefficients for group 3 are fixed equal to 0. Solid: $w_{1}$, dashed: $w_{2}$. Thin: group 1, thick:  group 2.}
\label{fig:gamma3}
\end{figure}
\end{landscape}

\section{Discussion}
This study provides evidence against standard clustering methods and insight into the strengths and weaknesses of the proposed methods to cluster by shape. K-means on original data finds groups based mainly on level. K-means on spline coefficients generalizes the original K-means to work with irregularly sampled data, but does not improve performance when the number of observations is small. The finite mixture model with independence also fails to detect shapes when shape and level are not perfectly correlated. Surprisingly, the mixture model with exponential correlation was able to detect the three shape groups when measurement error was small, but it fails to select the correct number of groups.

Of the two methods previously thought to detect shape, PAM with the correlation measure failed on all accounts in this simulation. Not only did it choose the incorrect number of shape groups, it could not find the shape groups when forced to have the correct number of shape groups. Quotient differences, on the other hand, worked well when measurement error was small. It was able to choose the correct number of groups and find the shape groups.

The K-means on the derivative spline coefficients is an attempt to use the derivative for clustering while accommodating irregularly sampled data. For data with small measurement error, this method performed well. However, under larger measurement error, it performed worse than difference quotients. This primarily was due to the small number of observation times over which the data are smoothed for each individual. I imagine the performance would improve with more observations over time.  The multilayer mixture model worked well when the data fit into the shape group, level component framework such that the components are well separated. Unfortunately, this is not a realistic for most data sets. Vertically shifting the data prior to fitting a mixture model outperformed the rest of the methods. Measurement error and overlap had little impact on the results. It consistency found the three shape groups even when the covariance was not modeled perfectly.

These results are limited as I focused on only one generating framework. The shapes are straight lines with very few observations. While this is restrictive, it does not drastically deviate from real data sets with long range trends. An interesting extension of this simulation would be to change the slopes of the lines to further explore how the separation impacts the proposed models. I kept the data generating correlation simple by including a random intercept on top of independent errors. One could try a more complex structure to detect the impact of covariance misspecification like in Chapter \ref{chap:misspecify}. Despite these limitations, the conclusions should generalize to similar situations.