Least-squares bilinear clustering of three-way data

A least-squares bilinear clustering framework for modelling three-way data, where each observation consists of an ordinary two-way matrix, is introduced. The method combines bilinear decompositions of the two-way matrices with clustering over observations. Different clusterings are defined for each part of the bilinear decomposition, which decomposes the matrix-valued observations into overall means, row margins, column margins and row–column interactions. Therefore up to four different classifications are defined jointly, one for each type of effect. The computational burden is greatly reduced by the orthogonality of the bilinear model, such that the joint clustering problem reduces to separate problems which can be handled independently. Three of these sub-problems are specific cases of k-means clustering; a special algorithm is formulated for the row–column interactions, which are displayed in clusterwise biplots. The method is illustrated via an empirical example and interpreting the interaction biplots are discussed. Supplemental materials for this paper are available online, which includes the dedicated R package, lsbclust.


Introduction
Multiway data, a generalization of the familiar two-way samples-by-variables data matrix, is becoming more common in a variety of fields. In computer vision applications, images are stored as three-way arrays, or third-order tensors, with rows and columns representing pixel locations and the third way (tubes) representing different color channels (e.g., Krizhevsky et al. 2012). Videos constitute fourth-order tensors, since they capture a sequence of such images over time (e.g., Abu-El-Haija et al. 2016). In the social sciences, a marketing research survey asking multiple individuals to rate several products on various characteristics using a Likert scale generates a three-way array of rating scores (e.g., DeSarbo et al. 1982). Similar research designs are commonly used in sensometrics (e.g., Cariou et al. 2021). Other applications include: high-throughput molecular data in bioinformatics (e.g., Lonsdale et al. 2013); spectroscopic data in chemometrics (e.g., Faber et al. 2003;Bro 2006); and neuroimaging data, collected using electroencephalography (EEG) or functional magnetic resonance imaging (fMRI), for example (e.g., Genevsky and Knutson 2015).
A taxonomy of measurement data is given by Carroll and Arabie (1980), where a mode is defined as "a particular class of entities" associated with a data array, such as stimuli, subjects or scale items; a three-way array may have up to three modes. Kiers (2000) introduced standardized notation and terminology for multiway analysis, while Kroonenberg (2008) is devoted to multiway data analysis methodology.
Performing cluster analysis on one or more of the modes of a three-way data array X is an important problem in multiway data analysis. Statistical research into three-way clustering can be broadly divided into two streams. The application of finite mixture models for clustering three-way data has been studied in, amongst others, Basford and McLachlan (1985), Hunt and Basford (1999), Vermunt (2007), Viroli (2011), Meulders and De Bruecker (2018), Gallaugher and McNicholas (2020a, b). Whereas these rely on maximum likelihood estimation, another research stream have focused on nonparametric data approximation methods mainly via least-squares estimation. The papers of Vichi (1999), Rocci and Vichi (2005), Vichi et al. (2007), Papalexakis et al. (2013), Wilderjans and Ceulemans (2013), Llobell et al. (2019), Llobell et al. (2020) and Cariou et al. (2021) belong to this stream, as does our paper. An approach often used in this stream is to generalize a three-way data decomposition, such as Tucker's three-mode factor analysis (Tucker3; Tucker 1966) or the Candecomp/Parafac decomposition (CP; Carroll and Chang 1970;Harshman 1970;Hitchcock 1927), by adding clustering over one or more of the modes.
Here we propose a new method-called least-squares bilinear clustering (or lsbclust)-derived in a similar manner but from a different starting point. In contrast to starting with a three-way decomposition, we start with a bilinear (or biadditive) decomposition of the matrix slices X i (i = 1, . . . , I ) of X for the mode (by convention, the first mode) over which we want to cluster (e.g., Denis and Gower 1994). An example of a bilinear approximation of a matrix X i is where 1 denotes a column vector of ones. The terms on the right-hand side of (1) can be interpreted as the grand mean (m), row main effects (a), column main effects (b) and row-column interactions (cd ), respectively, and are subject to appropriate identification restrictions. By adding one or more sets of clusters to the effects on the right-hand side of Eq. (1), and by decomposing all matrix slices at the same time, the lsbclust method ensures that the estimated effects are equal for observations in the same cluster. The choice of bilinear decomposition has two interesting features. First, it allows the results to be displayed graphically using biplots (Gower et al. 2011;Gower and Hand 1996;Gabriel 1971), and other standard graphs. Second, it provides the possibility to have a different set of clusters for each of the four types of effects. This can be useful in cases where interesting differences can be expected not only with respect to the row-column interactions, but also with respect to the grand mean and main effects on the right-hand side of (1). We will discuss this aspect further in Sect. 3.
The main contribution of this paper is the introduction of the lsbclust model and loss function, together with an algorithm for estimating the parameters. We also provide an implementation of this algorithm and its related graphical procedures in an R (R Core Team 2020) package lsbclust (available on the Comprehensive R Archive Network, or CRAN), while supplemental materials illustrate its use. The method can be used as a complement to or replacement of other three-way data analysis procedures.
The remainder of the paper is structured as follows. Section 2 introduces the basic model and loss function used, which Sect. 3 augments with clustering to arrive at the lsbclust formulation. Section 4 shows how to simplify the loss function, and Sect. 5 discusses our algorithm for estimating the parameters. Enhancements to the basic method are discussed in Sect. 6, including using different bilinear models and aspects of biplot construction. The results of a simulation study is reported in Sect. 7. An empirical example is presented in Sects. 8, and 9 concludes.
Before we introduce lsbclust in more detail, we provide more detail on related methods in Sect. 1.1.

