Since first introduced by (Lazarsfeld 1950a,b), LCA has been undergoing significant development (e.g., Goodman, 1974; Vermunt and Magidson, 2004). Typically, LCA assumes that the population that is being studied is a mixture of two or more mutually exclusive subsgroups (i.e., classes). The subgroup memberships of individuals are unobserved and lead to different distributions of the observed variables. Let x_i = (x_i1, …, x_iv, …, x_iV), where 1 ≤ v ≤ V, denote the observed responses of individual i (1 ≤ i ≤ N) to V polytomous items/variables, where N indicates the number of individuals in the sample. Invoking the conditional independence assumption, we have that for an individual i that belongs to subgroup k (1 ≤ k ≤ K), the probability that he/she/they have the response pattern x_i is

where P(x_ij|c_i = k) is the probability that a subgroup k individual score x_ij for item j. Let λ = (λ₁, …, λ_K) be a vector that consists of the probabilities that an individual belongs to each of the K subgroups. The multivariate density of the individual i can then be specified as a mixture distribution,

Then the overall likelihood L is the product of the contribution of each individual,

Estimates of model parameters are generally obtained through maximum likelihood (ML)-based estimation procedures, such as the expectation–maximization (EM) algorithm.

The k-means clustering is probably the most widely-applied unsupervised machine learning technique, and k-medians clustering can be viewed as one of its variants. Unlike LCA, k-means and k-medians do not assume any statistical models but assign individuals to the pre-specified number of subgroups through minimizing the sum of within-cluster variation. Specifically, with k-means, the variation within each cluster is computed as the sum of the squared Euclidean distance between each individual and its cluster center,

where x̄i represents the means of all individuals within subgroup k. On the other hand, the k-medians clustering uses the medians as cluster centers and is therefore more robust to outliers. A smaller W_k indicates that the individuals in the subgroup k are more homogeneous. Note that while the Euclidean distance is often considered the default of k-means and k-medians, other distance measures, such as the Manhattan distance, can also be used.

This minimization problem of k-means and k-medians is solved through applying the iterative algorithm shown in Algorithm 1. The algorithm starts with randomly assigning individuals to subgroups and updates the subgroup memberships based on individuals' distance to the cluster centers. The algorithm stops when the membership assignments stop changing. This optimization algorithm may stop at the local optimum rather than the global optimum. Therefore, in practice, one needs to run the algorithm multiple times with different initial configurations and then select the clustering solution with the smallest sum of within-cluster variation.

In addition to AIC and BIC for LCA, and MRPC, CH pseudo-F statistic and Silhouette method for k-means and k-medians, there are other methods for determining the optimal number of subgroups while tackling the clustering tasks.

Sample-size adjusted BIC. The sample-size adjusted BIC (SABIC; Sclove, 1987) is another popular information criterion and is used for choosing the number of subgroups when LCA is applied. SABIC is similar to BIC; however, the penalty on N is reduced,

Gap statistic. For k-means and k-medians, the Gap statistic (Tibshirani et al., 2002) is another option for determining the optimal number of subgroups. The Gap statistic compares the curve of log(W[K]) to the reference curve obtained from data uniformly distributed over a rectangle containing the data. Then the optimal number of subgroups is chosen as the K associated with the largest difference between the two curves.

To simulate data, we followed the process used in previous simulation studies (e.g., Brusco, 2004; Brusco et al., 2017; Dimitriadou et al., 2002), and manipulated six factors, including (1) the number of individuals N, (2) the number of subgroups K, (3) the number of true variables V, (4) the number of masking variables P, (5) the error level in the subgroup structure ϵ, and (6) the densities of subgroups.

Motivated by previous studies, we generated item responses in the following three steps:

Sample R code for data generation were included in online Supplementary material.

Number of individuals N. We considered three levels of the number of individuals, N = 250, 500, and 1,000. These numbers are generally viewed as moderately large/large sample sizes for LCA, k-means and k-medians clustering.

Number of subgroups K. While generating values of true variables, we considered three numbers of subgroups, K = 2, 3, and 4. The cluster numbers were chosen based on previous simulation studies (e.g., Brusco, 2004; Brusco et al., 2017; Dimitriadou et al., 2002).

Number of true variables V. In this simulation study, we considered three numbers of true variables (i.e., variables that define the subgroup structure), V = 10, 20, and 30. These numbers covered a wide range the lengths of psychological assessments and educational survey questionnaires.

Number of masking variables P. In practice, it is possible that not all variables in a measure define the subgroup structure. Therefore, in addition to conditions where there exist no masking variable (i.e., P = 0), we considered two levels of the number of masking variables, P = 5 and 10.

Level of error ϵ. Following previous simulation studies, we considered three levels of errors in Step 2 of the data generation process, ϵ = 5%, 10%, and 20%.

