poLCA : An R Package for Polytomous Variable Latent Class Analysis

poLCA is a software package for the estimation of latent class and latent class regression models for polytomous outcome variables, implemented in the R statistical computing environment. Both models can be called using a single simple command line. The basic latent class model is a ﬁnite mixture model in which the component distributions are assumed to be multi-way cross-classiﬁcation tables with all variables mutually independent. The latent class regression model further enables the researcher to estimate the eﬀects of covariates on predicting latent class membership. poLCA uses expectation-maximization and Newton-Raphson algorithms to ﬁnd maximum likelihood estimates of the model parameters.


Introduction
Latent class analysis is a statistical technique for the analysis of multivariate categorical data.When observed data take the form of a series of categorical responses-as, for example, in public opinion surveys, individual-level voting data, studies of inter-rater reliability, or consumer behavior and decision-making-it is often of interest to investigate sources of confounding between the observed variables, identify and characterize clusters of similar cases, and approximate the distribution of observations across the many variables of interest.Latent class models are a useful tool for accomplishing these goals.
The latent class model seeks to stratify the cross-classification table of observed (or, "manifest") variables by an unobserved ("latent") unordered categorical variable that eliminates all confounding between the manifest variables.Conditional upon values of this latent variable, responses to all of the manifest variables are assumed to be statistically independent; an assumption typically referred to as "conditional" or "local" independence.The model, in effect, probabilistically groups each observation into a "latent class," which in turn produces expectations about how that observation will respond on each manifest variable.Although the model does not automatically determine the number of latent classes in a given data set, it does offer a variety of parsimony and goodness of fit statistics that the researcher may use in order to make a theoretically and empirically sound assessment.
Because the unobserved latent variable is unordered categorical, the latent class model is actually a type of finite mixture model.The component distributions in the mixture are cross-classification tables of equal dimension to the observed table of manifest variables, and, following the assumption of conditional independence, the frequency in each cell of each component table is simply the product of the respective class-conditional marginal frequencies (the parameters estimated by the latent class model are the proportion of observations in each latent class, and the probabilities of observing each response to each manifest variable, conditional on latent class).A weighted sum of these component tables forms an approximation of the distribution of cases across the cells of the observed table.Observations with similar sets of responses on the manifest variables will tend to cluster within the same latent classes.The model may also be fit to manifest variables that are ordinal, but they will be treated as nominal.In practice, this does not usually restrict analyses in any meaningful way.
An extension of this basic model permits the inclusion of covariates to predict latent class membership.Whereas in the basic model, every observation has the same probability of belonging to each latent class prior to observing the responses to the manifest variables, in the more general latent class "regression" model, these prior probabilities vary by individual as a function of some set of independent (or, "concomitant") variables.
Examples of latent class models in political science include McCutcheon (1985); Feick (1989); Breen (2000); Hill and Kriesi (2001a,b); Blaydes and Linzer (2006).The latent class model is similar to the latent trait model widely used for the estimation of voter and legislator "ideal points" (see, for example, Clinton, Jackman, and Rivers (2004)) in that both models assume the presence of an underlying, unobserved latent variable to explain patterns among observed variables.Yet whereas the latent trait model assumes that the latent variable is continuous, the latent class model assumes that the latent variable is categorical.Moreover, unlike the ideal point model, the latent class model requires no assumptions about respondent utility functions, utility maximization, or rationality.
poLCA is the first R package to enable the user to estimate latent class models for manifest variables with any number of possible outcomes, and it is the only package that estimates latent class regression models with covariates (Linzer and Lewis 2007;R Development Core Team 2007).The two other R commands that currently exist to estimate latent class modelsthe lca command in package e1071, and the gllm command in package gllm-can only estimate the basic model for dichotomous outcome variables.
Note that there is occasionally some confusion over the term "latent class regression" (LCR); in practice it can have two meanings.In poLCA, LCR models refer to latent class models in which the probability of latent class membership is predicted by one or more covariates.In other contexts, however, LCR is used to refer to regression models in which the dependent variable is partitioned into latent classes as part of estimating the regression model.It is a way to simultaneously fit more than one regression to the data when the latent data partition is unknown.The regmix command in package fpc will estimate this other type of LCR model, as will the flexmix command in package flexmix.Because of these terminology issues, the LCR models estimated using poLCA are sometimes termed "latent class models with covariates" or "concomitant-variable latent class analysis," both of which are accurate descriptions of this model.

Latent class models
The basic latent class model is a finite mixture model in which the component distributions are assumed to be multi-way cross-classification tables with all variables mutually independent.This model was originally proposed by Lazarsfeld (1950) under the name "latent structure analysis."Chapter 13 in Agresti (2002) details the connection between latent class models and finite mixture models.