Related methods
The most prominent multiway data analysis methods, which include the tucker3 and CP decompositions, seek to generalize the singular value decomposition (SVD) of a matrix to the multiway case. Kolda and Bader (2009) provide an excellent review of these tensor decomposition methods and their applications across diverse fields. De Silva and Lim (2008) discusses mathematical aspects of generalizing the Eckart-Young theorem (Stewart 1993;Eckart and Young 1936;Schmidt 1907) to higher-order tensors, while Kiers and Van Mechelen (2001) discuss and illustrate the application of the tucker3 decomposition.
Three-mode factor analysis (Tucker 1966;Kroonenberg and de Leeuw 1980), commonly referred to as Tucker3, is the most general widely-used three-way method. Let X J ,K I = X 1 X 2 · · · X I , and let the component matrices B, C and D be lowdimensional columnwise orthonormal configurations for the first, second and third ways of X with dimensions P, Q and R respectively. Also, let H : P × Q × R be the so-called core array, which gives the interactions between the elements of the component matrices. The Tucker3 model has loss function where b i p denotes an element of B, H p : Q × R is a slice of H along the first mode, · is the Frobenius norm and ⊗ is the Kronecker product. Interpreting a Tucker3 solution involves interpreting the three component matrices and the core.
No clustering is performed, and (2) is typically minimized by alternating least squares (ALS; Kroonenberg and de Leeuw 1980). The CP decomposition (Carroll and Chang 1970;Harshman 1970) is a restricted version of Tucker3 where P = Q = R and H is replaced by E. Here E contains elements e pqr = 1 if p = q = r and e pqr = 0 otherwise. This implies that each component (column) in B is related to a single component in C and D, and vice versa, while in Tucker3 all components in one mode are related to all other components in the other modes. The CP loss function is Here diag(b i ) is the diagonal matrix with the ith row of B on the diagonal, and E P,P P = E 1 E 2 · · · E p with E p being the pth matrix slice of E.
Neither of these methods employ clustering along any of the modes, but Rocci and Vichi (2005) formulate such a variant of Tucker3, namely T3clus. The idea is to replace B by the indicator matrix G : I × U , which simply indicates which of U clusters each of the I observations belongs to. The core array H now represents the three-way array of cluster centroids in the reduced component space. The T3clus loss function is Another related approach is the clusterwise CP method of Wilderjans and Ceulemans (2013)-CPclus hereafter. The loss function for this approach is In contrast to T3clus, here the component matrices are cluster-specific. Table 1 gives a summary of the loss functions of these three-way decompositions, together with the lsbclust loss function for the interaction effects (where J denotes See the text for references and an explanation of the notation used. For lsbclust, we include only the part of the loss function which relates to the row-column interactions. Here J denotes a centring matrix of the appropriate size a centring matrix of the appropriate size). From the table it is clear that each method approximates the matrix slices X i (i = 1, . . . , I ) in different ways. Tucker3 and Candecomp/Parafac are symmetric with respect to all three modes, but do not involve clustering. In case all X i are double-centred, T3clus approximates X i in cluster u by CH u D , as compared to C u D u for lsbclust. T3clus therefore requires the clusterwise approximation of X i to lie in the same row and column subspaces across clusters, while lsbclust has no such restrictions. The lsbclust loss function is clearly a constrained version of CPclus which replaces diag(b i ) with the identity matrix. This reduces the number of parameters and simultaneously enforces the same interpretation for all members of the same cluster. The additional constraint means that a single biplot can be used to interpret the interaction effects in a cluster, instead of requiring one plot per cluster member. In the special case where X contains indicator matrices for each subject, lsbclust has a strong connection to the latent-class bilinear multinomial logit (lcbml) model of van Rosmalen et al. (2010). Their model was developed specifically to deal with response styles when analyzing two-way self-report survey data. More details on the lc-bml model, and its connection to lsbclust, are given in the supplemental materials.
The next section discusses the basic building blocks of lsbclust.

Basic model and loss function
Consider real-valued data collected in a three-way array X. For the marketing research survey example, the entry x i jk of X is the rating by person i of product j on characteristic k, with i = 1, . . . , I ; j = 1, . . . , J ; and k = 1, . . . , K , respectively. Let X i denote the J × K matrix of responses for subject i, which is also the ith (frontal) slice of the array X. By convention, we single out the first way of X by treating the X i as observations, but the method could also be applied to the second or third way of X. We derive the proposed lsbclust formulation by augmenting (with clusters; Sect. 3) the following least-squares loss function: Here · is the Frobenius norm. Moreover, 1 K is the length-K vector of ones; below we will drop the subscript-as in 1-for simplicity. Equation (6) approximates each X i using the same bilinear (or biadditive) model. This model contains a grand mean m, the row and column main effects a and b respectively, and row-column interactions CD .
The matrices C and D are low-dimensional representations of the rows (products) and columns (characteristics) of the X i matrices respectively, but after adjusting for main effects. Representing the interaction effects using such inner products permit these to be displayed in biplots if the dimensionality of C and D are low enough so that displays can be made (Gower et al. 2011;Gower and Hand 1996;Gabriel 1971). Therefore, the dimensionality is typically set to two, although values up to min{J , K } − 1 are possible. To ensure uniqueness of the model, the usual sum-to-zero constraints a 1 = b 1 = 0 and 1 C = 1 D = 0 must be imposed. Additionally, the columns of C and D are required to be orthogonal (for more information, see Denis and Gower 1994, as well as Sect. 6.3). Model (6) has an analytical solution.

Capturing heterogeneity with clusters
Equation (6) applies the same parameters to the entire data set. To capture potential heterogeneity, we allow for the existence of a prespecified number of clusters in a single mode between which the parameters in the bilinear decomposition of the X i matrices may vary. Moreover, since the bilinear decomposition itself has four different sets of parameters, we introduce four sets of clusters-one for each type of parameter. This allows the model to recognize that, for example, although X i and X i (i = i ) are similar with respect to main effects, they could differ in the interaction effects (or vice versa).
Let G (o) be the I × R matrix of cluster memberships for the grand mean effect m, which has g (o) ir = 1 if person i belongs to cluster r and g (o) ir = 0 otherwise (r = 1, 2, . . . , R). Similarly, G (r) is the I × S matrix of cluster memberships for the row effects, G (c) the I × T matrix of cluster memberships for the column effects, and G (i) the I × U matrix of cluster memberships for the interaction effects. Now, by incorporating the clustering, the least-squares loss function becomes with the summation over all observations (I ) and clusters (R, S, T and U ).
Note that it is assumed implicitly that all low-rank decompositions are of the same rank P for all clusters. For identifiability, the same sum-to-zero constraints now apply to each cluster-specific set of parameters, i.e., D u 1 = 0 for all u = 1, . . . , U .
The columns of all C u and D u matrices are orthogonal. Equation (7) allows for a total of RST U clusters by clustering each X i on four different types of effects at the same time. However, we show next that the joint clustering problem can be simplified into four separate clustering problems, significantly reducing the computational complexity since only R + S + T + U clusters will need to be found.

