Cluster Analysis

Loading the data

Paste columns directly from the sheet (including the header row) – select them in Excel, copy, and paste here:

Variables and settings

When to use?

When the question is about people and not about variables: are there distinct types among the respondents – passive consumers vs. content creators, students who need support vs. those who do not – without such groups having been defined in advance? Wendler and Gröttrup (2021, p. 623) define cluster analysis as a collection of multivariate methods that "are used to identify groups of 'similar' objects", and list among the examples customer segmentation in marketing and identifying "subgroups of students with special needs" by background and performance (p. 624). The most important distinction, in their Figure 7.1: cluster analysis groups the rows (the participants), and factor analysis groups the columns (the items). If the question is which items measure the same thing – confirmatory factor analysis and the unsupervised learning page (PCA) and the exploratory factor analysis calculator; if the groups are known in advance and the question is whether they differ – analysis of variance. There is no predicted variable here and no "right answer": this is unsupervised learning, and "since the group membership is not known in advance, the algorithm cannot learn the characteristics of the groups" (Wendler & Gröttrup, 2021, p. 627).

The unsupervised learning page surveys all the clustering methods (k-means, hierarchical, Gaussian mixtures, DBSCAN, TwoStep) and the question of whether there are clusters at all. This page runs the most common of them in questionnaires: Ward followed by k-means, on quantitative variables. Categorical variables (e.g., sex, field of study) do not enter the distance computation here; they are used afterwards to validate and describe the clusters.

The logic: two steps

Step one – agglomerative hierarchical clustering. Start with each participant as a separate cluster, and at each step merge the two closest clusters, until all are in one cluster (Hastie et al., 2009, p. 520; Everitt et al., 2011, p. 73). Each stage is a partition, and the whole tree is a dendrogram: "the height of each node is proportional to the value of the intergroup dissimilarity between its two daughters", and a horizontal cut of the tree gives a partition into clusters (Hastie et al., 2009, p. 521). What does "closest" mean when it comes to groups? That is the linkage (Everitt et al., 2011, pp. 76-77, Figure 4.3): single (the closest pair), complete (the farthest pair), average (all pairs) – and Ward, in which the two clusters merged are those whose merger gives the smallest "increase in the total within-cluster sum of squares… over all variables" (Ward, 1963, as cited in Everitt et al., 2011, p. 78, Table 4.1; and likewise Wendler & Gröttrup, 2021, p. 648, Table 7.24). Ward is the linkage that suits the objective of k-means: both minimize the same thing – the sum of squared distances from the centers within the clusters. Single linkage tends to "chain" participants into long clusters, and complete linkage – to compact clusters whose members are sometimes closer to another cluster (Hastie et al., 2009, pp. 524-525).

Step two – k-means. Hierarchical clustering has a built-in drawback: "the assignment of an object to a certain cluster can never be changed" (Wendler & Gröttrup, 2021, p. 628, Table 7.1) – a participant merged early with a neighbor stays there even if in the end they are closer to the center of another cluster. k-means corrects this: given centers, assign each participant to the nearest center; given an assignment, compute each center as the mean; and repeat until no assignment changes (Hastie et al., 2009, p. 510, Algorithm 14.1). Each step reduces the sum of squared distances, so convergence is assured – "however, the result may represent a suboptimal local minimum", and it depends on the starting point (ibid.). Here the two steps meet: Géron (2022, Chapter 9, "Centroid initialization methods") – "if you happen to know approximately where the centroids should be (e.g., if you ran another clustering algorithm earlier), then you can set" them as the starting point. And likewise Everitt and colleagues (2011, p. 122): the initial partition "might be the result of a previous application of some other clustering method, perhaps one of the hierarchical techniques". Kassambara (2017, p. 163) has a name for the method – "hierarchical k-means", a hybrid method "for improving k-means results", born from the fact that "the final k-means clustering solution is very sensitive to the initial random selection of cluster centers", and its three steps are exactly these: hierarchical clustering and cutting the tree into k clusters, the center of each one, and k-means from these centers – where "k-means will improve the initial partitioning", so it "can be slightly different" from the final one. The page does exactly this: it cuts Ward's tree into k clusters, computes the centroid of each one, and runs k-means from them. The result is deterministic – same data, same result, in every run – and it does not have the random initialization that Wendler and Gröttrup (2021, p. 677) warn about in Modeler's K-Means, which starts from the first record so that "the result depends to some extent on the order of the objects". On the other hand, a single starting point is also a limitation: "the results from an optimization method can be dramatically affected by the choice of initial partition", and "it is advisable to run an optimization algorithm several times with different initial partitions"; and if different initial partitions give very different groupings, this is usually a sign that the number of clusters is wrong, or that there is no clear evidence of clusters (Marriott, 1982, as cited in Everitt et al., 2011, p. 122). The page reports two small checks in this direction – the comparison between the two steps and the single-move check (below); a run from many random initializations – in R, with nstart.