Subgroup probabilities. We considered three types of cluster probabilities, one of which represents cases with equal probability and two for unequal probability conditions. Specifically, in the equal probability conditions, all subgroups are of the same size so that the probability that an individual belong to one of the K subgroups in 1/K. In the first type of unequal probability conditions, the first subgroup had a probability of 0.6, and the other K−1 subgroups have the same probability of (1 − 0.4)/(K−1). In the second type of unequal probability conditions, the first subgroup had a probability of 0.1, and other clusters have the probability of (1 − 0.1)/(K−1).

In sum, these six manipulated factors led to a total of 3 × 3 × 3 × 3 × 3 × 3 = 729 simulation conditions. For each of the conditions, we simulated and analyzed 50 data sets in R (R Core Team, 2023).

We analyzed each simulated data set using the three approaches of LCA, k-means and k-medians clustering, extracting 1 to 6 subgroups, and chose the optimal number of subgroups based on different methods. Specifically, we performed LCA using the PoLCA package (Drew and Jeffrey, 2011) and extracted the AIC and BIC indices. For k-means, we used the R package cluster (Maechler et al., 2013) for clustering and computed the MRPC, CH pseudo-F statistic, and Silhouette statistic based on the outputs. For k-medians clustering, we used the Kmedians package (Godichon-Baggioni, 2022) for subgroup extraction and obtained the CH pseudo-F statistic and Silhouette statistic using the cluster package. For all three approaches and numbers of subgroups, the numbers of random starting values/configurations were fixed at 100, the max numbers of iterations were set to 1,000, unless otherwise specified defaults were used.

In this study, we focused on the key indicator of subgroup partition recovery to evaluate the three clustering approaches and their associated methods for choosing the number of subgroups. Following previous simulation studies, for each simulation condition, we computed the Adjusted Rand Index (ARI; Hubert, 1974; Steinley, 2004) between the true subgroup memberships and predicted partitions for each combination of clustering approaches and methods for choosing the subgroup number using the Kmedians package (Azzalini and Menardi, 2014). The ARI index quantifies the agreement between two partitions and ranges from 0 to 1, with a larger value indicating a higher level of agreement. We then performed the analysis of variance (ANOVA) to investigate how the ARI was impacted by the clustering approaches and the manipulated factors. We considered main effects of the six manipulated factors and all two-way interactions. To better inform about the performances of these approaches, we also plotted the determined number of subgroups for each combination of clustering approaches and methods of choosing subgroup number.

As shown in Table 4, the ARIs of the seven combinations of clustering approaches and methods for selecting the optimal subgroup number (i.e., Methods) were significantly different (F = 375.6, df = 6, p < 0.01, partial η² = 0.31). The post hoc Tukey's HSD test indicated that LCA outperformed the k-means and k-medians clustering, with BIC led to a greater ARI than AIC. Additionally, k-means coupled with the MRPC performed significantly better than other four methods, including k-means and k-medians coupled with the CH Pseudo-F statistic and the Silhouette method.

Besides, five out of the six manipulated simulation factors showed significant main effects on the ARI, including the number of clusters K (F = 1352.7, df = 2, p < 0.01, partial η² = 0.35), the number of true variables V (F = 299.2, df = 2, p < 0.01, partial η² = 0.11), the number of masking variables P (F = 99.8, df = 2, p < 0.01, partial η² = 0.04), the error level E (F = 568.2, df = 2, p < 0.01, partial η² = 0.18), and the cluster density D (F = 65.9, df = 2, p < 0.01, partial η² = 0.03). Lastly, ten two-way interactions were significant; however, the effect sizes were all smaller than 0.04.

Figure 1 presents the aggregated ARI values by three out of the six manipulated factors with effect sizes that were greater than 0.1, including the number of clusters K, the number of true variables V, and the error level E. A general trend was that the ARI value decreases as K increases, V decreases, and E decreases. For instance, as shown in Figure 1b, when V increases, the ARI value of k-means coupled with MRPC improved from 0.7 to 0.93. Another example is that, as shown in Figure 1c, when E increases, the ARI of k-medians coupled with CH Pseudo-F decreased from 0.81 to 0.66. It is also worth noting that the LCA approach had strong subgroup recovery performance in V = 10 and E = 0.2 conditions, among which the lowest ARI value was 0.81.

Figure 2 shows boxplots of the determined numbers of subgroups by combinations of clustering approaches and associated methods, aggregated over all simulation factors other than K. As shown in Figure 2a, when K = 2, LCA coupled with BIC and k-means and k-medians tend to choose 2 as the optimal number of sub-groups; while LCA coupled with AIC overestimated K. When the true number of subgroups were 3 or 4, LCA coupled with BIC was still able to select the correct subgroup number; however, k-means and k-medians coupled with MRPC and CH Pseudo-F tend to underestimated K and selected 2 as the number of subgroups, while LCA coupled with AIC favored more subgroups.