Decomposing the loss function
To simplify the exposition, note that mathematically we can drop the sum-to-zero constraints from the formulation by introducing centering matrices of the form for some positive integer c controlling the size of the matrix (which we duly omit below for simplicity). This is done by redefining the terms in the summation in (7) as such that the sum-to-zero constraints on the columns of C u and D u , and on a s and b t , are automatically enforced. For example, estimating the parameters in Ja s is equivalent to estimating a s subject to 1 a s = 0. To simplify the notation, we redefine A = [Ja 1 · · · Ja S ] to avoid writing JA. The matrices B, C and D are also redefined analogously. We proceed to simplify (9) by first expanding the double-centred JX i J into separate terms. Then, by adding two additional centering operators for the row and column means, X i can be rewritten as We can now associate each of the terms in (10) with the corresponding terms in the model (9): It can be shown (see the appendix in the supplemental materials) that the decomposition (11) is orthogonal such that This equality follows from the fact that all the cross-products are zero. Importantly, the orthogonality leads to a great simplification in the clustering, since now the loss function (7) equals The main consequence of (13) is that the joint clustering reduces to separate clusterings on the grand means, row main effects, column main effects and interactions respectively. This implies that each of these four parts can be treated independently under this loss function and model. It also gives mathematical justification for the procedure whereby all X i are first double-centred to remove the overall, row and column margins, and then analyzed by minimizing L (i) G (i) , C, D to study the rowcolumn interactions. If the researchers are interested in the grand mean, row or column marginal effects, these can be analyzed separately by minimizing the corresponding loss functions. These are frequently omitted when only the row-column interactions Algorithm 1: ALS algorithm for minimizing L (i) G (i) , C, D 1. Set k = 0 and 0 = I . Randomly initialize G (i) as G 0 . 2. While k > 0: (a) Compute cluster sizes I 1 , . . . , by assigning all observations to the closest (approximated) cluster mean: (d) Compute the number of reassignments 3. Output G k for G (i) , as well as C and D.
are of interest; otherwise, the researcher may prefer to jointly cluster on more than one of these effects at the same time (see Sect. 6.1).
Next, we discuss computational aspects of the proposed lsbclust method.

Estimation algorithms
Due to the form of the loss function (13), we can treat each of the components separately. Conveniently, the loss functions are specific instances of the well-known k-means loss function (e.g., Everitt et al. 2011). They differ only with respect to the data matrix Y : I × d (say) on which k-means cluster analysis is to be applied, which can respectively be defined as follows: Hence optimizing L (o) G (o) , m , L (r) G (r) , A and L (c) G (c) , B can resort to standard methods for k-means on the overall mean, row margins and column margins respectively. Also, there are a variety of tools available for selecting R, S and T .
Minimizing L (i) G (i) , C, D however requires a custom algorithm. The main steps of our algorithm, which is based on block-relaxation methods (see, for example, de Leeuw 1994), is outlined in Algorithm 1. It iterates over optimizing C and D in Step 2b while keeping G (i) fixed, and vice versa in Step 2c. This approach is also known as alternating least squares (ALS).
Convergence is reached when no reassignment of a single observation to a different cluster will reduce the value of the loss function. This corresponds to k = 0 in Algorithm 1. The algorithm is guaranteed to converge monotonically, but only to some accumulation points which are usually local minima. It must be initialized by a starting configuration for G (i) . To increase the likelihood of locating the global minimum, it is advisable to use multiple (random) starting values for G (i) .
We now describe the derivation of the key steps of Algorithm 1.
Step 2b, where L (i) G (i) , C, D is minimized over C and D for fixed G (i) , relies on the decomposition Here iu is the cardinality of cluster u, and X u = 1 Since only the final term in (14) depends on C and D, it is sufficient to minimize this term only. Let the singular value decomposition (SVD) of JX u J be U u u V u , where U u and V u are orthonormal and u diagonal. The Eckart-Young theorem (Eckart and Young 1936;Schmidt 1907) establishes the best rank-P least-squares approximation of JX u J as the truncated SVD U u u L P V u , where Multiplication by L P sets all singular values except the first P equal to zero. Consequently, we can update C and D using where α denotes the diagonal matrix containing the singular values to the power α.
The parameter 0 ≤ α ≤ 1 is typically taken to be 0.5, but can be set by the user to improve the interpretability of the graphical output. See Sect. 6.3 for a description of the heuristic rule used in our implementation. Now in Step 2c, G (i) is updated while regarding C and D as fixed. The updated G (i) is constructed by simply assigning each i to the cluster with the closest mean, hence minimizing L (i) G (i) , C, D for each individual in a greedy manner. This entails assigning observation i to cluster A few details remain. When an empty cluster is encountered, that cluster is reinitialized using the worst-fitting observation across all clusters. Label-switching is countered by reassigning cluster labels between iterations when necessary. Finally, we prefer reporting the minimized value of L (i) G (i) , C, D after division by In the next section, we briefly discuss alternative model formulations.

Enhancements
Here we first consider accommodating alternative bilinear model specifications in the lsbclust loss function (7). Thereafter, we discuss restricted model formulations for the interaction effects that reduce the number of parameters and aid interpretation. Finally, we briefly consider scaling options used in biplot construction. Some practical guidelines on performing model selection are provided in Sect. 8, where an application is discussed.

Using other bilinear decompositions
In certain situations, a different bilinear model may be more appropriate than the one used in the loss function (7). For example, the researcher may not want to treat the grand, row and column means separately, in which case a more appropriate loss function would simply be with the only constraint required being orthogonality of the columns of C u and D u . Yet a situation may arise where the row means (but not the column means) are themselves interesting, leading instead to the following loss function: Here the centering matrix enforces sum-to-zero constraints on the columns of D u , which allows a s to be estimated. In fact, a variety of different bilinear models can be specified by dropping one or more constraints, and hence cluster sets, from the formulation (7). An exhaustive list of the nine possibilities are given in Table 2, with Model 9 being the original formulation in (7). With the exception of Model 6, these models are orthogonal as in Sect. 4. Hence the algorithms in Sect. 5 also apply for these bilinear models (except Model 6), after minor adjustments to account for the different centering options.
Note that Model 6 is not orthogonal and is only included for completeness In Table 2, and in our software, we characterize the different models using a vector of four binary indicators, δ = (δ 1 , δ 2 , δ 3 , δ 4 ). Each of these correspond to the presence or absence of one of the centering matrices in the model formulation. Specifically, δ 1 = 1 and δ 2 = 1 indicate centering the columns of C u and D u respectively (as in JC u and JD u ), while δ 3 = 1 and δ 4 = 1 in addition implies respectively centering the column and row means-as in Jb t and Ja s . Additionally, it is only possible to have δ 3 = 1 if δ 1 = 1; an equivalent relationship holds between δ 4 and δ 2 .