Terminology and model definition
Suppose we observe J polytomous categorical variables (the "manifest" variables), each of which contains K j possible outcomes, for individuals i = 1...N .The manifest variables may have different numbers of outcomes, hence the indexing by j.Denote as Y ijk the observed values of the J manifest variables such that Y ijk = 1 if respondent i gives the kth response to the jth variable, and Y ijk = 0 otherwise, where j = 1 . . .J and k = 1 . . .K j .
The latent class model approximates the observed joint distribution of the manifest variables as the weighted sum of a finite number, R, of constituent cross-classification tables.R is fixed prior to estimation on the basis of either theoretical reasons or model fit; this issue is addressed in greater detail in Section 2.4 below.Let π jrk denote the class-conditional probability that an observation in class r = 1 . . .R produces the kth outcome on the jth variable.Within each class, for each manifest variable, therefore, K j k=1 π jrk = 1.Further denote as p r the R mixing proportions that provide the weights in the weighted sum of the component tables, with r p r = 1.The probability that an individual i in class r produces a particular set of J outcomes on the manifest variables, assuming local independence, is the product (1) The probability density function across all classes is the weighted sum The parameters estimated by the latent class model are p r and π jrk .
Given estimates pr and πjrk of p r and π jrk , respectively, the posterior probability that each individual belongs to each class, conditional on the observed values of the manifest variables, can be calculated using Bayes' formula: Recall that the πjrk are estimates of outcome probabilities conditional on class r.
It is important to remain aware that the number of independent parameters estimated by the latent class model increases rapidly with R, J, and K j .Given these values, the number of parameters is R j (K j − 1) + (R − 1).If this number exceeds either the total number of observations, or one fewer than the total number of cells in the cross-classification table of the manifest variables, then the latent class model will be unidentified.

Parameter estimation
poLCA estimates the latent class model by maximizing the log-likelihood function with respect to p r and π jrk , using the expectation-maximization (EM) algorithm (Dempster, Laird, and Rubin 1977).This log-likelihood function is identical in form to the standard finite mixture model log-likelihood.As with any finite mixture model, the EM algorithm is applicable because each individual's class membership is unknown and may be treated as missing data (McLachlan and Krishnan 1997;McLachlan and Peel 2000).
The EM algorithm proceeds iteratively.Begin with arbitrary initial values of pr and πjrk , and label them pold r and πold jrk .In the expectation step, calculate the "missing" class membership probabilities using Eq. 3, substituting in pold r and πold jrk .In the maximization step, update the parameter estimates by maximizing the log-likelihood function given these posterior P(r|Y i ), with as the new prior probabilities and as the new class-conditional outcome probabilities (see Everitt and Hand (1981); Everitt (1984)).In Eq. 6, πnew jr is the vector of length K j of class-r conditional outcome probabilities for the jth manifest variable; and Y ij is the N × K j matrix of observed outcomes Y ijk on that variable.The algorithm repeats these steps, assigning the new to the old, until the overall log-likelihood reaches a maximum and ceases to increment beyond some arbitrarily small value.
poLCA takes advantage of the iterative nature of the EM algorithm to make it possible to estimate the latent class model even when some of the observations on the manifest variables are missing.Although poLCA does offer the option to listwise delete observations with missing values before estimating the model, it is not necessary to do so.Instead, when determining the product in Eq. 1 and the sum in the numerator of Eq. 6, poLCA simply excludes from the calculation any manifest variables with missing observations.The priors are updated in Eq. 3 using as many or as few manifest variables as are observed for each individual.
Depending on the initial values chosen for pold r and πold jrk , and the complexity of the latent class model being estimated, the EM algorithm may only find a local maximum of the log-likelihood function, rather than the desired global maximum.For this reason, it is always advisable to re-estimate a particular model a couple of times when using poLCA, in an attempt to find the global maximizer to be taken as the maximum likelihood solution.