Wendler and Gröttrup (2021, p. 628, Table 7.1) summarize the advantages and disadvantages of the two families, and the two-step method takes from each what it is good at: from the hierarchical one – a tree that shows the structure at all levels and a "deterministic (not random)" result, and from k-means – the flexibility of "reassigning objects to other clusters". The comparison between the two steps is also a small check: the page reports how many participants changed cluster in the second step and the adjusted Rand index (ARI) between the two partitions – 1 is full agreement. This is a check of the algorithm on the same data, not of the stability of the solution in another sample – on that, see the Assumptions and warnings section.

Standardization and distance

The distance. Both methods work on the Euclidean distance between participants – the square root of the sum of squared differences on all variables, "the most commonly used distance measure" (Everitt et al., 2011, p. 49, Equation 3.4) – and on its square, which is "more sensitive to outliers" (Wendler & Gröttrup, 2021, pp. 630, 637). From this follows the important question on the page: what enters the distance, and with what weight. Hastie and colleagues (2009, p. 506) state that "specifying an appropriate dissimilarity measure is far more important in obtaining success with clustering than choice of clustering algorithm".

Standardization – yes, usually. In Euclidean distance, a variable with a large range dominates: hours per day (0-12) will "drown out" a scale score (1-5), and the clusters will in fact be clusters of hours. Therefore, as Everitt and colleagues (2011, p. 67) write, "it would not make sense to treat, say, weight measured in kilograms, height measured in metres, and anxiety rated on a four-point scale as equivalent", and "the solution most often suggested… is simply to standardize each variable to unit variance (variance 1) prior to any analysis"; Wendler and Gröttrup (2021, pp. 656-657) – standardizing the numeric variables "is very important for the clustering procedure… to make the values comparable, they have to be standardized"; and Géron (2022, Chapter 9, "Limits of k-means") – "it is important to scale the input features before you run k-means, or the clusters may be very stretched". The page's default is z-scores (mean 0 and standard deviation 1 for each variable), so that each variable has equal weight in the distance.

Standardization – not always. Hastie and colleagues (2009, p. 506) warn against this rule: giving all variables equal influence, "although this may seem reasonable, and is often recommended, it can be highly counterproductive". If one variable discriminates well between the groups and another is noise, standardization increases the weight of the noise, and "will tend to obscure the two well-separated groups to the point where a clustering algorithm cannot uncover them" – in their Figure 14.5 two well-separated groups disappear after standardization. Everitt and colleagues (2011, pp. 66-67) phrase the same thing: standardization is a special case of weighting, and "when variables are standardized prior to analysis, the investigator assumes that the importance of a variable decreases with increasing variability"; weights based on the standard deviations "appear ineffective"; dividing by the variable's range is preferable to z-scores in many applications (Milligan & Cooper, 1988, as cited there); and it is better to standardize by the within-group variance rather than by the overall variance – but the groups are exactly what is being sought. And when all variables are already on the same scale (e.g., four scale scores 1-5), standardization changes the weights according to the standard deviations in the sample – not always what is wanted. The page offers z-scores or the original units; standardization by range is not on the page. The decision depends on the variables and the question ("there can be no substitute for careful thought in the context of each problem", ibid.); on the page, check or uncheck "Standardize the variables", and the output states what was chosen. When standardization is off and the ranges of the variables differ greatly – the page warns.

Likert scale. The distance assumes that the difference between 1 and 2 equals the difference between 4 and 5. Wendler and Gröttrup (2021, p. 635, Table 7.8) note that on a Likert scale "the distances should be interpreted with caution, or not at all", and that the question is controversial in the literature. By the site's convention (the Variables and measurement scales page), a multi-item scale score – a mean of several items – is treated as an interval variable, and a single Likert item – is not. Hence the recommendation: cluster by scale scores (the Cronbach's alpha and scale score page), not by single items.

How many clusters?

This is the hard decision in the analysis, and there is no significance test for it. Hastie and colleagues (2009, p. 518) distinguish between two situations: in segmentation, "k is usually defined as part of the problem" (a company with k salespeople); and when the clustering is a description of "natural groupings", k is unknown and must be estimated from the data. The page runs the two-step solution for every k from 1 up to the maximum you chose, and prints three criteria side by side:

The elbow. The sum of squared distances from the centers within the clusters, W_k, always decreases as k grows – at the extreme, each participant is a cluster of their own and W = 0 (Hastie et al., 2009, p. 519). But if there really are k* groups, the decrease will be sharp up to k* and moderate after it, because beyond it every additional cluster only splits a natural group in two; look for a "kink" in the curve – "as with other aspects of clustering procedures, this approach is somewhat heuristic" (Hastie et al., 2009, p. 519); Everitt and colleagues (2011, p. 126) add that it "can be very subjective", and "large" is often "a function of the user's prior expectations". In the table it is also shown as explained variance, 1 − W_k/W₁ – the proportion of the total variance that the partition explains (the equivalent of R²; Wendler & Gröttrup, 2021, p. 649, Table 7.25).

The silhouette. For each participant: a – the average distance to the members of their own cluster, b – the average distance to the members of the next nearest cluster, and the silhouette s = (b − a)/max(a, b), between 1− and 1+. Rousseeuw (1987, p. 56), who defined the measure: s close to 1 – the participant is "well-clustered"; around 0 – an "intermediate case", where "it is not clear at all" to which of the two clusters it belongs; close to 1− – it can almost be concluded that it has been "misclassified". For a participant who is alone in a cluster the silhouette is 0 – a choice that is "of course arbitrary, but the value zero appears to be most neutral" (ibid.). The measure suits distances on a ratio scale (like Euclidean distance) and when "one is seeking compact and clearly separated clusters" (p. 55), and the average silhouette "can be used to select an 'appropriate' number of clusters" (p. 53). The average silhouette over all participants is a measure of the quality of the partition, and the k with the highest average is a candidate; for the cutoffs there are two sources that mostly agree: Kaufman and Rousseeuw see a "reasonable classification" in an average silhouette above 0.5 (Kaufman & Rousseeuw, 1990, as cited in Everitt et al., 2011, p. 129), and the SPSS Modeler documentation, as cited by Wendler and Gröttrup (2021, p. 649), grades up to .25 poor, between .25 and .50 fair, above .50 good (IBM, 2015, as cited there). The silhouette is not defined for k = 1, so it cannot say "there are no clusters here".

The gap statistic (Tibshirani et al., 2001, as cited in Hastie et al., 2009, p. 519) "compares the curve log W_k to the curve obtained from data uniformly distributed over a rectangle containing the data", that is, from data with no clusters. The gap is the difference between the two; "it is essentially an automatic way of locating the aforementioned 'kink'", and it "also works reasonably well when the data fall into a single cluster, and in that case will tend to estimate the optimal number of clusters to be one – a scenario where most other competing methods fail" (ibid.). The rule (equation 14.39 there): the smallest k at which the gap is not lower than the gap at k + 1 minus its standard error. The page generates the uniform data with R's generator (seed 2026 by default, 100 samples), so that the result is identical to that of clusGap in the cluster package with the same seed (R guide below). It is the only one of the three criteria that can answer "one": Everitt and colleagues (2011, p. 129) stress that it "allows the quality of a single-cluster solution to be assessed", that is, to ask whether there is any evidence of distinct clusters at all – which is why it is the default for choosing k in the output.

The general recommendation is "not to rely on a single rule for selecting the number of groups, but to synthesize the results of several techniques" (Everitt et al., 2011, p. 130) – which is why the table shows all three. When the criteria disagree – and this happens – the decision passes to the researcher: Géron (2022, chapter 9) shows an example in which the silhouette slightly prefers k=4, and he chooses 5 because the clusters are of similar size; and Wendler and Gröttrup (2021, p. 650) note that "the clustering methods do not provide any help in finding an appropriate description of each cluster" – a solution that cannot be described in words is not a good solution, even if its indices are high. The √(n/2) rule of thumb (Mardia et al., 1979, as cited in Wendler & Gröttrup, 2021, p. 649) is too crude for choosing; it appears in the reading box only as a reference point.

How to read the output?

Table 1 – choosing the number of clusters. For each k: W_k, the explained variance, the average silhouette, the gap and its standard error. The row that the gap statistic chooses is marked, and so is the k with the highest silhouette. Figure 1 shows the three curves side by side.

Number of clusters for the solution. Below Table 1, choose k. The default is the gap statistic's choice, and if it chose 1 – the k with the highest silhouette, with a warning. Changing the choice recomputes everything that follows it, without repeating the gap statistic.