Common row or column coordinates
The formulation in (7) can contain a large number of parameters. This is mainly because the interaction approximations C u D u require the row and column representations, C u and D u respectively, to be different for each cluster. We can counter this by restricting either the rows or columns to have a common representation across clusters, which has the added benefit of making the biplots based on C u and D u easier to interpret. We see this interpretability as the main reason to elect such a constraint; in a practical application this must be weighed against the resulting reduction in goodnessof-fit.
Consequently, we allow the following three options for modelling the interactions: (I) C u D u : both C u and D u are specific to the interaction cluster (as above); or (II) C 1 D u : a common row representation C 1 for all interaction clusters, and a differential column representation D u for each interaction cluster; or (III) C u D 1 : differential row representations C u but a common column representation D 1 .
If the alternative specifications (II) or (III) are used, C u D u in (7) should be replaced by C 1 D u or C u D 1 respectively. These restricted formulations also require adjustments to Algorithm 1. To facilitate this, we define the following block matrices by stacking either row-or column-wise: The final term in (14) can then be rewritten in the respective formulations as Step 2b in Algorithm 1 is subsequently amended to perform a single SVD on either JX (c) to update C 1 and D * under (II), or on X (r ) J to update C * and D 1 under (III). From these, updates to JD u or JC u are derived for u = 1, . . . , U by extracting the relevant block matrices from D * or C * respectively.

Biplot interpretability
Under case (I), where there is no requirement for the interaction decompositions to be similar across clusters, it can aid interpretation to rotate the configurations so that the biplot axes lie more or less in the same direction. For any orthogonal matrix Q u , it holds for the inner product matrices that C u D u = C u Q u D u Q u , and hence these are invariant to orthogonal rotations. The problem of finding orthogonal matrices Q u , u = 1, 2, . . . , U , such that either the row or column configurations match each other as closely as possible is known as the generalized orthogonal Procrustes problem (Gower 1975;Gower and Dijksterhuis 2004). We apply this by default in our software implementation.
Besides this adjustment, two types of scalings can be used to make the biplot displays more attractive, namely so-called αand λ-scaling. First, since our choice of α in (16) does not change the inner product approximations, we are free to choose it such that the resulting biplots are easy to interpret. In our software implementation we use as a heuristic method the value of α which maximizes the minimum Euclidean distance over all row and column points to the origin. Alternatively, the user can choose any other quantile of these distances, such as the median, or specify the desired value of α explicitly.
Second, note that for matrices A and B it holds that AB = (λA)(B /λ), so that λ can also be freely chosen. Following Gower et al. (2011, Section 2.3.1), we choose λ such that the average squared Euclidean distances from the two sets of points represented by the rows of the matrices in (16) to the origin are equal. For case (I) in (16), for example, this amounts to choosing The appendix in the supplemental materials provides information on goodness-offit indices to quantify the quality of the within-cluster approximations.
Next, we report the results of a stimulation study.

Simulation study
A simulation study was conducted to assess cluster membership recovery of lsbclust. Since the separate clustering problems for the overall mean, row means and column means are simply k-means problems, we focus on reporting the results for the row-column interactions. Simulation studies for ordinary k-means can be found, for example, in Milligan (1980). The results of lsbclust are compared to that of T3clus, as well as to two different versions of k-means clustering. The first version of k-means, which we will call vectorized k-means or vecKmeans, merely vectorizes the matrix slices by stringing them out as long vectors and then applies ordinary k-means to the matrix having these vectors as rows. The second variation first obtains the best P-dimensional approximation to each matrix slice using the SVD, and then applies vecKmeans. This variant we call dimKmeans. Both these methods are to be compared to the lsbclust interaction clustering, hence the matrix slices X i are double-centred before being submitted to these procedures.
Simulated data from the lsbclust model are constructed according to the following steps: 1. Simulate the cluster membership matrix G (i) given the required number of clusters U and the number of observations I . The proportion of observations attributed to each class, π u , u = 1, . . . , U , must be specified. Similarly, generate G (o) , G (r) , G (c) with R, S and T clusters respectively. The same cluster membership probabilities are used for these clusters. 2. Simulate overall means m r , r = 1, . . . , R, from a standard normal (Gaussian) distribution, as well as row means a s , s = 1, . . . , S, and column means b t , t = 1, . . . , T . These are also drawn from standard normal distributions, and subsequently centred. 3. Simulate matrices X u , u = 1, . . . , U , representing the cluster means in (14), as follows: (a) Generate two random orthogonal matrices using Stewart (1980)'s method, representing the rows and columns in P-dimensional space. Column-centre both these matrices to arrive at U u and V u , say. If case (a) or (b) applies, either the same U u = U or V u = V is used for all clusters.
(b) Generate P singular values as random numbers on [0. 5,5], and sort these in decreasing order in the vector γ u . If case (a) or (b) applies, the same singular values γ u = γ are used for all clusters. (c) Construct u = diag γ u and consequently X u = U u u V u .
4. Finally, construct the simulated X by using the relevant simulated cluster means of each type for each matrix slice, and by adding random noise simulated from a normal (Gaussian) distribution with zero mean and standard error σ . Here σ controls the signal-to-noise ratio: larger σ makes it harder to estimate the parameters used in the simulation steps.
In our simulation study, the following factors were varied: • This translates to a 3 2 ×2 3 design, with 72 conditions. For simplification, the following were kept fixed: the model simulated from, namely model (9)-see Table 2; the size of the matrix slices, J = K = 8; and the number of clusters, R = S = T = U = 5.
For each of the 72 conditions, we generate 100 parameter sets. In turn, for each of these parameter sets, we generate 50 randomly sampled data sets, resulting in 5000 simulated data sets per condition, or 360,000 in total. The results can be assessed both on clustering quality and estimation accuracy, where the latter includes clustering quality as well as parameter recovery. Here we discuss clustering quality only.

