Maximum Likelihood Estimation in Latent Class Models For Contingency Table Data

Statistical models with latent structure have a history going back to the 1950s and have seen widespread use in the social sciences and, more recently, in computational biology and in machine learning. Here we study the basic latent class model proposed originally by the sociologist Paul F. Lazarfeld for categorical variables, and we explain its geometric structure. We draw parallels between the statistical and geometric properties of latent class models and we illustrate geometrically the causes of many problems associated with maximum likelihood estimation and related statistical inference. In particular, we focus on issues of non-identifiability and determination of the model dimension, of maximization of the likelihood function and on the effect of symmetric data. We illustrate these phenomena with a variety of synthetic and real-life tables, of different dimension and complexity. Much of the motivation for this work stems from the"100 Swiss Francs"problem, which we introduce and describe in detail.


Introduction
Latent class (LC) or latent structure analysis models were introduced in the 1950s in the social science literature to model the distribution of dichotomous attributes based on a survey sample from a populations of individuals organized into distinct homogeneous classes on the basis of an unobservable attitudinal feature. See Anderson (1954), Gibson (1955), Madansky (1960) and, in particular, Henry and Lazarfeld (1968). These models were later generalized in Goodman (1974), Haberman (1974), Clogg and Goodman (1984) as models for the joint marginal distribution of a set of manifest categorical variables, assumed to be conditionally independent given an unobservable or latent categorical variable, building upon the then recently developed literature on log-linear models for contingency tables. More recently, latent class models have been described and studied as a special cases of a larger class of directed acyclic graphical models with hidden nodes, sometimes referred to as Bayes nets, Bayesian networks, or causal models, e.g., see Lauritzen (1996), Cowell et al. (1999), Humphreys and Titterington (2003) and, in particular, Geiger et al. (2001). A number of recent papers have established fundamental connections between the statistical properties of latent class models and their algebraic and geometric features, e.g., see Smith (1998, 2005), Smith and Croft (2003), Rusakov and Geiger (2005), Watanabe (2001) and Garcia et al. (2005).
Despite these recent important theoretical advances, the basic statistical tasks of estimation, hypothesis testing and model selection remain surprisingly difficult and, in some cases, infeasible tasks, even for small latent class models. Nonetheless, LC models are widely used and there is a "folklore" associated with estimation in various computer packages implementing algorithms such as EM for estimation purposes, e.g., see Uebersax (2006a,b).
The goal of this article is two-fold. First, we offer a simplified geometric and algebraic description of LC models and draw parallels between their statistical and geometric properties. The geometric framework enjoys notable advantages over the traditional statistical representation and, in particular, offers natural ways of representing singularities and non-identifiability problems. Furthermore, we argue that the many statistical issues encountered in fitting and interpreting LC models are a reflection of complex geometric attributes of the associated set of probability distributions. Second, we illustrate with examples, most of which quite small and seemingly trivial, some of the computational, statistical and geometric challenges that LC models pose. In particular, we focus on issues of non-identifiability and determination of the model dimension, of maximization of the likelihood function and on the effect of symmetric data. We also show how to use symbolic software from computational algebra to obtain a more convenient and simpler parametrization and for unravelling the geometric features of LC models. These strategies and methods should carry over to more complex latent structure models, such as in Bandeen-Roche et al. (1997).
In the next section, we describe the basic latent class model and introduce its statistical properties and issues, and we follow that, in Section 3, with a discussion of the geometry of the models. In Section 4, we turn to our examples exemplifying identifiability issues and the complexity of the likelihood function, with a novel focus on the problems arising from symmetries in the data. Finally, we present some computational results for two real-life examples, of small and very large dimension, and remark on the occurrence of singularities in the observed Fisher information matrix.