Standard error estimation
poLCA estimates standard errors of the estimated class-conditional response probabilities πjrk and the mixing parameters pr using the empirical observed information matrix (Meilijson 1989), which, following McLachlan and Peel (2000, 66), equals where s(Y i ; Ψ) is the score function with respect to the vector of parameters Ψ for the ith observation, evaluated at the maximum likelihood estimate Ψ; where θ ir = P(r|Y i ) is the posterior probability that observation i belongs to class r (Eq.3).
The covariance matrix of the parameter estimates is then approximated by the inverse of I e ( Ψ; Y ).
Because of the sum-to-one constraint on the π jrk across each manifest variable, it is useful to reparameterize the score function in terms of log-ratios φ jrk = ln(π jrk /π jr1 ) for given outcome variable j and class r.Then, for the lth response on the hth item in the qth class, Likewise, denoting ω r = ln(p r /p 1 ), then for the log-ratio corresponding to the qth mixing parameter, s(Y i ; ω q ) = θ iq − p q . (10) To transform the covariance matrix of these log-ratios back to the original units of π and p, we apply the delta method.For the response probabilities, let g(φ jrk ) = π jrk = e φ jrk / l e φ jrl .Taking as VAR( φ) the submatrix of the inverse of I e ( Ψ; Y ) corresponding to the φ parameters, then where g (φ) is the Jacobian consisting of elements and h = j and l = k.
For the mixing parameters, similarly let h(ω r ) = p r = e pr / q e pq .Taking as VAR(ω) the submatrix of the inverse of I e ( Ψ; Y ) corresponding to the ω parameters, then where h (ω) is the Jacobian consisting of elements Standard errors of each parameter estimate are equal to the square root of the values along the main diagonal of covariance matrices VAR(π) and VAR(p).

Model selection and goodness of fit criteria
One of the benefits of latent class analysis, in contrast to other statistical techniques for clustered data, is the variety of tools available for assessing model fit and determining an appropriate number of latent classes R for a given data set.In some applications, the number of latent classes will be selected for primarily theoretical reasons.In other cases, however, the analysis may be of a more exploratory nature, with the objective being to locate the best fitting or most parsimonious model.The researcher may then begin by fitting a complete "independence" model with R = 1, and then iteratively increasing the number of latent classes by one until a suitable fit has been achieved.
Adding an additional class to a latent class model will increase the fit of the model, but at the risk of fitting to noise, and at the expense of estimating a further 1 + j (K j − 1) model parameters.Parsimony criteria seek to strike a balance between over-and under-fitting the model to the data by penalizing the log-likelihood by a function of the number of parameters being estimated.The two most widely used parsimony measures are the Bayesian information criterion, or BIC (Schwartz 1978) and Akaike information criterion, or AIC (Akaike 1973).Preferred models are those that minimize values of the BIC and/or AIC.Let Λ represent the maximum log-likelihood of the model and Φ represent the total number of estimated parameters.Then, AIC = −2Λ + 2Φ and BIC = −2Λ + Φ ln N.
poLCA calculates these parameters automatically when estimating the latent class model.The BIC will usually be more appropriate for basic latent class models because of their relative simplicity (Lin and Dayton 1997;Forster 2000).
Calculating Pearson's χ 2 goodness of fit and likelihood ratio chi-square (G 2 ) statistics for the observed versus predicted cell counts is another method to determine how well a particular model fits the data (Goodman 1970).Let q c denote the observed number of cases in each of the C = K j cells of the cross-classification table of the manifest variables.Let Q denote the expected number of cases in each cell under a given model.The cth cell (where c = 1 . . .C) corresponds to one particular sequence of J outcomes on the manifest variables.Taking the πjrk corresponding only to those outcomes, πjrk .
The two test statistics are then calculated as Like the AIC and BIC, these statistics are outputted automatically after calling poLCA.
Generally, the goal is to select models that minimize χ 2 or G 2 without estimating excessive numbers of parameters.Note, however, that the distributional assumptions for these statistics are not met if many cells of the observed cross-classification table contain very few observations.Indeed, common practice holds that no fewer than 10-20% of the cells should contain fewer than five observations if either chi-square test is to be used.

Latent class regression models
The latent class regression model generalizes the basic latent class model by permitting the inclusion of covariates to predict individuals' latent class membership (Dayton and Macready 1988;Hagenaars and McCutcheon 2002).This is a so-called "one-step" technique for estimating the effects of covariates, because the coefficients on the covariates are estimated simultaneously as part of the latent class model.An alternate estimation procedure that is sometimes used is called the "three-step" model: estimate the basic latent class model, calculate the predicted posterior class membership probabilities using Eq. 3, and then use these values as the dependent variable(s) in a regression model with the desired covariates.However, as demonstrated by Bolck, Croon, and Hagenaars (2004), the three-step procedure produces biased coefficient estimates.It is preferable to estimate the entire latent class regression model all at once.
Covariates are included in the latent class regression model through their effects on the priors p r .In the basic latent class model, it is assumed that every individual has the same prior probabilities of latent class membership.The latent class regression model, in contrast, allows individuals' priors to vary depending upon their observed covariates.