Clustering quality
Cluster recovery is measured by the adjusted Rand index (ARI; Hubert and Arabie 1985), which in our case quantifies the similarity between the actual, simulated clustering and that recovered by an algorithm. It improves on simple cluster agreement by adjusting for the chance of a randomly chosen pair of subjects falling in the same class. The ARI takes a value of one when the cluster recovery is perfect, and zero when the algorithm performs similarly to random assignment. The ARI can also take negative values, which indicate worse performance than random assignment. We first report the performance profiles of the ARI, which assesses how well each method performs relative to all the other methods. Thereafter, we consider the absolute performance of the methods. to express the performance of each algorithm relative to that of the best performing algorithm on each particular data set, and then to plot the cumulative distribution of this relative performance. Denote by D the set of data sets and by C the set of algorithms. Suppose that p d,c is the ARI for d ∈ D and c ∈ C. The performance ratio r d,c is then defined as the ratio of the best performing ARI on data set d to p d,c : Typically, the best performing method has a performance ratio of one, with other methods having larger performance ratios, indicating how close these methods came to the best method. However, for the ARI, there is one caveat: the ARI may be negative or even exactly zero, in which case (21) do not work as intended. We circumvent this problem by adding a small positive constant to both the numerator and denominator in (21).
The performance profile for algorithm c is simply the empirical cumulative distribution function (ecdf) of the performance ratios. This can be calculated as where | · | denotes the cardinality of a set. Therefore, P c (ν) is simply the empirical probability of algorithm c having an performance ratio of at most ν, and by extension P c (1) is the empirical probability that algorithm c achieves the best performance. We calculate performance profiles for each combination of the five factors manipulated in the simulation study. The results are shown in Figs. 1 and 2. Both of these figures pertain to P = 2 dimensions, with Figs. 1 and 2 purporting to data sets with I = 100 and I = 500 observations respectively. For the sake of brevity, we omit the figures for P = 5 dimensions: in these cases, all algorithms achieve close to optimal results.
The results in Figs. 1 and 2 show that lsbclust generally outperforms the other methods, since its performance profiles in general are larger than that of the other methods. The next best method is vecKmeans, with T3clus coming in last. It should be expected that lsbclust should perform well, since it was used to generate the data. vecKmeans is quite competitive when interaction decomposition (c) is used (where neither component matrices are fixed across clusters). When the actual model do include restrictions on the interaction decomposition, lsbclust performs much better than the other methods.
In terms of factors manipulated in the study, the error standard deviation (σ ), the sample size I and the interaction decomposition is most important. Whether the clusters are balanced or not has very little bearing on the results.
Having assessed the relative performance of the methods, we turn our attention to the absolute ARI achieved on the simulated data. Table 3 reports the average ARI for the different methods. We first calculate the average ARI for the 50 data sets generated for each set of parameters, and then average Table 3 The mean adjusted Rand indices for the different methods, averaged over data and parameter sets The results here considers only P = 2 dimensions and balanced clusters: for P = 5 clusters, all methods achieved nearly optimal ARI on most data sets, there is little difference in results for unbalanced cluster sizes the resultant 100 average ARI's over the 100 different parameter sets. The table only includes results for P = 2 dimensions, and for the case where the cluster sizes are balanced. The latter does not affect the results significantly, so has been omitted. The former does have an important but uninteresting effect: when P = 5, all methods achieved near optimal ARI. Overall, an increase in the error standard deviation degrades the model performance the most. With enough samples, lsbclust can however still achieve decent clustering performance when σ is high.

Adjusted Rand index
In the next section, we consider an illustrative empirical example.