Latent Class Models for Contingency Tables
Consider k categorical variables, X 1 , . . . , X k , where each X i takes value on the finite set , R D is the vector space of of k-dimensional arrays of the format d 1 × . . . × d k , with a total of d = i d i entries. The cross-classification of N independent and identically distributed realizations of (X 1 , . . . , X k ) produces a random integer-valued vector n ∈ R D , whose coordinate entry n i i ,...,i k corresponds to the number of times the label combination (i 1 , . . . , i k ) was observed in the sample, for each (i 1 , . . . , i k ) ∈ D. The table n has a Multinomial d (N, p) distribution, where p is a point in the (d − 1)-dimensional probability simplex ∆ d−1 with coordinates p i 1 ,...,i k = P r {(X 1 , . . . , X k ) = (i 1 , . . . , i k )} , (i 1 , . . . , i k ) ∈ D.
Let H be an unobservable latent variable, defined on the set [r] = {1, . . . , r}. In its most basic version, also known as the naive Bayes model, the LC model postulates that, conditional on H, the variables X 1 , . . . , X k are mutually independent. Specifically, the joint distributions of X 1 , . . . , X k and H form the subset V of the probability simplex ∆ dr−1 consisting of points with coordinates where λ h is the marginal probability P r{H = h} and p (h) l (i l ) is the conditional marginal probability P r{X l = i l |H = h}, which we assume to be strictly positive for each h ∈ [r] and (i 1 , . . . , i k ) ∈ D.
The log-linear model specified by the polynomial mapping (1) is a decomposable graphical model (see, e.g, Lauritzen, 1996) and V is the image set of a homeomorphism from the parameter space Θ ≡ θ : θ = (p (h) 1 (i 1 ) . . . p (h) k (i k ), λ h ), (i 1 , . . . , i k , h) ∈ D × [r] = i ∆ d i −1 × ∆ r−1 , so that global identifiability is guaranteed. The remarkable statistical properties of this type of model and the geometric features of the set V are well understood. Statistically, equation (1) defines a linear exponential family of distributions, though not in its natural parametrization. The maximum likelihood estimates, or MLEs, of λ h and p (h) l (i l ) exist if and only if the minimal sufficient statistics, i.e., the empirical joint distributions of (X i , H) for i = 1, 2, . . . , k, are strictly positive and are given in closed form as rational functions of the observed two-way marginal distributions between X i and H for i = 1, 2, . . . , k. The log-likelihood function is strictly concave and the global maximum is always attainable, possibly on the boundary of the parameter space. Furthermore, the asymptotic theory of goodness-of-fit testing is fully developed. The statistical problem arises because H is latent and unobservable.
Geometrically, we can obtain the set V as the intersection of ∆ dr−1 with an affine variety (see, e.g., Cox et al., 1996) consisting of the solutions set of a system of r i d i 2 homogeneous squarefree polynomials. For example, when k = 2, each of these polynomials take the form of quadric equations of the type with i 1 = i 1 , i 2 = i 2 and for each fixed h. Equations of the form (2) are nothing more than conditional odds ratio of 1 for every pair (X i , X i ) given H = hand, for each given h, the coordinate projections of the first two coordinates of the points satisfying (2) trace the surface of independence inside the simplex ∆ d−1 . The strictly positive points in V form a smooth manifold whose dimension is r i (d i − 1) + (r − 1) and whose co-dimension corresponds to the number of degrees of freedom. The singular points in V all lie on the boundary of the simplex ∆ dr−1 and identify distributions with degenerate probabilities along some coordinates. The singular locus of V can be described similarly in terms of stratified components of V, whose dimensions and co-dimensions can also be computed explicitly. Under the LC model, the variable H is unobservable and the new model H is a r-class mixture over the exponential family of distributions prescribing mutual independence among the manifest variables X 1 , . . . , X k . Geometrically, H is the set of probability vectors in ∆ d−1 obtained as the image of the marginalization map from ∆ dr−1 onto ∆ d−1 which consists of taking the sum over the coordinate corresponding to the latent variable. Formally, H is made up of of all probability vectors in ∆ d−1 with coordinates satisfying the accounting equations (see, e.g., Henry and Lazarfeld, 1968) Despite being expressible as a convex combination of very well-behaved models, even the simplest form of the LC model (3) is far from well-behaved and, in fact, shares virtually none of the properties of the standard log-linear models (1) described above. In particular, latent class models described by equations (3) do not define exponential families, but instead belong to a broader class of models called stratified exponential families (see Geiger et al., 2001), whose properties are much weaker and less well understood. The minimal sufficient statistics for an observed table n are the observed counts themselves and we can achieve no data reduction via sufficiency. The model may not be identifiable, because for a given p ∈ ∆ d−1 defined by (3), there may be a subset of Θ, known as the non-identifiable space, consisting of parameter points all satisfying the same accounting equations. The non-identifiability issue has in turn considerable repercussions for the determination of the correct number of degrees of freedom for assessing model fit and, more importantly, on the asymptotic properties of standard model selection criteria (e.g. likelihood ratio statistic and other goodness-of-fit criteria such as BIC, AIC, etc), whose applicability and correctness may no longer hold.
Computationally, maximizing the log-likelihood can be a rather laborious and difficult task, particularly for high dimensional tables, due to lack of concavity, the presence of local maxima and saddle points, and singularities in the observed Fisher information matrix. Geometrically, H is no longer a smooth manifold on the relative interior of ∆ d−1 , with singularities even at probability vectors with strictly positive coordinates, as we show in the next section. The problem of characterizing the singular locus of H and of computing the dimensions of its stratified components (and of the tangent spaces and tangent cones of its singular points) is of statistical importance: singularity points of H are probability distributions of lower complexity, in the sense that they are specified by lower-dimensional subsets of Θ, or, loosely speaking, by less parameters. Because the sample space is discrete, although the singular locus of H has typically Lebesgue measure zero, there is nonetheless a positive probability that the maximum likelihood estimates end up being either a singular point in the relative interior of the simplex ∆ d−1 or a point on the boundary. In both cases, standard asymptotics for hypothesis testing and model selection fall short.