Terminology and model definition
Denote the mixing proportions in the latent class regression model as p ri to reflect the fact that these priors are now free to vary by individual.It is still the case that r p ri = 1 for each individual.To accommodate this constraint, poLCA employs a generalized (multinomial) logit link function for the effects of the covariates on the priors (Agresti 2002).
Let X i represent the observed covariates for individual i. poLCA arbitrarily selects the first latent class as a "reference" class and assumes that the log-odds of the latent class membership priors with respect to that class are linear functions of the covariates.Let β r denote the vector of coefficients corresponding to the rth latent class.With S covariates, the β r have length S + 1; this is one coefficient on each of the covariates plus a constant.Because the first class is used as the reference, β 1 = 0 is fixed by definition.Then, Following some simple algebra, this produces the general result that The parameters estimated by the latent class regression model are the R − 1 vectors of coefficients β r and, as in the basic latent class model, the class-conditional outcome probabilities π jrk .Given estimates βr and πjrk of these parameters, the posterior class membership probabilities in the latent class regression model are obtained by replacing the p r in Eq. 3 with the function p r (X i ; β) in Eq. 11: The number of parameters estimated by the latent class regression model is equal to R j (K j − 1) + (S + 1)(R − 1).The same considerations mentioned earlier regarding model identifiability also apply here.

Parameter estimation
The latent class regression model log-likelihood function is identical to Eq. 4 except that the function p r (X i ; β) (Eq.11) takes the place of p r : To find the values of βr and πjrk that maximize this function, poLCA uses a modified EM algorithm with a Newton-Raphson step, as set forth by Bandeen-Roche, Miglioretti, Zeger, and Rathouz (1997).This estimation process begins with initial values of βold r and πold jrk that are used to calculate posterior probabilities P(r|X i ; Y i ) (Eq. 12).The coefficients on the concomitant variables are updated according to the formula where D β is the gradient and D 2 β the Hessian matrix with respect to β.The πnew jrk are updated as These steps are repeated until convergence, assigning the new parameter estimates to the old in each iteration.The formulas for the gradient and Hessian matrix are provided in Bandeen-Roche et al. (1997).
Because all of the concomitant variables must be observed in order to calculate p ri (Eq.11), poLCA listwise deletes cases with missing values on the X i before estimating the latent class regression model.However, missing values on the manifest variables Y i can be accommodated in the latent class regression model, just as they were in the basic latent class model.
Note that when employing this estimation algorithm, different initial parameter values may lead to different local maxima of the log-likelihood function.To be more certain that the global maximum likelihood solution has been found, the poLCA function call should always be repeated a handful of times.

Standard error estimation
For latent class models with covariates, standard errors are obtained just as for models without covariates: using the empirical observed information matrix (Eq.7).First, we generalize the score function (Eq.8) so that As before, θ ir denote posterior probabilities.Since this function is no different than Eq. 8 in terms of the π parameters, the score function s(X i , Y i ; φ hql ) = s(Y i ; φ hql ) (Eq. 9), and the covariance matrix VAR(π) may be calculated in precisely the same way as for models without covariates.
Now, however, the priors p ri are free to vary by individual as a function of some set of coefficients β, as given in Eq. 11.Letting q index classes and s index covariates, The standard errors of the coefficients β are equal to the square root of the values along the main diagonal of the submatrix of the inverse of the empirical observed information matrix corresponding to the β parameters.(Note that when the model has no covariates, X i = 1 and p iq = p q (that is, the priors do not vary by individual), so Eq. 17 reduces to Eq. 10 as expected.) To obtain the covariance matrix of the mixing parameters p r , which are the average value across all observations of the priors p ir , we apply the delta method.Let where h (β) is a Jacobian with elements Although poLCA does not automatically compute p-values for the regression coefficients, these may be calculated by the user, from the coefficients and standard errors, in the usual manner.

Using poLCA
The poLCA package makes it possible to estimate a wide range of latent class models in R using a single command line, poLCA.Also included in the package is the command poLCA.simdata,which enables the user to create simulated data sets that match the data-generating process assumed by either the basic latent class model or the latent class regression model.This functionality is useful for testing the poLCA estimator and for performing Monte Carlo-style analyses of latent class models.

Data input and sample data sets
Data are input to the poLCA function as a data frame containing all manifest and concomitant variables (if needed).The manifest variables must be coded as integer values starting at one for the first outcome category, and increasing to the maximum number of outcomes for each variable.If any of the manifest variables contain zeros, negative values, or decimals, poLCA will produce an error message and terminate without estimating the model.The input data frame may contain missing values.
poLCA also comes pre-installed with five sample data sets that are useful for exploring different aspects of latent class and latent class regression models.carcinoma: Dichotomous ratings by seven pathologists of 118 slides for the presence or absence of carcinoma in the uterine cervix.Source: Agresti (2002, 542).
cheating: Dichotomous responses by 319 undergraduates to questions about cheating behavior.Also each student's GPA, which is useful as a concomitant variable.Source: Dayton (1998, 33 and 85).
election: Two sets of six questions with four responses each, asking respondents' opinions of how well various traits describe presidential candidates Al Gore and George W. Bush.Also potential covariates vote choice, age, education, gender, and party ID.Source: The National Election Studies (2000).
gss82: Attitudes towards survey taking across two dichotomous and two trichotomous items among 1202 white respondents to the 1982 General Social Survey.Source: McCutcheon (1987, 30).
These data sets may be accessed using the data(name ) command.Examples of models and analyses using the sample data sets are included in the internal documentation for each.