Figure 2 – the dendrogram of the first step, with the cut line for the chosen k. The vertical axis is the merge height (in Ward distances): high merges, above a row of low merges, are candidates for natural clusters (Hastie et al., 2009, p. 521). SPSS shows the same tree on a rescaled 0-25 scale ("Rescaled Distance Cluster Combine"; Wendler & Gröttrup, 2021, p. 645) – the shape is identical. In a sample of more than 120 participants the leaves are not drawn one by one: the tree is shown above the 40 groups at its bottom, with the number of participants in each group.

Table 2 – size and quality of the clusters: the number of participants and their percentage, the average silhouette in each cluster and the number of participants with a negative silhouette (who are closer to another cluster), and the size in Ward's cluster before the second step. In the note: how many participants changed cluster in the second step, and the adjusted Rand index between the two steps.

Table 3 – the cluster profile: for each variable – the mean and standard deviation in each cluster and in the whole sample, in the original units, and η² – the proportion of the variable's variance that is between the clusters. Figure 3 shows the same means in standard scores, so that you see at a glance on what each cluster is higher or lower than the mean. From this derives the name you will give each cluster, and that is your interpretation, not a result of the computation (Wendler & Gröttrup, 2021, p. 650).

Figure 4 – the silhouette: a bar for each participant, sorted within each cluster – a "knife shape" for each cluster. The height of the shape is the number of participants in the cluster and its width is their silhouettes; a cluster most of whose members are below the overall average (the dashed line) is a weak cluster (Géron, 2022, Chapter 9). The findings paragraph – in English; the "Interpreting the result" box – do the criteria agree?, the quality of the clusters, what moved between the steps, and what to do now. Downloading the data with the cluster column – the same file with an additional column, for description and validation on the other pages.

Resampling stability – on click. Below the silhouette plot there is a "Compute stability" button (100 pairs and samples by default, with the same seed as the gap statistic). Table 4 gives for each k the distribution of the adjusted Rand index between the solutions of two bootstrap samples – minimum, quartiles, median and maximum – and Figure 5 shows it in a box plot, with a point for each pair. Table 5 moves from the whole partition to the chosen solution: for each cluster – the mean Jaccard similarity to the most similar cluster in the bootstrap samples, how many times it "dissolved" and how many times it was "recovered", and an interpretation according to Hennig (2023). Changing k recomputes only Table 5. Once stability has been computed, it also enters the findings paragraph and the reading box. The methods – in the next section.

Assumptions and warnings

The algorithm always returns clusters. "Such clustering methods always come up with k groups, whatever the data" (Rousseeuw, 1987, p. 55): k-means will return three clusters even from completely uniform data, without any structure; and "hierarchical methods impose hierarchical structure whether or not such structure actually exists in the data" (Hastie et al., 2009, p. 522). That is why the gap statistic is the default for choosing k, and why the reading box warns when the silhouette is weak.

Stability: will the solution replicate in another sample? Leisch (2016, p. 637), in the chapter on stability assessment in the handbook by Hennig and colleagues, distinguishes between two sources of randomness in every partition: the sample – a random subset of the population – and the algorithm, which "is known to impose structure on data, even if they are not well structured" (Dolnicar & Leisch, 2010, as cited there); therefore "using a single partition resulting from a single computation of a particular algorithm puts researchers at risk of using a random, rather than a reliable solution". The two-step method eliminates the randomness of the algorithm (there is no random initialization), but not that of the sample. Stability is especially important when decisions are made based on the partition – in market segmentation it is one of six criteria for evaluating segments (Wedel & Kamakura, 1998, as cited there, p. 639) – and it is not the same as the existence of separated clusters: in Figure 28.1 four k-means clusters on uniform data in a square are very stable, and on a circle – not stable at all (p. 638). The way to check is resampling: draw bootstrap samples from the data, cluster each one, and compare the partitions (pp. 639-640); the chapter presents three methods – prediction strength (Tibshirani & Walther, 2005), per-cluster stability with the Jaccard index (Hennig, 2007), and ARI between pairs of bootstrap samples (Dolnicar & Leisch, 2010; all as cited there, pp. 641-643) – and concludes that "resampling makes a big difference, while the exact method is not so important" (p. 650). k can also be chosen by stability, but with caution: "it often cannot be ruled out that a partition other than the true one also has high stability" (p. 639). The page runs, on click, the third method and the second.