Geometric Description of Latent Class Models
In this section, we give a geometric representation of latent class models, summarize existing results and point to some of the relevant mathematical literature. For more details, see Garcia et al. (2005) and Garcia (2004).
The latent class model defined by (3) can be described as the set of all convex combinations of all r-tuple of points lying on the surface of independence inside ∆ d−1 . Formally, let be the map that sends the vectors of marginal probabilities into the k-dimensional array of joint probabilities for the model of complete independence. The set S ≡ σ(∆ d 1 −1 × . . . × ∆ d k −1 ) is a manifold in ∆ d−1 known in statistics as the surface of independence and in algebraic geometry (see, e.g. Harris, 1992) as (the intersection of ∆ d−1 with) the Segre embedding of , the dimension of the corresponding decomposable model of mutual independence. The set H can then be constructed geometrically as follows. Pick any combination of r points along the hyper-surface S, say p (1) , . . . , p (r) , and determine their convex hull, i.e. the convex subset of ∆ d−1 consisting of all points of the form h p (h) λ h , for some choice of (λ 1 , . . . , λ r ) ∈ ∆ r−1 . The coordinates of any point in this new subset satisfy, by construction, the accounting equations (3). In fact, the closure of the union of all such convex hulls is precisely the latent class model H. In algebraic geometry, H would be described as the intersection of ∆ d−1 with the r-th secant variety of the Segre embedding mentioned above. Example 3.1. The simplest example of a latent class model is for a 2 × 2 table with r = 2 latent classes. The surface of independence, i.e. the intersection of the simplex ∆ 3 with the Segre variety, is shown in Figure 1. The secant variety for this latent class models is the union of all the secant lines, i.e. the lines connecting any two distinct points lying on the surface of independence. Figure  1 displays three such secant lines. It is not to hard to picture that the union of all such secant lines is the enveloping simplex ∆ 3 and, therefore, H fills up all the available space (for formal arguments, see Catalisano et al., 2002, Proposition 2.3).
The model H, thought of as a portion of the r-th secant variety to the Segre embedding, is not a smooth manifold. Instead, it is a semi-algebraic set (see, e.g., Benedetti, 1990), clearly singular on the boundary of the simplex, but also at strictly positive points along the (r − 1)st secant variety, (both of Lebesgue measure zero). This means that the model is singular at all points in H which satisfy the accounting equations with one or more of the λ h 's equal to zero. In Example 3.1 above, the surface of independence is a singular locus for the latent class model. From the statistical viewpoint, singular points of H correspond to simpler models for which the number of latent classes is less than r (possibly 0). As usual, for these points one needs to adjust the number of degrees of freedom to account for the larger tangent space.
Unfortunately, we have no general closed-form expression for computing the dimension of H and the existing results only deal with specific cases. Simple considerations allow us to compute an upper bound for the dimension of H, as follows. As Example 3.1 shows, there may be instances for which H fills up the entire simplex ∆ d−1 , so that d − 1 is an attainable upper bound. Counting the number of free parameters in (3), we can see that this dimension cannot exceed r i (d i −1)+r −1, c.f., (Goodman, 1974, page 219). This number, the standard dimension, is the dimension of the fully observable model of conditional independence. Incidentally, this value can be determined mirroring the geometric construction of H as follows (c.f., Garcia (2004)). The number r i (d i − 1) arises from the choice of r points along the i (d i − 1)-dimensional surface of independence, while the term r − 1 accounts for the number of free parameters for a generic choice of (λ 1 , . . . , λ r ) ∈ ∆ r−1 . Therefore, we conclude that the dimension of H is bounded by a value known in algebraic geometry as the expected dimension the variety H. Cases of latent class models with dimension strictly smaller than the expected dimension have been known for a long time, however. In the statistical literature, Goodman (1974) noticed that the latent class models for 4 binary observable variables and a 3-level latent variable, whose expected dimension is 14, has dimension 13. In algebraic geometry, secant varieties with dimension smaller than the expected dimension (4) are called deficient (e.g., see Harris, 1992). In particular, Exercise 11.26 in Harris (1992) gives an example of deficient secant variety, which corresponds to a latent class model for a 2-way table with a latent variable taking on 2 values. In this case, the deficiency is 2, as is demonstrated below in equation (5). The true or effective dimension of a latent class model, i.e. the dimension of the semi-algeraic set H representing it, is crucial for establishing identifiability and for computing correctly the number of degrees of freedom. In fact, if a model is deficient, then that the pre-image of each probability array in H arising from the accounting equations is a subset (in fact, a variety) of Θ called the non-dentifiable subspace, with dimension exactly equal to the deficiency itself. Therefore, a deficient model is non-identifiable, with adjusted degrees of freedom equal to number of degrees of freedom for the observable graphical model plus the value of the deficiency.
Theoretically, it is possible to determine the effective dimension of H by computing the maximal rank of the Jacobian matrix for the polynomial mapping from Θ into H given coordinatewise by (3). In fact, Geiger et al. (2001) showed that this value is equal to the dimension of H almost everywhere with respect to the Lebsegue measure, provided the Jacobian is evaluated at strictly positive parameter points θ. These symbolic evaluations, however, require the use of symbolic software which can only handle small tables and models, so that, in practice, computing the effective dimension of a latent class model is computationally difficult and often unfeasible.
Recently, in the algebraic-geometry literature, Catalisano et al. (2002Catalisano et al. ( , 2003 have obtained explicit formulas for the effective dimensions of some secant varieties which are of statistical interest. In particular, they show that for k = 3 and r ≤ min{d 1 , d 2 , d 3 }, the latent class model has the expected dimension and is identifiable. On the other hand, assuming Finally, under the same conditions, H is identifiable when 1 2 i (d i − 1) + 1 ≥ max{d k , r}. Obtaining bounds and results of this type is highly non-trivial and is an open area of research.
In the remainder of the paper, we will focus on simpler latent class models for tables of dimension k = 2 and illustrate with examples the results mentioned above. For latent class models on two-way tables, there is an alternative, quite convenient way of describing H by representing each p in ∆ d−1 as a d 1 × d 2 matrix and by interpreting the map σ as a vector product. In fact, each point p in S is a rank one matrix obtained as p 1 p 2 , where p 1 ∈ ∆ d 1 −1 and p 2 ∈ ∆ d 1 −2 are the appropriate marginal distributions of X 1 and X 2 . Then, the accounting equations for a latent class models with r-level become i.e. the matrix p is a convex combination of r rank 1 matrices lying on the surface of independence. Therefore all points in H are non-negative matrices with entries summing to one and with rank at most r. This simple observation allows one to compute the effective dimension of H for 2way table as follows. In general, a real valued d 1 × d 2 matrix has rank r or less if and only if the homogeneous polynomial equations corresponding to all of its (r + 1) × (r + 1) minors all vanish. Provided k < min{d 1 , d 2 }, on R d 1 × R d 2 , the zero locus of all such such equations form a determinantal variety of co-dimension (d 1 − r)(d 2 − r) (Harris, 1992, Proposition 12.2) and hence has dimension r(d 1 +d 2 )−r 2 . Subtracting this value from the expected dimension computed above, and taking into account the fact that all the points lie inside the simplex, we obtain This number is also the difference between the dimension of the (fully identifiable, i.e. of expected dimension) graphical model of conditional independence X 1 and X 2 given H, and the deficient dimension of the latent class model obtained by marginalizing over the variable H.
The study of higher dimensional tables is still an open area of research. The mathematical machinery required to handle larger dimensions is considerably more complicated and relies on the notions higher-dimensional tensors, rank tensors and non-negative rank tensors, for which only partial results exist. See Kruskal (1975), Cohen and Rothblum (1993) and Strassen (1983) for details. Alternatively, Mond et al. (2003) conduct an algebraic-topological investigation of the topological properties of stochastic factorization of stochastic matrices representing models of conditional independence with one hidden variable and Rhodes (2006, 2007) explore an overlapping set of problems framed in the context of trees with latent nodes and branches.
The specific case of k-way tables with 2 level latent variables is a fortunate exception, for which the results for 2-way tables just described apply. In fact, Landsberg and Manivel (2004) show that that these models are the same as the corresponding model for any two-dimensional table obtained by any "flattening" of the d 1 × . . . × d k -dimensional array of probabilities p into a two-dimensional matrix. Flattening simply means collapsing the k variables into two new variables with f 1 and f 2 levels, and re-organizing the entries of the k-dimensional tensor p ∈ ∆ d−1 into a f 1 × f 1 matrix accordingly, where, necessarily, f 1 + f 2 = i d i . Then, H is the determinantal variety which is the zero set of all 3 × 3 sub-determinants of the matrix obtained by any such flattening. The second example in Section 4.1 below illustrates this result.

Examples Involving Synthetic Data
We further elucidate the non-identifiability phenomenon from the algebraic and geometric point of view, and the multi-modality of the log-likelihood function issue using few, small synthetic examples. In particular, in the "100 Swiss Frank" problem below, we embark on a exhaustive study of a table with symmetric data and describe the effects of such symmetries on both the parameter space and the log-likelihood function. Although this example involves one of the simplest cases of LC models, these tables already exhibit considerable statistical and geometric complexity.

Effective Dimension and Polynomials
We show how it is possible to take advantage of the polynomial nature of equations (3) to gain further insights into the algebraic properties of distributions obeying latent class models. All the computations that follow were made in SINGULAR (Greuel et al., 2005) and are described in details, along with more examples, in Zhou (2007). Although in principle symbolic algebraic software allows one to compute the set of polynomial equations that fully characterize LC models and their properties, this is still a rather difficult and costly task that can be accomplished only for smaller models.
The accounting equations (3) determine a polynomial mapping f : Θ → ∆ d−1 given by so that the latent class model can be analytically defined as the image of this map, i.e. H = f (Θ). Then, following the geometry-algebra dictionary principle (see, e.g., Cox et al., 1996), the problem of computing the effective dimension of H can in turn be geometrically cast as a problem of computing the dimension of the image of a polynomial map. We illustrate how this representation offers considerable advantages with some small examples. Consider a 2 × 2 × 2 table with r = 2 latent classes. From Proposition 2.3 in Catalisano et al. (2002), the latent class models with 2 classes and 3 manifest variables are identifiable. The standard dimension, i.e. the dimension of the parameter space Θ is r i (d i − 1) + r − 1 = 7, which coincides with the dimension of the enveloping simplex ∆ 7 . Although this condition implies that the number of parameters to estimate is no larger than the number of cells in the table, a case which, if violated, would entail non-identifiability, it does not guarantee that the effective dimension is also 7. This can be verified by checking that the symbolic rank of the Jacobian matrix of the map (6) is indeed 7, almost everywhere with respect to the Lebesgue measure. Alternatively, one can determine the dimension of the non-identifiable subspace using computational symbolic algebra. First, we define the ideal of polynomials determined by the 8 equations in (6) in the polynomial ring in which the (redundant) 16 indeterminates are the 8 joint probabilities in ∆ 7 and the 3 pairs of marginal probabilities in ∆ 1 for the observable variables, and the marginal probabilities in ∆ 1 for the latent variable. Then we use implicization (Cox et al., 1996, Chapter 3) to eliminate all the marginal probabilities and to study the Groebner basis of the resulting ideal in which the indeterminates are the joint probabilities only. There is only one element in the basis, which gives the trivial condition for probability vectors. This implies the map (6) is surjective, so that H = ∆ 7 and the effective dimension is also 7, showing identifiability, at least for positive distributions.
Next, we consider the 2×2×3 table with r = 2. For this model Θ has dimension 9 and the image of the mappings (6) is ∆ 9 . The symbolic rank of the associated Jacobian matrix is 9 as well and the model is identifiable. The image of the polynomial mapping determined by (6) is the variety associated to the ideal whose Groebner basis consists of the trivial equaiton  (7) where the subscript symbol "+" indicates summation over that coordinate. In turn, the zero set of the above determinants coincide with the determinantal variety specified by the zero set of all 3 × 3 minors of the 3×4 matrix  which is a flattening of the 2 × 2 × 3 array of probabilities describing the joint distribution for the latent class model under study. This is in accordance with the result in Landsberg and Manivel (2004) of mentioned above. Now, the determinantal variety given by the vanishing locus of all the 3 × 3 minors of the matrix (8) is the latent class model for a 3 × 4 table with 2 latent classes, which, according to (5), has deficiency equal to 2. The effective dimension of this variety is 9, computed as the standard dimension, 11, minus the deficiency. Then,the effective dimension of the model we are interested is also 9 and we conclude that the model is identifiable. Table 1 summarizes some of our numerical evaluations of the different notions of dimension for a different LC models. We computed the effective dimensions by evaluating with MATLAB the numerical rank of the Jacobian matrix, based on the simple algorithm suggested in Geiger et al. (2001) and also using SINGULAR, for which only computations involving small models were feasible. Table 1: Different dimensions of some latent class models. The Complete Dimension is the dimension d − 1 of the envoloping probability simplex ∆ d−1 . See also Table 1 in Kocka and Zhang (2002).
For the basic latent class model, the standard dimension of Θ = ∆ 3 × ∆ 3 × ∆ 1 is 2(3 + 3) + 1 = 13 and, by (5), the deficiency is 2. Thus, the model is not identifiable and the pre-image of each point p ∈ H by the map (6) is a 2-dimensional surface in Θ. To keep the notation light, we write α ih for p . . , 4 and α (h) and β (h) for the conditional marginal distribution of X 1 and X 2 given H = h, respectively. The accounting equations for the points in H become and the log-likelihood function, ignoring an irrelevant additive constant, is It is worth emphasizing, as we did above and as the previous display clearly shows, that the observed counts are minimal sufficient statistics. Alternatively, we can re-parametrize the log-likelihood function using directly points in H rather the points in the parameter space Θ. Recall from our discussion in section 3 that, for this model, the 4 × 4 array p is in H if and only if each 3 × 3 minor vanishes. Then, we can write the log-likelihood function as where p * ij is the 3 × 3 sub-matrix of p obtained by erasing the ith row and the jth column. Although the first order optimality conditions for the Lagrangian corresponding to the parametrization (11) are algebraically simpler and can be given the form of a system of a polynomial equations, in practice, the classical parametrization (10) is used in both the EM and the Newton-Raphson implementations in order to compute the maximum likelihood estimate of p. See Goodman (1979), Haberman (1988), and Redner and Walker (1984) for more details about these numerical procedures.

Global and Local Maxima
Using both EM and Newton-Raphson algorithm with several different starting points, we found 7 local maxima of the log-likelihood function, reported in Table 2. The global maximum was found experimentally to be −20.8074 + const., where const. denotes the additive constant stemming from the multinomial coefficient. The maximum is achieved by the three tables of fitted values Table 2 a). The remaining four tables are local maximum of −20.8616 + const., close in value to the actual global maximum. Using SINGULAR (see (Greuel et al., 2005)), we checked that the tables found satisfy the first order optimality conditions (11). After verifying numerically the second order optimality conditions, we conclude that those points are indeed local maxima. Furthermore, as indicated in Pachter and Sturmfels (2005), the log-likelihood function also has a few saddle points.
A striking feature of the global maxima in Table 2 is their invariance under the action of the symmetric group on four elements acting simultaneously on the row and columns. Different symmetries arise for the local maxima. We will give an explicit representation of these symmetries under the classical parametrization (10) in the next section.
Despite the simplicity and low-dimensionality of the LC model for the Swiss franc problem and the strong symmetric features of the data, we have yet to provide a purely mathematical proof that the three top arrays in Table 2 correspond to a global maximum of the likelihood function. We view the difficulty and complexity of the 100 Swiss Francs problem as a consequence of the inherent difficulty of even small LC models and perhaps an indication that the current theory has still many open, unanswered problems. In Section 6, we present partial results towards the completion of the proof.

Unidentifiable Space
It follows from equation (5) that the non-identifiable subspace is a two-dimensional subset of Θ. We give an explicit algebraic description of this space, which we will then use to obtain interpretable plots of the profile likelihood.
Firstly, we focus on the three global maxima in Table 2 a). By the well-known properties of the EM algorithm (see, e.g., Pachter and Sturmfels, 2005, Theorem 1.15), if the vector of parameters θ is a stationary point in the maximization step of the EM algorithm, then θ is a critical point and  (13), when evaluated over the ball in R 3 of radius 3, centered at the origin. The inner box is the unit cube [0, 1] 3 . hence a good candidate for a local maximum. Using this observation, it is possible to show (see Zhou, 2007) that any point in Θ satisfying the equations is a stationary point. Notice that the first four equations in (12) require α (h) and β (h) to each have the first and second pairs of coordinates identical, for h = 1, 2. The equation (12) defines a 2-dimensional surface in Θ. Using SINGULAR, we can verify that, holding, for example, α 11 and β 11 fixed, determines all of the other parameters according to the equations Using elimination (see Cox et al., 1996, Chapter 3) to remove all the variables in the system except for λ 1 , we are left with one equation 80λ 1 α 11 β 11 − 20λ 1 α 11 − 20λ 1 β 11 + 6λ 1 − 1 = 0.
Without the constraints for the coordinates of α 11 , β 11 and λ 1 to be probabilities, (13) defines a twodimensional object in R 3 , depicted in Figure 2. Notice that the axes do not intersect this surface, so that zero is not a possible value for α 11 , β 11 and λ 1 . Because the non-identifiable space in Θ is 2-dimensional, equation (13) actually defines a bijection between α 11 , β 11 and λ 1 and the rest of the parameters. Then, the intersection of the surface (13) with the unit cube [0, 1] 3 , depicted as a red box in Figure 2, is the projection of the whole non-identifiable subspace into the 3-dimensional unit cube where α 11 , β 11 and λ 1 live. Figure 3 displays two different views of this projection.
The preceding arguments hold unchanged if we replace the symmetry conditions in the first two lines of equation (12) with either of these other two conditions, requiring different pairs of coordinates to be identical, namely and where h = 1, 2. By our computations, the non-identifiable surfaces inside Θ corresponding each to one of the three pairs of coordinates held fixed in equations (12), (14) and (15), produce the three distinct tables of maximum likelihood estimates reported in Table 2 a). Figure 3 shows the projection of the non-identifiable subspaces for the three MLEs in Table 2 a) into the three dimensional unit cube for λ 1 , α 11 and β 11 . Although each of these three subspaces are disjoint subsets of Θ, their lower dimensional projections comes out as unique. By projecting onto the different coordinates λ 1 , α 11 and β 21 instead, we obtain two disjoint surfaces for the first, and second and third MLE, shown in Figure 4. Table 3 presents some estimated parameters using the EM algorithm. Though these estimates are hardly meaningful, because of the non-identifiability issue, they show the symmetry properties we pointed out above and implicit in equations (12), (14) and (15), and they explain the invariance under simultaneous permutation of the fitted tables. In fact, the number of global maxima is the number of different configurations of the 4 dimensional vectors of estimated marginal probabilities with two identical coordinates, namely 3. This phenomenon, entirely due to the strong symmetry in the observed table (9), is completely separate from the non-identrifiability issues, but just as problematic.
By the same token, we can show that vectors of marginal probabilities with 3 identical coordinates also produce stationary points for the EM algorithms. This type of stationary points trace surfaces inside Θ which determine the local maxima of Table 2 b). The number of these local maxima corresponds, in fact, to the number of possible configurations of 4-dimensional vectors with 3 identical coordinates, namely 4. Figure 5 depicts the lower dimensional projections into λ 1 , α 11 and β 11 of the non-identifiable subspaces for the first MLE in Table 2 a), the first three local maxima and the last local maxima in Table 2 b).
We can summarize our finding as follows: the maxima in Table 2 define disjoint 2-dimensional surfaces inside the parameter space Θ, the projection of one of them depicted in Figure 3. While non-identifiability is a structural feature of these models which is independent of the observed data, the multiplicity and invariance properties of the maximum likelihood estimates and the other local maxima is a phenomenon cause by the symmetry in the observed table of counts.