poLCA command line options
To specify a latent class model, poLCA uses the standard, symbolic R model formula expression.The response variables are the manifest variables of the model.Because latent class models have multiple manifest variables, these variables must be"bound"as cbind(Y1,Y2,Y3,...) in the model formula.For the basic latent class model with no covariates, the formula definition takes the form The 1 instructs poLCA to estimate the basic latent class model.For the latent class regression model, replace the 1 with the desired function of covariates, as, for example: Further assistance on formula specification in R can be obtained by entering ?formula at the command prompt.
To estimate the specified latent class model, the default poLCA command is: > poLCA(formula, data, nclass=2, maxiter=1000, graphs=FALSE, tol=1e-10, na.rm=TRUE, probs.start=NULL) At minimum, it is necessary to enter a formula (as just described) and a data frame (as described in the previous subsection).The remaining options are: nclass: The number of latent classes to assume in the model; R in the above notation.Setting nclass=1 results in poLCA estimating the loglinear independence model (Goodman 1970).The default is two.
maxiter: The maximum number of iterations through which the estimation algorithm will cycle.If convergence is not achieved before reaching this number of iterations, poLCA terminates and reports an error message.The default is 1000, but this will be insufficient for certain models.
graphs: Logical, for whether poLCA should graphically display the parameter estimates at each stage of the updating algorithm.The default is FALSE, as setting this option to TRUE slows down the estimation process.
tol: A tolerance value for judging when convergence has been reached.When the oneiteration change in the estimated log-likelihood is less than tol, the estimation algorithm stops updating and considers the maximum log-likelihood to have been found.The default is 1 × 10 −10 which is a standard value; this option will rarely need to be invoked.

poLCA output
The poLCA function returns an object containing the following elements: y: A data frame of the manifest variables.
x: A data frame of the covariates, if specified.
N: Number of cases used in the model.

Nobs: Number of fully observed cases (less than or equal to N).
probs: A list of matrices containing the estimated class-conditional outcome probabilities πjrk .Each item in the list represents one manifest variable; columns correspond to possible outcomes on each variable, and rows correspond to the latent classes.
probs.se:Standard errors of the estimated class-conditional response probabilities, in the same format as probs.

P:
The respective size of each latent class; equal to the estimated mixing proportions pr in the basic latent class model, or the mean of the priors in the latent class regression model.

P.se:
The standard errors of P.
posterior: An N ×R matrix containing each observation's posterior class membership probabilities.
predclass: A vector of length N of predicted class memberships, by modal assignment.
predcell: A table of observed versus predicted cell counts for cases with no missing values.

llik:
The maximum value of the estimated model log-likelihood.
numiter: The number of iterations required by the estimation algorithm to achieve convergence.
coeff: An (S + 1) × (R − 1) matrix of estimated multinomial logit coefficients βr , for the latent class regression model.Rows correspond to concomitant variables X. Columns correspond to the second through Rth latent classes; see Eq. 11.
coeff.se:Standard errors of the coefficient estimates, in the same format as coeff.
coeff.V: Covariance matrix of the coefficient estimates.
Chisq: Pearson Chi-square goodness of fit statistic.
time: Length of time it took to estimate the model.

npar:
The number of degrees of freedom used by the model (that is, the number of estimated parameters).
resid.df:The number of residual degrees of freedom, equal to the lesser of N and ( j K j )− 1, minus npar.
eflag: Logical, error flag.True if estimation algorithm needed to automatically restart with new initial parameters, otherwise false.A restart is caused in the event of either a noninvertible Hessian matrix, or computational/rounding errors that result in nonsensical parameter estimates.If one of these errors occurs, poLCA outputs an error message to alert the user.
probs.start:A list of matrices containing the class-conditional response probabilities used as starting values in the EM estimation algorithm.If the algorithm needed to restart (see eflag), this contains the starting values used for the final, successful, run of the estimation algorithm.
Selected items from this list are outputted automatically once the specified latent class model has been estimated.