Stability of the whole partition (Table 4 and Figure 5). In the method of Dolnicar and Leisch (2010, as cited in Leisch, 2016, section 28.2.3), two bootstrap samples of the sample size are drawn, each is clustered, all the original participants are assigned to the nearest center in each of the two, and the two assignments are compared with the adjusted Rand index; this is repeated many times, and each pair gives an independent value. The method is "in the middle" between the other two: it measures the whole partition, but allows some clusters, perhaps small ones, to be unstable (p. 642). Look at the distribution, in a box plot, and "values close to one indicate high stability" (p. 643). In Leisch's artificial example, at the correct number of clusters all the indices equal 1, and at a wrong number the distribution shifts to the left, and is sometimes bimodal (p. 645). On the page, the same pairs are used for all values of k, so that the comparison between values of k is on the same samples.

Stability of each cluster (Table 5). In Hennig's method (2007, as cited in Leisch, 2016, Section 28.2.2) the chosen partition is examined, not k in general (p. 646): in each bootstrap sample, take the participants that were sampled, each one once, cluster them, and for each original cluster find the most similar cluster in the sample by Jaccard similarity – the number of shared participants divided by the number of participants in both together. The mean similarity across the samples is the cluster's stability measure. Hennig (2023), in the documentation of the clusterboot function he wrote in the fpc package, gives rules of thumb: a valid, stable cluster should have a mean of .75 or more; between .60 and .75 the cluster indicates a pattern in the data, but which participants belong to it is highly doubtful; below .60 – do not trust it; and "highly stable" – .85 or more. A value of .50 or less in a single sample is considered a "dissolution" of the cluster, and this has a theoretical justification (Hennig, 2008, as cited in García-Escudero et al., 2016, p. 668). And in the same documentation: stability is not the only criterion for validity, and even a very rigid method can give stable clusters that are not valid.

In the sample file (100 pairs, seed 2026) Table 4 shows exactly what the chapter describes. For k = 3 the median ARI is .95 and the lower quartile .94 – stable. For k = 2 the median is .87, but the lower quartile .33 – the distribution is bimodal, that is, two different partitions into two clusters compete with each other. Leisch describes such a distribution as a sign of "several local maxima": a local optimum in the original data can be the global optimum in a bootstrap sample (p. 643). From k = 4 the median declines: .76, .62, and down to .50 at k = 8. In Table 5, the three clusters have mean Jaccard values of .98, .98 and .96 – highly stable, and none of them dissolved in any sample. And another warning from the chapter: k can be chosen by stability, but with caution – "it often cannot be ruled out that a partition other than the true one also has high stability" (p. 639). In the tourism survey in the chapter, the partition into 5 was the most stable, but its segments had no clear profile, so it was not chosen (p. 648); and in the partition into 9 some of the clusters were more stable than others, and the stable ones were the interesting ones (p. 649).

The shape of the clusters. k-means and Ward prefer "spherical" clusters of similar size: "k-means does not behave very well when the clusters have varying sizes, different densities, or nonspherical shapes" (Géron, 2022, Chapter 9, "Limits of k-means"); Ward "tends to find spherical clusters of similar size" and is "sensitive to outliers" (Everitt et al., 2011, p. 78, Table 4.1), "often appears to work well but may impose a spherical structure where none exists" (p. 84), and in simulation studies it succeeded when the clusters were of similar size and failed when they were of different sizes (p. 83); in Ward "the cluster sizes are approximately equal" (Wendler & Gröttrup, 2021, p. 648), and k-means "tends to produce clusters of similar size" (p. 677). A small, truly distinct group may be swallowed up; for elliptical clusters – a Gaussian mixture (the unsupervised learning page, section 4).

Outliers. One extreme participant pulls a center toward it (a mean) and increases the sum of squares; in single linkage outliers remain as single participants, and in complete linkage they are swallowed into clusters (Wendler & Gröttrup, 2021, p. 648, Table 7.24). A cluster of one or two participants in the solution is almost always an outlier, not a type – check on the Outliers and assumptions page.

Practical limitations. Missing values: listwise deletion on the clustering variables (the page reports how many were dropped; for alternatives – Hastie et al., 2009, p. 507, and the Missing values page). Size: up to 3,000 participants, because hierarchical clustering stores all the distances between pairs of participants; the gap statistic runs the whole analysis on each reference set, so with 3,000 participants and 100 samples it takes about a minute (the number of samples can be reduced). The result is deterministic, but when there are identical participants or exactly equal distances (common with Likert items), the order of equal-height merges is determined by the order of the rows, as in R. The variables must be quantitative: a categorical variable with numeric coding (1 = male, 2 = female) enters as a number and makes the distance meaningless – the page warns.

What the page does not compute. Other linkages (single, complete, average), other distances (Manhattan, cosine), categorical and mixed variables (TwoStep with log-likelihood), Gaussian mixtures and BIC, DBSCAN, significance tests between clusters (see the common mistakes), stability with other resampling methods (prediction strength, subsamples, adding noise), and automatic external validation – download the cluster column and check it on the chi-square and analysis of variance pages.