Plotting the Log-likelihood Function
Having determined that the non-identifiable space is 2-dimensional and that there are multiple maxima, we proceed with some plots of the profile log-likelihood function. To obtain a non-trivial surface, we need to consider three parameters. Figures 9 and 7 display the surface and contour plot of the profile log-likelihhod function for α 11 and α 21 when α 31 is one of the fixed parameters. Both Figures show clearly the different maxima of the log-likelihood function, each lying on the top of "ridges" of the log-likelihood surface which are placed symmetrically with respect to each others. The position and shapes of these ridges reflect, once again, the invariance properties of the estimated probabilities and parameters.

Further Remarks and Open Problem
We conclude this section with some observations and pointers to open problems.
One of the interesting aspects we came across while fitting the table (9) was the proximity of the values of the local and global maxima of the log-likelihood function. Furthermore, although these values are very close, the fitted tables corresponding to global and local maxima are remarkably different. Even though the data (9) are not sparse, we wonder about the effect of cell sizes. Figure  8 show the same profile log-likelihood for the table (9) multiplied by 10000. While the number of global and local maxima, the contour plot and the basic symmetric shape of the profile loglikelihood surface remain unchanged after this rescaling, the peaks around the global maxima have become much more pronounced and so has the difference between of the values of the global and local maxima.
We have studied at a number of variations of table (9), focussing in particular on the symmetric data. We report only some of our results and refer to Zhou (2007) for a more extensive study. Table  4 shows the values and number of local and global maxima for a the 6 × 6 version of (9). As for the 4 × 4 case, we notice strong invariance features of the various maxima of the likelihood function and a very small difference between the value of the global and local maxima. Fitting the same model to the table     1 2 2 2 2 1 2 2 2 2 1 2 2 2 2 1     we found 6 global maxima of the likelihood function, which give as many maximum likelihood estimates, all obtainable via simultaneous permutation of rows and columns of the table     7/4 7/4 7/4 7/4 7/4 7/4 7/4 7/4 7/4 7/4 7/6 7/3 7/4 7/4 7/3 7/6     , log-likelihood = −77.2927 + const.
Based on the various cases we have investigated, we have the following conjecture, which we verified computationally up to dimension k = 50:

Example: Michigan Influenza
Despite the reduced dimensionality of this problem and the large sample size, we report on the instability of the Fisher scoring algorithm implemented in the R package gllm, e.g., see Espeland (1986). As the algorithm cycles through, the evaluations of Fisher information matrix become increasing ill-conditioned and eventually produce instabilities in the estimated coefficients and in the standard errors. These problems disappear in the modified Newton-Raphson implementation, originally suggested by Haberman (1988), based on an inexact line search method known in the convex optimization literature as the Wolfe conditions. Erosheva (2002) and Erosheva et al. (2007) analyze an extract from the National Long Term Care Survey in the form of a 2 16 contingency table that contains data on 6 activities of daily living (ADL) and 10 instrumental activities of daily living (IADL) for community-dwelling elderly from 1982, 1984, 1989, and 1994 survey waves. The 6 ADL items include basic activities of hygiene and personal care (eating, getting in/out of bed, getting around inside, dressing, bathing, and getting to the bathroom or using toilet). The 10 IADL items include basic activities necessary to reside in the community (doing heavy housework, doing light housework, doing laundry, cooking, grocery shopping, getting about outside, travelling, managing money, taking medicine, and telephoning). Of the 65,536 cells in the table, 62384 (95.19%) contain zero counts, 1729 (2.64%)contain counts of 1, 499 (0.76%) contain counts of 2. The largest cell count, corresponding to the (1, 1, . . . , 1) cell, is 3853.  Erosheva (2002) and Erosheva et al. (2007) use an individual-level latent mixture model that bears a striking resemblance to the LC model. Here we report on analyses with the latter.