Reordering the latent classes
Because the latent classes are unordered categories, the numerical order of the estimated latent classes in the model output is arbitrary, and is determined solely by the start values of the EM algorithm.If probs.start is set to NULL (the default) when calling poLCA, then the function generates the starting values randomly in each run.This means that repeated runs of poLCA will typically produce results containing the same parameter estimates (corresponding to the same maximum log-likelihood), but with reordered latent class labels.
To manually "fix" the order of the estimated latent classes, run poLCA once, and extract the outputted list of probs.start.
> lc <-poLCA(f,dat,nclass=2) > probs.start<-lc$probs.start Then rearrange the order of the rows (corresponding to the latent classes) in each matrix of that list, to match the desired ordering of the outputted latent classes.For example, suppose you have estimated a two-class model and wish to reverse the class labels in the output.

Recognizing and avoiding local maxima
A well-known drawback of the EM algorithm is that depending upon the initial parameter values chosen in the first iteration, the algorithm may only find a local, rather than the global, maximum of the log-likelihood function (McLachlan and Krishnan 1997).To avoid these local maxima, a user should always call poLCA at least a couple of times to attempt to locate the estimated model parameters that correspond to the model with the global maximum likelihood.
We demonstrate this using a basic three-class latent class model to analyze the four survey variables in the gss82 data set included in the poLCA package.
> data(gss82) > f <-cbind(PURPOSE,ACCURACY,UNDERSTA,COOPERAT)~1 We estimate this model 500 times, and after each function call, we record the maximum log-likelihood and the estimated population sizes of the three types of survey respondent.
Following McCutcheon (1987), from whom these data were obtained, we label the three types ideal, skeptics, and believers.Among other characteristics, the ideal type is the most likely to have a good understanding of surveys, while the believer type is the least likely.
Results of this simulation are reported in Table 1.Of the five local maxima of the loglikelihood function that were found, the global maximum was obtained in only approximately half of the trials.At the global maximum, the ideal type is estimated to represent 62.1% of the population, with another 17.2% skeptics and 20.7% believers.In contrast, the second-most frequent local maximum was also the lowest of the local maxima, and the parameter estimates corresponding to that "solution" are substantially different: 29.7% ideal types, 53.3% skeptics, and 17.0% believers.This is why it is essential to run poLCA multiple times until you can be reasonably certain that you have found the parameter estimates that produce the global maximum likelihood solution.

Creating simulated data sets
The command poLCA.simdatawill generate simulated data sets that can be used to examine properties of the latent class and latent class regression model estimators.The properties of the simulated data set are fully customizable, but poLCA.simdatauses the following default arguments in the function call.
probs: A list of matrices of dimension nclass × nresp, containing, by row, the classconditional outcome probabilities π jrk (which must sum to 1) for the manifest variables.Each matrix represents one manifest variable.If probs is NULL (default) then the outcome probabilities are generated randomly.
nclass: The number of latent classes, R. If probs is specified, then nclass is set equal to the number of rows in each matrix in that list.If classdist is specified, then nclass is set equal to the length of that vector.Otherwise, the default is two.
ndv: The number of manifest variables, J.If probs is specified, then ndv is set equal to the number of matrices in that list.If nresp is specified, then ndv is set equal to the length of that vector.Otherwise, the default is four.

nresp:
The number of possible outcomes for each manifest variable, K j , entered as a vector of length ndv.If probs is specified, then ndv is set equal to the number of columns in each matrix in that list.If both probs and nresp are NULL (default), then the manifest variables are assigned a random number of outcomes between two and five.
x: A matrix of concomitant variables, of dimension N × niv.If niv > 0 but x is NULL (default) then the concomitant variable(s) will be generated randomly.If both x and niv are entered, then then the number of columns in x overrides the value of niv.
niv: The number of concomitant variables, S. Setting niv = 0 (default) creates a data set assuming no covariates.If nclass=1 then niv is automatically set equal to 0. Unless x is specified, all covariates consist of random draws from a standard normal distribution and are mutually independent.
classdist: A vector of mixing proportions of length nclass, corresponding to p r .classdist must sum to 1. Disregarded if niv>1 because then classdist is, in part, a function of the concomitant variables.If classdist is NULL (default), then the p r are generated randomly.
missval: Logical.If TRUE then a fraction pctmiss of the observations on the manifest variables are randomly dropped as missing values.Default is FALSE.
pctmiss: The percentage of values to be dropped as missing, if missval=TRUE.If pctmiss is NULL (default), then a value between 5% and 40% is chosen randomly.
Note that in many instances, specifying values for certain arguments will override other specified arguments.Be sure when calling poLCA.simdatathat all arguments are in logical agreement, or else the function may produce unexpected results.
Specifying the list of matrices probs can be tricky; we recommend a command structure such as this for, for example, five manifest variables, three latent classes, and K j = (3, 2, 3, 4, 3).
probs: A list of matrices of dimension nclass × nresp containing the class-conditional outcome probabilities.
nresp: A vector containing the number of possible outcomes for each manifest variable.
b: A matrix containing the coefficients on the covariates, if used.
classdist: The mixing proportions corresponding to each latent class.
pctmiss: The percent of observations missing.
trueclass: A vector of length N containing the "true" class membership for each individual.
Examples of possible uses of poLCA.simdataare included in the poLCA internal documentation and may be accessed by entering ?poLCA.simdata in R. One example demonstrates that even when the "true" data generating process involves a series of covariates-so that each observation has a different prior probability of belonging to each class-the posterior probabilities of latent class membership can still be recovered with high accuracy using a basic model specified without covariates.A second example confirms that in data sets with missing values, the poLCA function produces consistent estimates of the class-conditional response probabilities π jrk regardless of whether the researcher elects to include or listwise delete the observations with missing values.