And how do you do it in SPSS?

In SPSS Statistics the two steps are two separate procedures under Analyze › Classify. The paths and options below are general information about the software, as on the unsupervised learning page, section 7; the page's sources deal with SPSS Modeler and not with Statistics. (1) Standardization: Analyze › Descriptive Statistics › Descriptives with Save standardized values as variables – variables named Z… are created. (2) Ward: Hierarchical Cluster on the Z variables; under the Method button – Ward's method with Squared Euclidean distance; under the Plots button – Dendrogram; under the Save button – Single solution with k. The Coefficients column in SPSS's Agglomeration Schedule table is not on the scale of the tree height in R and on the page, but the order of the merges should be the same order, except for equal-height merges. (3) k-means: K-Means Cluster on the same variables, with the number of clusters, and initial centers from a file (Read initial); the centers are computed as the means of the Z variables by Ward cluster (Data › Aggregate). Without initial centers, SPSS chooses them itself from the data, and the result may differ from the page's.

The syntax the page produces performs the three steps: DESCRIPTIVES /SAVE, CLUSTER with /METHOD WARD /MEASURE SEUCLID and /SAVE CLUSTER(k), and QUICK CLUSTER with /INITIAL – the Ward centers that the page computed, embedded as numbers – and with /METHOD KMEANS(NOUPDATE), that is, updating the centers only at the end of each pass over the data, as in the page's algorithm. The syntax was not run in SPSS. In PSPP, free software that reads SPSS syntax, only the first part (standardization) was tested, because it has no CLUSTER command and no /INITIAL. The silhouette, the gap statistic and the resampling stability checks are not computed by SPSS Statistics – for them, R or the page.

And how do you do it in R?

In base R: scale for standardization, hclust(dist(z), method = "ward.D2") for Ward (ward.D2 is Ward on Euclidean distances; the old ward.D assumes you have already squared them), cutree for cutting, and kmeans with centers = the Ward centers and algorithm = "Lloyd" – the page's algorithm (R's default is Hartigan-Wong, which "ensures that there is no single switch of an observation from one group to another group that will decrease the objective"; Hastie et al., 2009, p. 510. In the sample file, for k = 2 and 3 both algorithms give the same solution; for larger k Hartigan-Wong sometimes finds a solution with a slightly smaller W_k – and the page reports in the reading box how many participants' single move would reduce W_k). The same method as a ready-made function: hkmeans in the factoextra package (Kassambara, 2017, pp. 163-164) – Ward (ward.D2) followed by kmeans, but with R's defaults (Hartigan-Wong, and iter.max = 10); in the sample file it is identical to the page for k = 2 and 3, and for larger k gives a slightly smaller W_k. The cluster package gives silhouette and clusGap; with d.power = 2, spaceH0 = "original" and set.seed equal to the page's seed – the same gap and the same standard error as in Table 1. At the end of the script – the two stability checks of Tables 4-5: 100 pairs of bootstrap samples with ARI for each k, and clusterboot from the fpc package with the two-step solution as the clustering method; with the same seed they give exactly the page's values. The script was run as is on the sample file (the values in the comments). The site's rule for Hebrew in R: a UTF-8 file, every comment line starts with an English word, and English variable names (the English file).

# Cluster analysis – Professor Stat. This file is saved in UTF-8; every comment line starts with an English word
library(readxl)    # reading Excel
library(cluster)   # silhouette and gap statistic
library(mclust)    # adjusted Rand index (adjustedRandIndex)
library(fpc)       # stability of each cluster (clusterboot)

# Loading: the standard template – row = participant, column = variable; the sheet "Data"
d <- read_excel("C:/Rdata/DATA_cluster_example_eng.xlsx", sheet = "Data")
v <- c("passive_use", "social_comparison", "content_creation", "social_connection", "daily_hours")
x <- as.data.frame(d[complete.cases(d[, v]), v])           # listwise deletion (148 of 150)
z <- scale(x)                                             # standardization: z-scores (divisor n − 1)

# Step one: Ward on Euclidean distances – the tree, and the height of the last merges
hc <- hclust(dist(z), method = "ward.D2")
rev(tail(hc$height, 5))                                   # heights: 22.80 17.44 7.15 6.99 5.71 – a jump after the third