Application
The data comes from a brand positioning study where 187 consumers evaluated 10 car manufacturers on a set of 8 attributes (Bijmolt and van de Velden 2012). It was collected via an online survey using a representative Dutch sample from the Cen-tERpanel of Tilburg University in the Netherlands, with the aim of studying how consumers perceive different car brands. The data X is therefore of size 187 × 10 × 8. The 10 international car brands considered were: Citroën, Fiat, Ford, Opel, Peugeot, Renault, Seat, Toyota, Volkswagen and Volvo. Respondents rated each of these brands on 8 different attributes using a 10-point rating scale. For 6 out of the 8 items, namely Affordability, Attractiveness, Safety, Sportiness, Reliability and Features, a score of 10 is the most desirable outcome. However, for the items Size and Operating Cost, a score of 10 reflects small cars and those with high operating costs, respectively. Hence, higher ratings on these two attributes reflect increasingly negative sentiment, in contrast to the remaining six items.
We fit an lsbclust model with δ = (1, 1, 1, 1)-Model 9 in Table 2-so that the overall means, row means, column means and interactions are estimated separately. Also, we select P = 2 dimensions for display purposes, and we fix the coordinates of the 10 car brands across all interaction biplots-using case (II) from Sect. 6.2-to aid with interpretation.
Besides these choices, the number of clusters for each of the four components must be determined. For simplicity, we opt to do this separately for each of the four sets of clusters, although it is possible to make a joint selection. Selecting the number of clusters is a common problem, and many procedures and criteria have been proposed in the literature. Milligan and Cooper (1985), Hardy (1996) and Everitt et al. (2011) provide an assessment of some of these criteria and additional references. The simplest approach, and the one we use here for illustration, is probably the scree test (Cattell 1966). This method involves running the algorithm for several values of k and plotting the loss function against k. The user must then choose a value for k based on this so-called scree plot, such that the chosen k is close to an "elbow" in the plot. This indicates that adding additional groups to the analysis does not significantly increase how well the results describe the data. Here we fit lsbclust models for 1 to 15 clusters and inspect the resulting scree plots to select R, S, T and U . 1 Based on these plots, we selected R = 5, S = 5, T = 6 and U = 8 clusters. The number of row clusters was reduced to S = 5 from an initial choice of S = 8 to avoid clusters contain a single observation only. We note that these choices are subjective, should take into account the aims of the research and that alternative selection criteria can also be used. The number of random starts used for the interaction and k-means clustering were 100 and 1000 respectively.
The cluster sizes and centres (i.e., m r in Eq. (7)) for each of the five overall mean clusters are shown in Table 4. Noticeable here are that the small clusters O4 (8 observations) and O5 (3 observations) identified persons who invariably used very high and very low scores, respectively. These respondents obviously do not provide very interesting information in their answers. But since their corresponding row means, column means and interactions do not differ from the overall mean, they are merely assigned to the row, column and interaction segments containing negligible effects (see below). lsbclust has therefore been able to identify the 11 persons in clusters O4 and O5 who provide very little sensible information. Figure 3 displays the means of the eight car brand (row) clusters across all attributes. Effect sizes can be read off on the horizontal axis. The first two clusters, Segments R1 (98 observations) and R2 (59 observations), contain no pronounced large effects, which indicates that these consumers do not have strong, consistent opinions on any of the brands across all attributes. The 23 observations in Segment R3 do assign lower scores to Volvo and somewhat higher scores to Peugeot across all items. These effects amount to (approximately) -1.4 and 0.7 for Volvo and Peugeot, respectively. Even larger effects are observed in the smaller segments R4 (4 observations) and R5 (3 observations), but these comprise a relatively small part of the data.
The attribute (column) mean effects for all six clusters are displayed in Fig. 4. Here all segments contain at least 13 observations. Again, the largest cluster, Segment C1 (69 observations), contains no large effects. Segments C2 (35 observations), C3 (33 observations), and C4 (22 observations) are similar in that they display an inclination to assign higher scores on Reliability and Safety, and lower scores on Affordability, irrespective of the car brand being assessed. The magnitude of these effects vary greatly over these clusters though, and Segments C2 and C4 also display negative effects for Size. Segment C5 (15 observations) is dominated by a tendency to assign low scores  The most interesting results can be found among the interactions, which is where respondents distinguish between different car manufacturers on the measured attributes. Figs. 5 and 6 shows the biplots for the eight interaction segments. The car manufacturers are represented by points, and the attributes by arrows. The labels, points and arrows are shaded according to their goodness-of-fit, with well-fitting points being darker. All car brands, except Ford with a fit of only 0.09, fit reasonably wellsee Table 5. The locations of the car brands are fixed across all biplots to make them easier to interpret. It is immediately apparent that the French manufacturers (Peugeot, Citroën and Renault) are judged to be similar, while the German brands Opel and Volkswagen are also located close together. Volvo, the Swedish car manufacturer, is somewhat isolated towards the right of the biplots, in contrast with Fiat at the opposite side of the plot. Fiat and Toyota are judged to be somewhat similar to the French and German brands respectively. Seat in turn are most similar to Toyota. As a result of its low fit, Ford is hardly visible and lies near the origin.
The fit for the eight attributes vary per segment, and is summarized in Table 6. Typically only a subset of items fit well in each segment, and only the better-fitting ones are adorned with calibrated axes in Figs. 5 and 6. For any manufacturer, the estimated within-cluster interaction effects can be read off from the orthogonal projection of its representing point onto the respective biplot axes. For example, Volvo Biplots for the interaction clusters I1 to I4 detected in the cars data. The relative cluster sizes are displayed in the titles of the different panels. Each attribute is represented by a vector, and those which fit best in each panel also by calibrated axes. The car manufacturers are represented by points, which has identical locations across panels. The orthogonal projections of the car manufacturer points onto the attribute axes give the estimated mean effects. The colors and labels are faded according to how well they fit into the display: solid colors fit well and transparent ones fit badly scores approximately 0.25 points above that predicted by the overall mean, row mean and column mean on Sportiness in Segment I1, and Fiat score about 2 rating points above the overall and marginal effects on Affordability in the Segment I6. The overall variance accounted for in approximating the cluster means is 82.8% in two dimensions, with 64.1% and 18.7% attributed to dimensions 1 and 2, respectively. Hence two dimensions are a reasonable choice.
The effects in the interaction clusters should be interpreted as deviations from that expected from the overall and marginal means alone. Here are summaries of the interaction segments: • Segment I1 (32.6%) has no large effects. Hence for this sizable group of individuals, the overall and marginal means contain most of the information provided. • Segment I2 (12.3%) interprets Fiat and Seat to be more affordable than expected from the overall and marginal effects alone. The opposite applies to Volvo and Volkswagen. In terms of safety and reliability, however, the roles are reversed: Fiat and Seat score lower than expected, while Volvo and Volkswagen excel on these items. The effect sizes are roughly the same on the aforementioned three items. Fig. 6 Biplots for the interaction clusters I5 to I8 detected in the cars data • Segment I3 (12.3%) contrasts positive interaction effects on Affordability and Size with negative effects on Safety, Reliability, Attractiveness and Features, and vice versa. Taking into account that Size is reverse-coded, this segment considers Fiat and Seat to be more affordable but smaller cars. At the same time, they are considered less reliable, attractive and safe than especially Volvo and Volkswagen. This segment is not dissimilar to Segment I2, except that effect sizes are larger and the inclusion of Attractiveness and Features on the right side of the plot. • Segment I4 (10.7%) contains smaller effects than Segment I3. It distinguishes somewhat between Affordability and Size, in contrast to Segment I3. Whereas Fiat still better than expected in terms of Affordability, Seat is perceived to perform worse than Fiat on Size (after adjusting for the overall and marginal means). There is also now a much bigger difference between Volvo and Volkswagen in terms of size, with the latter scoring worse than expected compared to Renault, Peugeot and Citroën. • Segment I5 (10.2%) again interprets the left of the plot as better in terms of Affordability and worse in terms of Size. But it also associates this with lower operating costs (since the latter is reverse-coded). Volvo and Volkswagen are scores higher than expected in terms of safety, in contrast to Seat and Fiat. • Segment I6 (8.6%) contains large interaction effects on three pairs of items.
Renault, Citroën and Peugeot score higher on Attractiveness and Features than expected. This is in contrast to Seat and Volkswagen, who exhibit negative interaction effects with these items. Safety and Reliability again favour Volkswagen and Volvo, while Affordability and Size have similar but larger effects than in Segment I4. • Segment I7 (7.5%) have similar interaction effects for Affordability and Size as with Segments I4 and I5. But these individuals also consider Volkswagen, Volvo and Opel to have higher than expected operating costs, in contrast to Fiat, Renault and Citroën. The latter brands, together with Peugeot are also considered the most attractive, with Volkswagen being considered less attractive than expected using only the overall and main effects. • Segment I8 (5.9%) contains individuals who consider Renault, Fiat, Citroën and Peugeot to be more affordable than expected using only the overall and main effects. At the same time, these brands are considered to have higher expected operating costs and smaller cars. The direction of Attractiveness, Safety and Reliability indicate that the combination of these aspects are considered to correlate with higher prices.
Clearly, there are strong similarities but also interesting differences in how these groups of individuals interpreted the performance of the brands on the respective items after adjusting for overall and main effects. In the context of this application, these insights can provide valuable input to a brand positioning strategy, for example. Code reproducing our analysis of these data appears in the appendix of the supplemental materials.