Two examples
To illustrate the usage of the poLCA package, we present two examples: a basic latent class model and a latent class regression model, using sample data sets included in the package.

Basic latent class modeling with the carcinoma data
The carcinoma data from Agresti (2002, 542) consist of seven dichotomous variables that represent the ratings by seven pathologists of 118 slides on the presence or absence of carcinoma.
The purpose of studying these data is to model "interobserver agreement" by examining how subjects might be divided into groups depending upon the consistency of their diagnoses.It is straightforward to replicate Agresti's published results (Agresti 2002, 543) using the series of commands: The four-class model will typically require a larger number of iterations to achieve convergence.
Figure 1 shows a screen capture of the estimation of model lc3 with the graphs option set to TRUE and maxiter=400.In this case, the model has converged after 34 iterations.As Agresti describes, the three estimated latent classes clearly correspond to a pair of classes that are consistently rated negative (37%) or positive (44%), plus a third "problematic" class representing 18% of the population.In that class, pathologists B, E, and G tend to diagnose positive; C, D, and F tend to diagnose negative; and A is about 50/50.The full output from the estimation of model lc3 is given below.First, the estimated classconditional response probabilities πjrk are reported for pathologists A through G, with each row corresponding to a latent class, and each column corresponding to a diagnosis; negative in the first column, and positive in the second.Thus, for example, a slide belonging to the first ("negative") class has a 94% chance of being rated free from carcinoma by rater A, an 86% chance of the same from rater B, an 100% chance from rater C, and so forth.
Next, the output provides the estimated mixing proportions pr corresponding to the share of observations belonging to each latent class.These are the same values that appear in Figure 1.An alternate method for determining the size of the latent classes is to assign each observation to a latent class on an individual basis according to its model posterior class membership probability.Values using this technique are reported directly below the estimated mixing proportions.Congruence between these two sets of population shares often indicates a good fit of the model to the data.
The next set of results simply reports the number of observations, the number of fully observed cases (for data sets with missing values and na.rm=FALSE), the number of estimated parameters, residual degrees of freedom, and maximum log-likelihood.It is always worth checking to ensure that the number of residual degrees of freedom is non-negative; poLCA will output a warning message if this is the case.
Finally, poLCA outputs a number of goodness of fit statistics as described in Section 2.4.For the carcinoma data, the minimum AIC and BIC criteria both indicate that the three-class model is most parsimonious: with two classes, the AIC is 664.5 and the BIC is 706.1; with three classes, the AIC decreases to 633.4 and the BIC decreases to 697.1; and with four classes, the AIC increases again to 641.6 and the BIC increases to 727.5.

Latent class regression modeling with the election data
In the election data set, respondents to the 2000 American National Election Study public opinion poll were asked to evaluate how well a series of traits-moral, caring, knowledgable, good leader, dishonest, and intelligent-described presidential candidates Al Gore and George W. Bush.Each question had four possible choices: (1) extremely well; (2) quite well; (3) not too well; and (4) not well at all.