# Step two: for each k – cutting the tree, the centroids as initial centers, and k-means (Lloyd)
two_stage <- function(z, k, hc = hclust(dist(z), "ward.D2")) {
  w <- cutree(hc, k)
  centers <- apply(z, 2, function(col) tapply(col, w, mean))
  kmeans(z, centers = centers, iter.max = 100, algorithm = "Lloyd")
}
res <- lapply(2:8, function(k) two_stage(z, k, hc))
W   <- c(sum(z^2), sapply(res, function(m) m$tot.withinss))           # W_k, including k = 1
sil <- sapply(res, function(m) mean(silhouette(m$cluster, dist(z))[, "sil_width"]))
round(cbind(k = 1:8, W = W, R2 = 1 - W / W[1], silhouette = c(NA, sil)), 3)   # table: R² .365/.575/.605; silhouette .337/.371/.274 (k = 2/3/4)

# Gap statistic: the same function, 100 uniform reference samples within the range of each variable; set.seed = the page's seed
set.seed(2026)
gap <- clusGap(z, FUNcluster = function(x, k) two_stage(x, k), K.max = 8, B = 100,
               d.power = 2, spaceH0 = "original")
print(gap, method = "Tibs2001SEmax")                      # gap: k = 3; gap 0.404/0.607/0.863/0.814, SE ≈ 0.035

# The chosen solution: sizes, and profile in the original units
km <- res[[2]]                                            # k = 3
km$size                                                   # sizes: 59 51 38
round(aggregate(x, by = list(cluster = km$cluster), FUN = mean), 2)
summary(silhouette(km$cluster, dist(z)))$clus.avg.widths  # silhouette per cluster: .370 .369 .373
d$cluster <- NA; d$cluster[complete.cases(d[, v])] <- km$cluster   # the cluster column for the file

# The same solution in a ready-made function: hkmeans in the factoextra package (Kassambara, 2017, chapter 16) – Ward then k-means, but with Hartigan-Wong
# and without iter.max = 100; in the sample file identical for k = 2-3, and for k = 4 and above W slightly smaller
# factoextra::hkmeans(z, 3)$size                          # 59 51 38

# Resampling stability (Leisch, 2016, chapter 28) – like Tables 4-5 on the page, with the same seed
# (a) The whole solution (Dolnicar & Leisch, 2010, as cited there): two bootstrap samples, clustering in each for each k, assigning all the
# original participants to the nearest center in each of the two, and ARI between the two assignments; 100 pairs – the same pairs for each k
nearest <- function(z, centers) apply(z, 1, function(r) which.min(colSums((t(centers) - r)^2)))
set.seed(2026)
stab <- t(replicate(100, {
  s1 <- z[sample.int(nrow(z), replace = TRUE), ]
  s2 <- z[sample.int(nrow(z), replace = TRUE), ]
  sapply(2:8, function(k) adjustedRandIndex(nearest(z, two_stage(s1, k)$centers), nearest(z, two_stage(s2, k)$centers)))
}))
colnames(stab) <- paste0("k=", 2:8)
round(apply(stab, 2, quantile, c(.25, .5, .75)), 3)       # Table 4: median k=2 .868 (lower quartile .333), k=3 .945 (.939), k=4 .760, k=5 .621 and down to .503 at k=8
boxplot(stab, xlab = "k", ylab = "ARI")                   # Figure 5 (without the points)

# (b) Each cluster (Hennig, 2007, as cited there): clusterboot from the fpc package, with the two-step solution as the clustering function
hkm <- function(data, k) {
  cl <- two_stage(as.matrix(data), k)$cluster
  list(result = cl, nc = k, clusterlist = lapply(1:k, function(i) cl == i), partition = cl, clustermethod = "Ward + k-means")
}
cb <- clusterboot(z, B = 100, bootmethod = "boot", clustermethod = hkm, k = 3, seed = 2026, count = FALSE)
round(cb$bootmean, 3)                                     # Table 5, mean Jaccard: .981 .975 .963
cb$bootbrd; cb$bootrecover                                # dissolved (≤ .50) / recovered (> .75), out of 100: 0 0 0 / 100 100 100

Reporting the results

Cluster analysis involves many decisions by the researcher, so the report should make it possible to reproduce them: which variables entered and why; whether they were standardized; the distance and the method (Ward, followed by k-means from the Ward centers); how the number of clusters was chosen – the criteria and their values, even when they disagree; the size of each cluster; a profile table with the means and standard deviations; the average silhouette; whether stability was checked – the method, the number of samples and the result; and whether validation was done – on what the clusters differ in variables that did not take part in building them. APA (2020, section 6.43) asks to report "sufficient information to help the reader corroborate the analyses". The names of the clusters ("passive consumers") are interpretation, and are reported with the means on which they are based. The formatting rules (APA, 2020, Section 6.36): means and standard deviations with a zero before the decimal point, the average silhouette and η² without. The findings paragraph in the output is built this way, in English.