Data From the National Long Term Care Survey
We use both the EM and Newton-Raphson algorithms to fit a number of LC models with up to 20 classes, which can be shown to be all identifiable in virtue of Proposition 2.3 in Catalisano et al. (2002). Table 6 reports the maximal value of log-likelihood function and the value of BIC (the Bayesian Information Criterion), which seem to indicate that larger LC models with many levels are to be preferred. To provide a better sense of how well these LC models fit the data, we show in Table 7 the fitted values for the six largest cells, which, as mentioned, deviates considerably from most of the cell entries. We have also considered alternative model selection criteria such as AIC and modifications of it. AIC (with and without a 2nd order correction) points to k > 20! (An adhoc modification of AIC due to Anderson et al. (1994) for overdispersed data gives rather bizarre results.) The dimensionality of a suitable LC model for these data appears to be much greater than for the individual level mixture model in Erosheva et al. (2007). Because of its high dimensionality and remarkable degree of sparsity, this example offers an ideal setting in which to test the relative strengths and disadvantages of the EM and Newton-Raphson algorithm. In general, the EM algorithm, as a hill-climbing method, moves steadily towards solutions with higher value of the log-likelihood, but converges only linearly. On the other hand, despite its faster quadratic rate of convergence, the Newton-Raphson method tends to be very time and space consuming when the number of variables is large, and may be numerically unstable if the Hessian matrices are poorly conditioned around critical points, which again occurs more frequently in large problems (but also in small ones, such as the Michigan Influenza examples above).
For the class of basic LC models considered in this paper, the time complexity for one single step of the EM algorithm is O (d · r · i d i ), while the space complexity is O (d · r). In contrast, for the Newton-Raphson algorithm, both the time and space complexity are O d · r 2 · i d i . Consequently, for the NLTCS dataset, when r is bigger than 4, Newton-Raphson is sensibly slower than EM, and when r goes up to 7, Newton-Raphson needs more than 1G of memory. Another significant drawback of the Newton-Raphson method we experienced while fitting both the Michigan influenza and the NLTCS datasets is its potential numerical instability, due to the large condition numbers of the Hessian matrices. As remarked at the end of the previous section, following Haberman (1988), a numerically convenient solution is to modify the Hessian matrices so that they remain negative definite and then approximate locally the log-likelihood by a quadratic function. However, since the log-likelihood is neither concave and nor quadratic, these modifications do not necessarily guarantee an increase of the log-likelihood at each iteration step. As a result, the algorithm may experience a considerable slowdown in the rate of convergence, which we in fact observed with the NLTCS data. Table 8 shows the condition numbers for the true Hessian matrices evaluated at the numerical maxima, for various values of r. This table suggests that, despite full identifiability, the log-likelihood has a very low curvature around the maxima and that the log-likelihood may, in fact, look quite flat.
To elucidate this point and some of the many difficulties in fitting LC models, we show in Figure  9 the profile likelihood plot for the parameter α 12 in simplest LC model with r = 2. The actual profile log-likelihood is shown in red and is obtained as the upper envelop of two distinct, smooth curves, each corresponding to a local maxima of the log-likelihood. The location of the optimal value of α 12 is displayed with a vertical line. Besides illustrating multimodality, the log-likelihood function in this example is notable for its relative flatness around its global maximum.

On Symmetric Tables and the MLE
In this section, we show how symmetry in data allows one to symmetrize via averaging local maxima of the likelihood function and to obtain critical points that are more symmetric.
In various examples we looked at, these have larger likelihood than the tables from which they are obtained. We also prove that if the aforementioned averaging process always causes likelihood to go up, then among the 4 × 4 matrices of rank 2, the ones maximizing the log-likelihhod function for the 100 Swiss Francs problem (9) are given in Table 2 a). We will further simplify the notation and will write L for the matrix of observed counts and M for the matrix of MLEs.

Introduction and Motivation
A main theme in this section is to understand in what ways symmetry in data forces symmetry in the global maxima of the likelihood function. One question is whether our ideas can be extended at all to nonsymmetric data by suitable scaling. We prove that nonsymmetric local maxima will imply the existence of more symmetric points which are critical points at least within a key subspace and are related in a very explicit way to the nonsymmetric ones. Thus, if the EM algorithm leads to a local maximum which lacks certain symmetries, then one may deduce that certain other, more symmetric points are also critical points (at least within certain subspaces), and so check these to see if they give larger likelihood. There is numerical evidence that they do, and also a close look at our proofs shows that for "many" data points this symmetrization process is guaranteed to increase maximum likelihood, by virtue of a certain single-variable polynomial encoding of the likelihood function often being real-rooted.
Here is an example of our symmetrization process. Given the data 4 2 2 2 2 2 2 4 2 2 2 2 2 2 4 2 2 2 2 2 2 4 2 2 2 2 2 2 4 2 2 2 2 2 2 4 , one of the critical points located by the EM algorithm is 7/3 7/3 7/3 7/3 7/3 7/3 7/3 13/5 13/5 13/5 29/15 29/15 7/3 13/5 13/5 13/5 29/15 29/15 7/3 13/5 13/ in which two coordinates are averaged; however, since one of the values being averaged is zero, it is not so clear whether this should increase likelihood. However, repeatedly applying such symmetrization steps to this example, does converge to a local maximum. Now let us speak more generally. Let M be an n by n matrix of rank at most two which has row and column sums all equalling kn, implying (by results of Section 6.2) that we may write M i,j as k + e i f j where e, f are each vectors whose coordinates sum to 0. We are interested in the following general question: Question 6.1. Suppose a data matrix is fixed under simultaneously swapping rows and columns i, j.
Consider any M as above, i.e. with M i,j = k + e i f j . products also satisfied. Does e i > e j > 0, f i > f j > 0 (or similarly e i < e j < 0, f i < f j < 0 ) imply that replacing e i , e j each by e i +e j 2 and f i , f j each by