Models with one covariate
A reasonable theoretical approach might suppose that there are three latent classes of survey respondents: Gore supporters, Bush supporters, and those who are more or less neutral.Gore supporters will tend to respond favorably towards Gore and unfavorably towards Bush, with the reverse being the case for Bush supporters.Those in the neutral group will not have strong opinions about either candidate.We might further expect that falling into one of these three groups is a function of each individual's party identification, with committed Democrats more likely to favor Gore, committed Republicans more likely to favor Bush, and less intense partisans tending to be indifferent.We can investigate this hypothesis using a latent class regression model.
Begin by loading the election data into memory, and specifying a model with 12 manifest variables and PARTY as the lone concomitant variable.The PARTY variable is coded across seven categories, from strong Democrat at 1 to strong Republican at 7. People who primarily consider themselves Independents are at 3-4-5 on the scale.Next, estimate the latent class regression model and assign those results to object nes.party.> data(election) > f.party <-cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG, MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB)~PARTY > nes.party <-poLCA(f.party,election,nclass=3)By examining the estimated class-conditional response probabilities, we confirm that the model finds that the three groups indeed separate as expected, with 27% in the favor-Gore group, 34% in the favor-Bush group, and 39% in the neutral group.This example also illustrates a shortcoming of the χ 2 goodness of fit statistic, which is calculated to be over 34.5 billion.With only 1300 observations but nearly 17 million cells in the observed cross-classification table (that is, four responses to each of 12 questions, or 4 12 cells), the vast majority of the cells will contain zero cases.For models such as this, using the χ 2 statistic to assess model fit is not advised.
In addition to the information outputted for the basic model, the poLCA output now also includes the estimated coefficients βr on the covariates, and their standard errors.Here, the neutral group is the first latent class, the favor-Bush group is the second latent class, and the favor-Gore group is the third latent class.Following the terminology in Section 3.1, the log-ratio prior probability that a respondent will belong to the favor-Bush group with respect to the neutral group is ln(p 2i /p 1i ) = −3.82+ 0.79 × PARTY.Likewise, the log-ratio prior probability that a respondent will belong to the favor-Gore group with respect to the neutral group is ln(p 3i /p 1i ) = 1.16 − 0.57 × PARTY.Eq. 11 provides the formula for converting these log-ratios into predicted prior probabilities for each latent class.

===========================================
As expected, regardless of age, strong Democrats are very unlikely to belong to the Bushaffinity group, and strong Republicans are very unlikely to belong to the Gore-affinity group.However, it is interesting to observe that while strong Republicans in 2000 had extremely high levels of affinity for Bush at all ages, strong Democrats below the age of 30 tended to be just as (or more) likely to belong to the neutral group as to the Gore-affinity group.

Conclusion
The R package poLCA provides an easy-to-use framework for the estimation of latent class models and latent class regression models for the analysis of multivariate categorical data.The manifest variables may contain any number of possible (polytomous) outcomes.In the latent class regression model, covariates may be used to predict individual observations' latent class membership.poLCA also includes a tool for observing the iterative parameter estimation process and visualizing model results, and a function for creating simulated multivariate categorical data sets with an unobserved latent categorical structure.
poLCA is still undergoing active development.Planned extensions include flexibility to relax the assumption of local independence among selected manifest variables, explicit treatment of manifest variables as ordinal as well as nominal, incorporation of sampling weights, and accommodation of user-specified constraints on the class-conditional response probabilities π jrk as a way to simplify models, achieve model identifiability, test substantive hypotheses, and analyze model fit.Such constraints might, for example, require selected response probabilities to be set equal to one another across different classes, across manifest variables within classes, or equal to fixed constant values, as in Goodman (1974).This extension would also permit the estimation of so-called "simultaneous" latent class models across multiple groups where the groups are already known (or at least theorized) to exist in the data (Clogg and Goodman 1986).The researcher would then include in the model a manifest variable measuring this known categorization, and specify that one of the grouping variable response probabilities be fixed at 1.0 in each class.All future updates will be made available both through the online Comprehensive R Archive Network (http://cran.r-project.org) and the poLCA project Web site at http://userwww.service.emory.edu/~dlinzer/poLCA.

Figure 1 :
Figure 1: Estimation of the three-class basic latent class model using the carcinoma data; obtained by setting graphs=TRUE in the poLCA function call.Each group of red bars represents the conditional probabilities, by latent class, of being rated positively by each of the seven pathologists (labeled A through G).Taller bars correspond to conditional probabilities closer to 1 of a positive rating.

Figure 2 :
Figure 2: Predicted prior probabilities of latent class membership at varying levels of partisan self-identification. Results are from a three-class latent class regression model.

Figure 3 :
Figure 3: Predicted prior probabilities of latent class membership for strong Democrats (left) and strong Republicans (right) at ages 18-80.
jrk , to be used as the starting values for the EM estimation algorithm.Each matrix in the list corresponds to one manifest variable, with one row for each latent class, and one column for each possible outcome.The default is NULL, meaning that starting values are generated randomly.
na.rm: Logical, for how poLCA handles cases with missing values on the manifest variables.If TRUE, those cases are removed (listwise deleted) before estimating the model.If FALSE, cases with missing values are retained.(As discussed above, cases with missing covariates are always removed.)The default is TRUE.probs.start:A list of matrices of class-conditional response probabilities, π

Table 1 :
Results of 500 poLCA function calls for three-class model using gss82 data set.Five local maxima of the log-likelihood function were found.Estimated latent class proportions pr are reported for each respondent type at each local maximum.