We examined a parent-rated measurement instrument named the Vineland Adaptive Behavior Scales Third Edition (Vineland-3; Sparrow et al., 2016) in the SPARK data base. The Vineland-3 is used as a measure of adaptive skills and is often a part of diagnostic evaluations for intellectual and developmental disabilities. Participant data were drawn from the SPARK data repository that includes phenotype and genotype data for individuals with an autism spectrum diagnosis. Thus, participants included in the empirical illustration all had a diagnosis of an autism spectrum disorder (ASD) and were all under the age of 18.

In this analysis, we focused on one of the 13 subdomains in Vineland-3, the Domestic subdomain, which includes 30 items. Items in this sub-domain measure a child's ability to perform household tasks (e.g., [The individual] is careful when using sharp objects, for example, scissor, knives). Parents rated the children's behaviors on a 3-point Likert-type scale (0 = Never, 1 = Sometimes, and 2 = Usually or often). A sample of N = 500 individuals were randomly selected from complete observations for the analysis.

We performed LCA, k-means and k-medians clustering on the sample data in R (R Core Team, 2023). Specifically, for each of the three approaches, we ran a series of analyses with different numbers of subgroups, ranging from 1 to 10. Based on the outputs, we chose the optimal number of subgroups using various methods, including AIC and BIC for LCA, and MPRC, CH pseudo-F and Silhouette statistics for k-means and k-medians. Then we computed the agreements between partitions with the subgroups numbers determined by different methods.

Table 5 summarizes various indices for the three approaches with different numbers of subgroups. The first observation is that the number of subgroups selected based on different approaches varied. Based on LCA coupled with AIC, the number of subgroups was 9, and the selected number of subgroups was 3 when BIC was applied. On the other hand, the numbers of subgroups determined by k-means and k-medians coupled with MRPC, CH pseudo-F and Silhouette statistics were all 2.

Table 6 summarizes the ARIs between subgroup partitions determined based on different methods. The agreements between LCA coupled with AIC and other methods were low, with the ARI value ranging from 0.15 to 0.24. The agreements of LCA coupled with BIC and k-means and k-medians were moderate. The ARI between k-means and k-medians was 0.98, suggesting a high level of agreement between the two partitions.

As with other simulations, one of the limitations of the present study concerns the manipulated factors and the chosen levels of these manipulated factors. For instance, the sample sizes N we considered in this simulation (250, 500, and 1,000) are generally viewed as moderately large and large sample sizes in psychological research. Therefore, our results may not be generalizable to cases where the sample sizes are relatively small. We also designed only one subgroup structure for each combination of the number of true variables V and the number of subgroups K. In reality, the distances between subgroups tend to be different. In addition, we simulated responses to items with three categories, while many psychological assessments and education surveys consists of items with four or more response categories. Therefore, for future research, we suggest more investigations on the impacts of true subgroup structures and measurement scale properties on the performance of clustering approaches.

Second, while maintaining consistency with prior studies, the process we followed to generate data is prone to several issues. For instance, we introduced randomness in Steps 2 and 3 of the process, but duplicates may exist and result in numerical instability of the optimization algorithms. We advocate for more discussion in the psychological and educational measurement community on how to generate data that can better reflect the reality. In addition, when different approaches are to be compared via a simulation, the data generation process should be carefully designed so that the simulated data do not favor any of the approaches.

Last but not least, some other popular clustering approaches were not considered due to the scope of this study. These approaches include but are not limited to HCA, spectral clustering (e.g., Ng et al., 2002; Rohe et al., 2011), and autoencoders for clustering (e.g., Klingler et al., 2017; Zhang et al., 2022). In addition, to better align with existing studies (e.g., Brusco et al., 2017), we did not include all possible methods for determining the optimal number of subgroups for the three clustering approaches in our simulation and empirical data analysis, such as SABIC and the Gap statistic.

To enhance the understanding of LCA, k-means and k-medians when applied to psychological and educational research, we suggest expanding the scope of the simulation to account for various item types, true subgroup structures, and sample sizes. Future studies can also include factors related to psychometric properties of measurement scales, such as their reliabilities and item properties, in the simulation.

Secondly, to provide practitioners with a more comprehensive picture, we recommend comparing the three approaches considered in this study with other approaches for clustering tasks and other methods for choosing the subgroup numbers. Moreover, characteristics of the approaches, such as available software programs and the computation time, should be discussed in addition to their performance in recovering the subgroup partition. We also suggest that practitioners compare the solutions of different approaches (e.g., numbers of subgroups and assignment of subgroup memberships) when completing clustering tasks and choose the more interpretable results.