Conclusions
This paper introduces lsbclust, a modelling framework for three-way data, where one of the three ways is clustered over whilst the corresponding matrix slices are approximated by low-rank decompositions. The clustering is done simultaneously with respect to up to four different aspects of these matrix slices, namely the overall mean responses, the row means, the column means, and the row-column interactions. These are the four elements of the biadditive (or bilinear) model used to approximate each of the matrix slices. Which of these terms are included in the model depends on the choice of identifiability constraints, as parametrized by δ. We show that in eight out of nine unique choices for δ, the combination of the bilinear model and least-squares loss allows the four clustering problems to be treated independently. This important property greatly simplifies the complexity of the clustering problem, which also has positive implications for model selection and the interpretation of the results. The low-rank decompositions of the interaction cluster means lead to readily interpretable biplots which aid in the interpretation of the results. As illustrated in an application, lsbclust is a useful and natural alternative to more traditional three-way matrix decomposition methods such as parafac/candecomp and tuckals3. Our method uses a combination of well-known multivariate statistical methods, namely k-means cluster analysis, low-rank decompositions of two-way matrices as well as biplots, whereas traditional three-way methods require domainspecific expertise. Since least-squares loss functions are used, the problems can be treated very efficiently in software. Such software implementing lsbclust has been developed in the form of an eponymous R (R Core Team 2020) package. The package, lsbclust (Schoonees 2019), is available for download from the Comprehensive R Archive Network (CRAN, http://cran.r-project.org).
There are some points that require further research. The treatment of missing values have not been discussed, and should be investigated in the future. In terms of model selection, a wide variety of alternatives to the scree test can and should be investigated. There are a number of promising methods available in the literature, including using multiple criteria and taking a vote to determine the most attractive choice. We note that the rank of the low-rank decomposition can also be considered as a model selection step. Furthermore, it would be possible to add case weights to the methodology. An advantage of case weights is that it allows a mechanism for implementing the bootstrap (e.g. Efron and Tibshirani 1994) to assess the variability of any given solution.
Author Contributions All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Pieter Schoonees, Patrick Groenen and Michel van de Velden. The first draft of the manuscript was written by Pieter Schoonees and Patrick Groenen, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding None.
Availability of data and material All data required to reproduce the results in the manuscript are available in the lsbclust (Schoonees 2019) package for R (R Core Team 2020), which is freely available on the Comprehensive R Archive Network (CRAN).

Conflicts of interest None.
Code availability The reported analyses were conducted using the lsbclust package (Schoonees 2019). See Sect. D in the Supplemental Material for more information.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix:
The appendix contains additional materials cited in the main paper. The first section discusses orthogonality, as referenced in Sect. 4. The second section discusses goodness-of-fit indices for biplots, as mentioned in Sect. 6. Additional information on the lc-bml model is given in the third section. The final section shows how to reproduce the analysis of Sect. 8 using the dedicated R package. (PDF document)

R package lsbclust:
The R package lsbclust contains the code to perform the methods described in the article. The package also contains the data set used in the article. The most up-to-date version can be found on CRAN at http://cran.r-project.org. (GNU zipped tar file)

A Orthogonality
Here we discuss the orthogonality of the decomposition (13), using the notation of Sect. 6.1. Equation (13) represents a special case of the more general notation. For the decomposition (13) to be orthogonal, it must be shown that all six cross-products occurring among the terms in a generalized version of (10) are zero. We treat each of these cross-products in turn.
1. For the cross-product between the interaction term and the row term, it holds that tr J The last equality follows since when δ 2 = 1, J (δ 2 ) When δ 2 = 0, the equality is trivial. 2. For the cross-product between the interaction term and the column term, the result is analogous to the above, except that now the equality 1 J J J = 0 is used. 3. For the cross-product between the interaction term and the term for the overall mean, we have that This cross-product equals zero whenever at least one of the following is true: δ 1 = 1, δ 2 = 1 or δ * = 0. But whenever both δ 1 = δ 2 = 0, δ * = δ 1 δ 3 + δ 2 δ 4 − δ 1 δ 2 = 0 irrespective of δ 3 and δ 4 . Hence the cross-product always equals zero. 4. For the cross-product between the row and column terms, we have Deducing when the cross-product equals zero uses the same concepts as above, but when δ = (1, 1, 0, 0) none of these apply and the cross-product is not necessarily equal to zero.

B Fit diagnostics
The quality of the solution can be assessed using a variety of metrics. We focus on measures for investigating the quality of the interaction approximations. The interaction approximations-case (I), (II), or (III)-allow biplots to be used for visualizing the relationships between the J rating categories and K items for each of the clusters. Biplots generalize scatterplots of two variables to multiple variables (Gower and Hand 1996;Gower et al. 2011), and rely on low-rank inner product approximations. These are most useful when the number of dimensions is low: for example, P ∈ {1, 2, 3}. Constructing biplots for the interactions simply entails plotting the approximation of X u for each cluster. For example, under (III) the rating categories are represented in P dimensions by the rows in C u , while the items are represented in the same space by the rows of D u . The inner products between the pairs of rows in these matrices are rank-P approximations of the corresponding entries in X u . We defer discussion of the interpretation of these biplots to Sect. 8, where empirical examples are examined.
Goodness-of-fit measures for the biplots-and hence, the interactions-are based on the proportion of variation accounted for by the model. A fit value of one indicates perfect fit, while low fit values imply that a substantial amount of variation occurs in the subspace orthogonal to that identified by the model. These are calculated separately for each interaction cluster, as appropriate. Again, we must distinguish between the three cases (I), (II) and (III). We discuss first case (I) here. Cases (II) and (III) are briefly treated after. We do not explicitly include any eventual centering matrices to keep the notation simple.
Measures can be defined for (a) the overall fit, (b) the J rating categories in the rows of X i , as well as for (c) the K items in the columns of X i . For case (I), these are: (a) The overall quality of fit for the interactions within cluster u for P dimensions is This is just the proportion of the variation in the cluster mean explained by the model. Here tr A denotes the trace of the matrix A, which is just the sum of its diagonal elements. (b) The proportion of the variation explained by each of the rows (rating categories), also known as sample predictivities (Gower et al. 2011), is calculated as with each element bounded on [0, 1]. In this context, diag A denotes the diagonal matrix constructed from the main diagonal of A. (c) The column fit can be defined analogously for each of the K items as These quantities are also known as axis predictivities (Gower et al. 2011). Next, we briefly consider goodness-of-fit measures for each X i . The loss contribution for person i towards the interactions is defined as This gives an indication of badness-of-fit, and the sum over all persons gives the minimized value of L (i) G (i) , C, D . These loss contributions account for possible differences in origin, scale and/or rotation between a person's interactions and the modelled cluster mean C u D u . A more informative manner of presenting these loss contributions may be as percentage contributions to L (i) G (i) , C, D . An alternative measure of person fit which is bounded on [−1, 1] is given by This only takes into account differences in rotation and origin, and high values indicate good fit whilst negative values indicate poor fit. When the origins coincide, the quantity (32) can be interpreted as a product-moment correlation coefficient between Vec X i and Vec C u D u . The notation Vec A denotes the vector formed by concatenating the columns of a matrix A into a single vector. Finally, we briefly note fit diagnostics for cases (II) and (III). For case (II), we again have an orthogonal decomposition, namely The row fit can therefore be defined as For case (III), a similar decomposition is available, and we have for the rows, and similarly for the columns.

C The latent-class bilinear multinomial logit model
In this appendix, we briefly expand on the lc-bml model of van Rosmalen et al. (2010) mentioned in Sect. 1.1, which was developed to analyze survey data contaminated with response styles.
Response styles occur when respondents use rating scales heterogeneously (e.g., Schoonees et al. 2015;Baumgartner and Steenkamp 2001). The lc-bml model is a parametric finite mixture of multinomial logit models which models the responses to all items jointly. It simultaneously segments respondents into two types of clusters, namely response style and substantive item segments. Similarly to lsbclust, the lc-bml model produces biplots describing the relationship between the values and the rating categories within each item segment. The response styles are modelled as marginal effects for the rating categories. A nonparametric equivalent of the lc-bml model can be formulated within the lsbclust framework. The resulting model is, except for the inclusion of demographic variables in the lc-bml model, equivalent to the lc-bml model. It has the distinct advantage of being much faster to compute, as least-squares estimation and crisp clustering are used instead of maximum likelihood and finite mixture models.
Whereas lsbclust models entries in X directly, the lc-bml model focuses on modelling the probability of a certain response pattern across a number of items measured on a common rating scale. Mathematically, we have P(X i jk = 1) = S s=1 U u=1 π su P(X i jk = 1 | s, u), where {X i jk } are the random variables whose realizations are the entries in X. The mixing proportions, or a priori class membership probabilities, are denoted by π su , where s and u indexes the two types of latent classes analogous to the two types of clusters in lsbclust. The cluster-specific probabilities are modelled using multinomial logits such that where η jk | s,u is a segment-specific linear predictor. The basic form of the linear predictor is η jk | s,u = α j | s + γ jk | u .
Here α j | s is the attractiveness of rating category j under response style s, and γ jk | u captures the joint effect of rating category j for item k under interaction cluster u, after adjusting for α j | s . To reduce the large number of parameters and to facilitate the use of biplots for interpreting the results, Van Rosmalen et al. (2010) further restrict (33) to be of the form Here c 1 j and d k | u form the rows of C 1 and D u respectively. Identifiability restrictions are applicable to these parameters-see the original paper for more information.
The lc-bml model specification is completed by the likelihood function: This likelihood function is optimized using an EM-algorithm (Dempster et al. 1977), which requires significantly more computation time than our algorithm in Sect. 5. A nonparametric equivalent of the lc-bml model can be formulated within the lsbclust framework. The data array is constructed by transforming each observation into an indicator matrix, with the rows representing the respective rating categories and the columns the survey items. Each column contains a single one indicating which rating was used to answer that item. In effect, we therefore consider the rating scale as one of the modes in our three-way data set. Choosing Model 2 in Table 2 fits a model containing only row, or, in this context, response style effects, and interactions. Additionally, we use the case (II) from Sect. 6.2 (C 1 D u ) as a model for the interactions so that the coordinates for the rating categories are fixed across biplots. The resulting model is, except for the inclusion of demographic variables in the lc-bml model, equivalent to the lc-bml model. The similarity between lc-bml and lsbclust is apparent from comparing (34) and (7). As for lsbclust under case (I), the matrix C 1 contains the coordinates of the rating category effects across all interaction segments in a P-dimensional space, while D u contains the item coordinates for interaction cluster u in that same space.
Note that Van Rosmalen et al. (2010) also include effects for demographic variables in the linear predictor. This is not currently possible for lsbclust.
An empirical comparison of lc-bml and lsbclust, based on an analysis of the list-of-values data conducted in van Rosmalen et al. (2010), is available from the authors upon request.

D Reproducing the empirical example
All computations in this paper were carried out with the lsbclust package (Schoonees 2019) for the open-source statistical software environment R (R Core Team 2020). Version 1.1 of lsbclust is available for download from the Comprehensive R Archive Network (CRAN) and can be installed from within an R session with the command The cars data used in Section 4.1 is available as part of lsbclust as the dcars data set. This is a three-way array containing the data for 187 respondents. We can load the package and the data, and inspect the dimensions of dcars with the commands R> library("lsbclust") R> data("dcars") R> dim(dcars) [1] 10 8 187 The first observation is The object fit can now be queried and plotted to produce the results in the accompanying paper. To create the figures, which are returned as graphical objects created by the ggplot2 package (Wickham 2009), the following code can be used: • A plot containing the information in Table 4: R> plot(fit, type = "overall") The information in Tables 5and 6 can be accessed by the summary() method for the interactions as

R> summary(fit$interactions)
Note that the rows in Table 5 have been reordered such that the fit statistics are in descending order. The three persons in Segment O5 can be identified by looking at the cluster membership vector for the overall means. This can be accessed as part of the "overall" slot of our fitted object