Common mistakes

Forgetting to standardize variables in different units – and getting clusters of the variable with the large range. And standardizing without thinking – when a noisy variable gets the weight of a discriminating variable (Hastie et al., 2009, p. 506).

Clustering by dozens of single Likert items instead of by scale scores – distance on single items is controversial (Wendler & Gröttrup, 2021, p. 635), and many items of one scale give it multiple weight.

Putting in "everything there is", including background variables (sex, age) and outcome variables – and then "discovering" that the clusters differ on them. A variable that took part in building the clusters cannot validate them.

Choosing k by what looks good without reporting the criteria, or choosing the k with the highest silhouette when it is .20 – that is, when there is no clear structure.

F or t tests between the clusters on the clustering variables, and reporting that "the clusters differ significantly". The algorithm built the clusters so that they would differ precisely on these variables, so significance is guaranteed in advance and adds no information. The page prints η² as a descriptive measure, without p.

Interpreting the dendrogram as "the structure of the data". "Different hierarchical methods, as well as small changes in the data, can lead to quite different dendrograms" (Hastie et al., 2009, p. 521).

Settling for a single run. A single partition from a single sample may be "a random, rather than a reliable solution" (Leisch, 2016, p. 637); before reporting types – check resampling stability (Tables 4-5).

Seeing a cluster as a real group – a type that exists in the population – without validation: a cluster is a description of this sample, with these variables and this distance (the unsupervised learning page, section 6).

References

American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). https://doi.org/10.1037/0000165-000

Everitt, B. S., Landau, S., Leese, M., & Stahl, D. (2011). Cluster analysis (5th ed.). Wiley. https://doi.org/10.1002/9780470977811

García-Escudero, L. A., Gordaliza, A., Matrán, C., Mayo-Iscar, A., & Hennig, C. (2016). Robustness and outliers. In C. Hennig, M. Meilă, F. Murtagh, & R. Rocci (Eds.), Handbook of cluster analysis (pp. 653-677). Chapman & Hall/CRC.

Géron, A. (2022). Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow: Concepts, tools, and techniques to build intelligent systems (3rd ed.). O'Reilly Media.

Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical learning: Data mining, inference, and prediction (2nd ed.). Springer. https://doi.org/10.1007/978-0-387-84858-7

Hennig, C. (2023). fpc: Flexible procedures for clustering (R package version 2.2.11) [Computer software]. https://CRAN.R-project.org/package=fpc

Kassambara, A. (2017). Practical guide to cluster analysis in R: Unsupervised machine learning (Ed. 1). STHDA.

Leisch, F. (2016). Resampling methods for exploring cluster stability. In C. Hennig, M. Meilă, F. Murtagh, & R. Rocci (Eds.), Handbook of cluster analysis (pp. 637-652). Chapman & Hall/CRC.

Rousseeuw, P. J. (1987). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20, 53-65. https://doi.org/10.1016/0377-0427(87)90125-7

Wendler, T., & Gröttrup, S. (2021). Data mining with SPSS Modeler: Theory, exercises and solutions (2nd ed.). Springer. https://doi.org/10.1007/978-3-030-54338-9

Ward (1963), Kaufman and Rousseeuw (1990), Marriott (1982), Milligan and Cooper (1988) – as cited in Everitt et al. (2011); IBM (2015) and Mardia et al. (1979) – as cited in Wendler and Gröttrup (2021); Tibshirani et al. (2001) – as cited in Hastie et al. (2009) and in Everitt et al. (2011); Dolnicar and Leisch (2010), Wedel and Kamakura (1998), Tibshirani and Walther (2005) and Hennig (2007) – as cited in Leisch (2016); Hennig (2008) – as cited in García-Escudero et al. (2016). Hennig (2023) – the documentation of clusterboot (the help page in the package), which was read in the original. According to the site's rule, the list of sources includes only sources that were read in the original. Géron (2022) in a digital copy without page numbers – cited by chapter and section; Everitt et al. (2011), Hastie et al. (2009), Wendler and Gröttrup (2021), Rousseeuw (1987), Leisch (2016), García-Escudero et al. (2016) and Kassambara (2017) – by the printed page. The computation on the page is verified against R 4.3.3 (hclust, cutree, kmeans) and the packages cluster (silhouette, clusGap) and mclust (adjustedRandIndex), and the stability – against R's sample.int and fpc (clusterboot).