always increases the likelihood?
Remarks The weaker conditions e i > e j = 0 and f i > f j = 0 (resp. e i < e j = 0, f i < f j = 0) do not always imply that this replacement will increase likelihood. However, one may consider the finite list of possibilities for how many zeroes the vectors e and f may each have; an affirmative answer to Question 6.1 would give a way to find the matrix maximizing likelihood in each case, and then we could compare this finite list of maxima to find the global maximum. One may check that the example discussed in Chapter 1 of Pachter and Sturmfels (2005) gives another instance where this averaging approach leads quickly to what appears to be a global maximum. Namely, given the data matrix 4 2 2 2 2 4 2 2 2 2 4 2 2 2 2 4 and a particular starting point, the EM algorithm converges to the saddle point 1 48 4 2 3 3 2 4 3 3 3 3 3 3 3 3 3 3 , whose entries may be written as M i,j = 1/48(3 + a i b j ) for a = (−1, 1, 0, 0) and b = (−1, 1, 0, 0). Averaging −1 with 0 and 1 with the other 0 simultaneously in a and b immediately yields the global maximum directly by symmetrizing the saddle point, i.e. rather than finding it by running the EM algorithm repeatedly from various starting points. An affirmative answer to Question 6.1 would imply several things. It would yield a (positive) solution to the 100 Swiss francs problem, as discussed in Section 6.3. More generally, it would explain in a rather precise way how certain symmetries in data seem to impose symmetry on the global maxima of the maximum likelihood function. Moreover it would suggest good ways to look for global maxima, as well as constraining them enough that in some cases they can be characterized, as we demonstrate for the 100 Swiss francs problem. To make this concrete, one thing it would tell us for an n by n data matrix which is fixed by the S n action simultaneously permuting rows and columns in the same way, is that any probability matrix maximizing likelihood for such a data matrix will have at most two distinct types of rows.
We do not know the answer to this question, but we do prove that this type of averaging will at least give a critical point within the subspace in which e i , e j , f i , f j may vary freely but all other parameters are held fixed. Data also provides evidence that the answer to the question may very well be yes. At the very least, this type of averaging appears to be a good heuristic for seeking local maxima, or at least finding a way to continue to increase maximum likelihood beyond what it is at a critical point one reaches. Moreover, while real data is unlikely to have these symmetries, perhaps it could come close, and this could still be a good heuristic to use in conjunction with the EM algorithm.

Preservation of Marginals and Some Consequences
Proposition 6.3. Given data in which all row and column sums (i.e. marginals) are equal, then for M to maximize the likelihood function for this data among matrices of a fixed rank, row and column sums of M all must be equal.
We prove the case mentioned in the abstract, which should generalize by adjusting exponents and ratios in the proof. It may very well also generalize to distinct marginals and tables with more rows and columns.
Proof. Let R 1 , R 2 , R 3 , R 4 be the row sums of M . Suppose R 1 ≥ R 2 ≥ R 3 > R 4 ; other cases will be similar. Choose δ so that R 3 = (1 + δ)R 4 . We will show that multiplying row 4 by any 1 + with 0 < < min(1/4, δ/2) will strictly increase L, giving a contradiction to M maximizing L. The result for column sums follows by symmetry.
Let us write L(M ) for the new matrix M in terms of the variables x i,j for the original matrix M , so as to show that L(M ) > L(M ). The first inequality below is proven in Lemma 6.4.
Proof. It is equivalent to show Corollary 6.5. There exist vectors (e 1 , e 2 , e 3 , e 4 ) and (f 1 , f 2 , f 3 , f 4 ) such that 4 i=1 e i = 4 i=1 f i = 0 and M i,j = K + e i f j . Moreover, K equals the average entry size.
In particular, this tells us that L may be maximized by treating it as a function of just six variables, namely e 1 , e 2 , e 3 , f 1 , f 2 , f 3 , since e 4 , f 4 are also determined by these; changing K before solving this maximization problem simply has the impact of multiplying the entire matrix M that maximizes likelihood by a scalar.
Let E be the deviation matrix associated to M , where E i,j = e i f j .
Question 6.6. Another natural question to ask, in light of this corollary, is whether the matrix of rank at most r maximizing L is expressible as the sum of a rank one matrix and a matrix of rank at most r − 1 that maximizes L among matrices of rank at most r − 1.
Remarks When we consider matrices with fixed row and column sums, then we may ignore the denominator in the likelihood function and simply maximize the numerator.
Corollary 6.7. If M which maximizes L has e i = e j , then it also has f i = f j . Consequently, if it has e i = e j , then it also has f i = f j .
Proof. One consequence of having equal row and column sums is that it allows the likelihood function to be split into a product of four functions, one for each row, or else one for each column; this is because the sum of all table entries equals the sum of those in any row or column multiplied by four, allowing the denominator to be written just using variables from any one row or column. Thus, once the vector e is chosen, we find the best possible f for this given e by solving four separate maximization problems, one for each f i , i.e. one for each column. Setting e i = e j causes the likelihood function for column i to coincide with the likelihood function for column j, so both are maximized at the same value, implying f i = f j .
Next we prove a slightly stronger general fact for matrices in which rows and columns i, j may simultaneously be swapped without changing the data matrix: Proposition 6.8. If a matrix M maximizing likelihood has e i > e j > 0, then it also has f i > f j > 0.
Question 6.9. Does having a data matrix which is symmetric with respect to transpose imply that matrices maximizing likelihood will also be symmetric with respect to transpose?
Perhaps this could also be verified again by averaging, similarly to what we suggest for involutions swapping a pair of rows and columns simultaneously.

The 100 Swiss Francs Problem
We use the results derived to far so show how to reduce the 100 Swiss Francs problem to Question 6.1. Thus, an affirmative answer to Question 6.1 would provide a mathematical proof formally that the three tables in 2 a) are global maxima of the log-likelihood function for the basic LC model with r = 2 and data given in (9). Theorem 6.10. If the answer to Question 6.1 is yes, then the 100 Swiss francs problem is solved.
Proof. Proposition 6.3 showed that for M to maximize L, M must have row and column sums which are all equal to the quantity which we call R 1 , R 2 , R 3 , R 4 , C 1 , C 2 , C 3 , or C 4 at our convenience. The denominator of L may therefore be expressed as (4C 1 ) 10 (4C 2 ) 10 (4C 3 ) 10 (4C 4 ) 10 or as (4R 1 ) 10 (4R 2 ) 10 (4R 3 ) 10 (4R 4 ) 10 , enabling us to rewrite L as a product of four smaller functions using distinct sets of variables.
Note that letting S 4 simultaneously permute rows and columns will not change L, so let us assume the first two rows of M are linearly independent. Moreover, we may choose the first two rows in such a way that the next two rows are each nonnegative combinations of the first two. Since row and column sums are all equal, the third row, denoted v 3 , is expressible as xv 1 + (1 − x)v 2 for v 1 , v 2 the first and second rows and x ∈ [0, 1]. One may check that M does not have any row or column with values all equal to each other, because if it had one, then it would have the other, reducing to a three by three problem which one may solve, and one may check that the answer does not have as high of likelihood as 3 3 2 2 3 3 2 2 2 2 3 3 2 2 3 3 . Proposition 6.11 will show that if the answer to Question 6.1 is yes, then for M to maximize L, we must have x = 0 or x = 1, implying row 3 equals either row 1 or row 2, and likewise row 4 equals one of the first two rows. Proposition 6.12 shows M does not have three rows all equal to each other, and therefore must have two pairs of equal rows. Thus, the first column takes the form (a, a, b, b) T , so it is simply a matter of optimizing a and b, then noting that the optimal choice will likewise optimize the other columns (by virtue of the way we broke L into a product of four expressions which are essentially the same, one for each column). Thus, M takes the form a a b b a a b b b b a a b b a a since this matrix does indeed have rank two. Proposition 6.13 shows that to maximize L one needs 2a = 3b, finishing the proof.
Proposition 6.11. If the answer to Question 6.1 is yes, then row 3 equals either row 1 or row 2 in any matrix M which maximizes likelihood. Similarly, each row i with i > 2 equals either row 1 or row 2.
Proof. M 3,3 = xM 1,3 + (1 − x)M 2,3 for some x ∈ [0, 1], so M 3,3 ≤ max(M 1,3 , M 2,3 ). If M 1,3 = M 2,3 , then all entries of this column are equal, and one may use calculus to eliminate this possibility as follows: either M has rank one, and then we may replace column three by (c, c, 2c, c) T for suitable constant c to increase likelihood, since this only increases rank to at most two, or else the column space of M is spanned by (1, 1, 1, 1) T and some (a 1 , a 2 , a 3 , a 4 ) with a i = 0; specifically, column three equals (1/4, 1/4, 1/4, 1/4) + x(a 1 , a 2 , a 3 , a 4 ) for some x, allowing its contribution to the likelihood function to be expressed as a function of x whose derivative at x = 0 is nonzero, provided that a 3 = 0, implying that adding or subtracting some small multiple of (a 1 , a 2 , a 3 , a 4 ) T to the column will make the likelihood increase. If a 3 = 0, then row three is also constant, i.e. e 3 = f 3 = 0. But then, an affirmative answer to the second part of Question 6.1 will imply that this matrix does not maximize likelihood.
Suppose, on the other hand, M 1,3 > M 2,3 . Our goal then is to show x = 1. By Proposition 6.3 applied to columns rather than rows, we know that (1, 1, 1, 1) is in the span of the rows, so each row may be written as 1/4(1, 1, 1, 1) + cv for some fixed vector v whose coordinates sum to 0. Say row 1 equals 1/4(1, 1, 1, 1) + kv for k = 1. Writing row three as 1/4(1, 1, 1, 1) + lv, what remains is to rule out the possibility l < k. However, Proposition 6.8 shows that l < k and a 1 < a 3 together imply that swapping columns one and three will yield a new matrix of the same rank with larger likelihood. Now we turn to the case of l < k and a 1 ≥ a 3 . If a 1 = a 3 then swapping rows one and three will increase likelihood. Assume a 1 > a 3 . By Corollary 6.5, we have (e 1 , e 2 , e 3 , e 4 ) with e 1 > e 3 and (f 1 , f 2 , f 3 , f 4 ) with f 1 > f 3 . Therefore, if the answer to Question 6.1 is yes, then replacing e 1 , e 3 each by e 1 +e 3 2 and f 1 , f 3 each by f 1 +f 3 2 yields a matrix with larger likelihood, completing the proof. but then the fact that M maximizes L ensures d = f = h and c = e = g since L is a product of four expressions, one for each column, so that the second, third and fourth columns will all maximize their contribution to L in the same way. Since all row and column sums are equal, simple algebra may be used to show that all entries must be equal. However, we have already shown that such matrices do not maximize L.
Proposition 6.13. To maximize M requires a, b related by 2a = 3b.

Conclusions
In this paper we have reconsidered the classical latent class model for contingency table data and studied its geometric and statistical properties. For the former we have exploited tools from algebraic geometry and computation tools that have allowed us to display the complexities of the latent class model. We have focused on the problem of maximum likelihood estimation under LC models and have studied the singularities arising from symmetries in the contingency table data and the multiple maxima that appear to result from these. We have given an informal characterization of this problem, but a strict mathematical proof of the existence of identical multiple maxima has eluded us; we describe elements of a proof in a separate section.
We have also applied LC models to data arising in two real-life applications. In one, the model is quite simple and maximum likelihood estimation poses little problems, whereas in the other high-dimensional example various issues, computational as well as model-based, arise. From computational standpoint, both the EM and the Newton-Raphson algorithm are especially vulnerable to problems of multimodality and provide little in the way of clues regarding the dimensionality difficulties associated with the underlying structure of LC models. Furthermore, the seemingly singular behavior of the Fisher information matrix at the MLE that we observe even for well-behaved, identifiable models is an additional element of complexity.
Based on our work, we would advise practitioners to exercise caution in applying LC models, especially to sparse data. They have a tremendous heuristic appeal and in some examples provide a clear and convincing description of the data. But in many situations, the kind of complex behavior explored in this paper may lead to erroneous inferences.