randomLCA: An R Package for Latent Class with Random Effects Analysis

Ken J. Beath, Macquarie University

October 07, 2026

1 Introduction

Latent class models (Lazarsfeld and Henry 1968) are a method originally developed for sociology where they are used to identify clusters or sub-groups of subjects, based on multivariate binary observations, and as such are a form of finite mixture model. Their application has been further expanded into many areas, such as psychology, market research and medicine. Diagnostic classification using latent class methods has been applied in a number of areas, with early applications by Golden (1982) to dementia, Young (1982) to develop diagnostic criteria for schizophrenia and Rindskopf and Rindskopf (1986) for myocardial infarction. An advantage of using latent class analysis over other classification methods is that the classification is model based, allowing use of model selection techniques to determine which classification scheme is most appropriate. This compares with classifications developed from simple observation, that may give undue weight to one or more symptoms or outcomes. An example of the problems with this type of analysis is demonstrated by Nyholt et al. (2004) who used latent class analysis of headache symptoms to show that classification of migraine with and without aura as separate diagnoses is not supported. While latent class methods have been extended to any outcomes with a variety of distributions, binary is the most commonly used.

The assumption of latent class models is conditional or local independence, where the outcomes are independent conditional on the latent class. This assumes that subjects within a class are homogeneous. Where this does not apply, that is the true classes are heterogeneous, a consequence of this will be to increase the number of latent classes required in attempting to explain the heterogeneity, with possible consequent difficulty in interpretation. A solution is to incorporate random effects so that the outcomes are independent conditional on the latent class and random effect or effects.

There are a number of packages capable of fitting latent class models in R (R Core Team 2016). Two of these solely for fitting of latent class models are poLCA (Linzer and Lewis 2011) and BayesLCA (White and Murphy 2014). BayesLCA is particularly designed to perform Bayesian analyses, but also offers the choice of the EM algorithm and Variational Bayes, but has limited facilities for producing plots and summaries. poLCA is a more fully featured package which allows for polytomous outcomes and latent class regression, which are not available in randomLCA. The advantage of randomLCA over the other packages is that it will fit both standard latent class models and those incorporating random effects. This is important for use with diagnostic tests, as it allows for the variation of the test response between subjects, but may be also used to model heterogeneity in other applications, for example see Muthén (2006). Commercial software packages that also allow latent class with random effects are Mplus (Muthén and Muthén 2015) and Latent GOLD Syntax Module (Vermunt and Magidson 2013), both of which require the model to be defined using a symbolic language.

The purpose of this paper is to describe the randomLCA package. The remainder of the paper is organised as follows. Section 2 describes the models, starting with standard latent class and then continuing with the random effect extensions, including references allowing for further investigation. Section 3 describes three examples with explanation of how the features of the package may be used. Section 4 summarises the capabilities of the package and describes some areas in which the package could be extended.

2 Models

2.1 Latent class model

The basis of latent class analysis is that each subject is assumed to belong to one of a finite number of classes, with each class described by a set of parameters that define the distribution of outcomes or manifest variables for a subject, and is a form of finite mixture model (McLachlan and Peel 2000). Generally, the number of classes is unknown and must be determined from the data. For binary outcomes, the model is

\[ \begin{aligned} P\left(y_{i1}, y_{i2}, ..., y_{ik}\right|c_i=c)& = \prod_{j = 1}^k \pi_{cj}^{y_{ij}} \left(1-\pi_{cj}\right)^{1-y_{ij}}\\ \end{aligned} \]

where \(y_{ij}\) is the \(j\)th binary outcome for subject \(i\), \(\pi_{cj}\) is the probability of the \(j\)th outcome equal to 1 for a subject in class \(c\), \(k\) is the number of outcomes and \(c_i\) is the class corresponding to the \(i\)th subject. The marginal probability, obtained by summing over the classes, for each subject is \[ \begin{aligned} P\left(y_{i1}, y_{i2}, ..., y_{ik}\right)& = \sum_{c = 1}^C \eta_c\prod_{j = 1}^k \pi_{cj}^{y_{ij}} \left(1-\pi_{cj}\right)^{1-y_{ij}} \end{aligned} \] where \(\eta_c\) is the probability of a subject being in class \(c\) with \(\sum_{c = 1}^C{\eta_c} = 1\) and \(C\) is the number of classes. From this can be obtained the marginal likelihood.

A requirement for the estimates of the probabilities \(\pi_{cj}\) is that they be restricted to the interval zero to one, and that the \(\eta_c\) sum to one, something that did not always occur with the original methods used for analysis. A solution developed by formann1978 and formann1982 is to use a logistic (or alternatively probit) transformation, allowing unconstrained estimation of parameters but correctly restricting the probabilities to between zero and one. This can be obtained for the logistic using the following relations \(\pi_{cj} = {e^{a_{cj}}}/\left(1+e^{a_{cj}}\right)\) and \(\eta_c = {e^{\theta_c}}/{\sum_{\ell = 1}^C e^{\theta_\ell}}\). Hence we estimate the \(a_{cj}\) and \(\theta_c\), rather than \(\pi_{cj}\) and \(\eta_c\). Similar equations apply for the probit transformation.

When used as a classification algorithm the model does not simply return the most likely class for each subject but returns a probability of class membership, based on the observed outcomes. The posterior probability for a set of given observed outcomes can be obtained from Bayes theorem. Given the observed outcomes \(y_{i1}, y_{i2}, \cdots, y_{ik}\) then the probability that the subject is in class \(d\) is: \[ \begin{aligned} P\left(d| y_{i1}, y_{i2}, ..., y_{ik}\right)&=\frac{ \eta_d P\left(y_{i1}, y_{i2}, ..., y_{ik}\right|d)}{\sum_{c = 1}^C \eta_c P\left(y_{i1}, y_{i2}, ..., y_{ik}\right|c)}\\ & = \frac{ \eta_d\prod_{j = 1}^k \pi_{dj}^{y_{ij}} \left(1-\pi_{dj}\right)^{1-y_{ij}}}{\sum_{c = 1}^C \eta_c\prod_{j = 1}^k \pi_{cj}^{y_{ij}} \left(1-\pi_{cj}\right)^{1-y_{ij}}}. \end{aligned} \]

2.2 Latent class with random effect model

A major difficulty with latent class models is the requirement for local independence or equivalently homogeneity of the outcome probabilities within each class. When the classes are heterogeneous the assumption that the manifest outcomes are independent, conditional on the latent class, does not apply. Uebersax (1999) describes the problems associated with conditional dependence as “to add spurious latent classes that are not truly present at the taxonic level” and Vacek (1985) showed that ignoring the conditional dependence produced biased estimates in the context of diagnostic testing. Pickles and Angold (2003) discuss the classification of diseases in psychology as categories or by severity, arguing that ``most forms of psychopathology (indeed, most forms of pathology of any sort) manifest both continuous and discontinuous relationships with other phenomena’’, and provide examples of where this may occur.

A solution to the problem of heterogeneity was developed by Qu et al. (1996) combining a latent class model with a random effect to explain the heterogeneity. The probabilities are transformed to the probit scale and a normally distributed random effect added for each subject, before transforming back to probabilities. An alternative to the probit scale is the logit scale. A model for latent class incorporating a random effect \(~N(0, 1)\) is: \[ \begin{aligned} P\left(y_{i1}, y_{i2}, ..., y_{ik}|c_i=c, \lambda_i=\lambda\right)& = \prod_{j = 1}^k \pi_{icj}^{y_{ij}} \left(1-\pi_{icj}\right)^{1-y_{ij}} \end{aligned} \] where, either, if a probit scaling of the random effect \[ \begin{aligned} \pi_{icj} = \Phi^{-1}\left({a_{cj}+b_{cj}\lambda_i}\right) \end{aligned} \] or, if a logistic scaling \[ \begin{aligned} \pi_{icj} = \frac{\exp \left({a_{cj}+b_{cj}\lambda_i}\right)}{1+\exp \left({a_{cj}+b_{cj}\lambda_i}\right)}. \end{aligned} \] and \(a_{cj}\) determines the conditional class probability for a value of zero for the random effect, and \(b_{cj}\) scales the random effect, and is usually known as the loading or discriminant. The loadings will generally be constrained to be equal between classes, and either the same loading for each outcome (\(b_{cj} = b\)) or independent loading for each outcome (\(b_{cj} = b_j\)). The marginal likelihood is obtained by integrating over the random effect and summing over the latent classes, which is then maximised to obtain the parameter estimates. Posterior class probabilities can be obtained as for the standard latent class using Bayes theorem.

2.3 Two level latent class with random effect model

The previously described models for latent class can be considered to be for a single time point. Where the outcomes are observed at multiple time points then consideration must be made for the correlation between time points. A method described for this is the mixed latent Markov model (Langeheine and Pol 1990). This assumes that the population is a mixture of latent Markov models which consists of a Markov chain which is a latent class at each time point. An alternative model is that of Beath and Heller (2009) which extends the model to two levels with an additional random effect modelling the correlation between outcomes at each time point. This has some similarities to the model by Muthén and Shedden (1999) for longitudinal normally distributed data.

We now define \(y_{ijt}\) as the \(j\)th binary outcome for subject \(i\) at time \(t\), \(\pi_{cjt}\) is the probability of the \(j\)th outcome equal to 1 for a subject in class \(c\) at time \(t\), \(k\) is the number of outcomes, \(T\) is the number of time points and \(c\) is the class corresponding to the \(i\)th subject. An additional random effect \(\tau_{t} \sim N\left(0, 1\right)\) is incorporated to model the additional correlation between outcomes at a time point, and \(\ell_c\) scales this to the appropriate variance.

\[ \begin{aligned} P\left(y_{i11},\ldots,y_{ikT}|c_i=c,\lambda_i=\lambda,\tau_{i1}=\tau_1,\ldots,\tau_{iT}=\tau_T\right) &=&\prod_{t=1}^T\prod_{j=1}^k \pi_{icjt}^{y_{ijt}} \left(1-\pi_{icjt}\right)^{1-y_{ijt}} \end{aligned} \] where \[ \begin{aligned} \pi_{icjt}=\Phi^{-1}\left(a_{cjt}+b_{cj}\left(\lambda_i+\ell_c\tau_{it}\right)\right) \end{aligned} \] or \[ \begin{aligned} \pi_{icjt}=\frac{\exp\left(a_{cjt}+b_{cj}\left(\lambda_i+\ell_c\tau_{it}\right)\right)}{1+\exp\left(a_{cjt}+b_{cj}\left(\lambda_i+\ell_c\tau_{it}\right)\right)}. \end{aligned} \]

Again, we usually constrain \(b_{cj}=b_j\) and \(\ell_c=\ell\).

2.4 Model selection

An important aspect of using latent class and latent class with random effect models is the choice of the number of classes, and whether inclusion of a random effect is required. This presents two difficulties. The first involves a test for a parameter on the boundary of the parameter space. For the number of classes, this is the proportion in a class that is zero, and for the random effect, that the random effect variance is zero. As a consequence asymptotic likelihood theory does not hold so other methods must be used. One method is the bootstrapped likelihood ratio test as described by McLachlan (1987). However, a second difficulty is that when models may include a random effect then comparisons must be made between non-nested models. As an example we may wish to compare a latent class model with two or more classes to a latent class with random effect model with a single class. This is equivalent to an item response theory (IRT) model, either the Rasch model when the loadings are constant and a two parameter logistic otherwise (Bartholomew et al. 2002). This is an important choice as it determines if the underlying latent variable is categorical or continuous (Muthén 2006).

The usual method used is an information criterion (Lin and Dayton 1997) with the two main ones that are used being the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). Nylund and Muthén (2007) showed using simulation that BIC is superior to AIC for selection in latent class models, and this is the method most often used in applied publications. With BIC the penalty is greater than for AIC and dependent on the number of observations, so will select models with a smaller number of classes. An alternative version of AIC, AIC3 with a penalty of 3 was shown to have better performance with latent class models by Dias (2006). This will select models with a complexity between BIC and AIC, and this was found with the examples. No evaluation of the criteria has been performed for latent class with random effects models, and it has seen little use in applied publications. An important point is that with any model selection it is desirable to make use of existing information. So for example, if it is already known there is at least 2 classes then the model choice should be restricted to these. In the paper I have used BIC for model selection, but as no research has been performed on model selection for random effect latent class models the other information criteria are provided.

2.5 Identifiability

A difficulty that is sometimes encountered in fitting latent class models is lack of identifiability, which occurs when the value of the maximum likelihood occurs for more than one unique set of parameter estimates, that is there is not a single global maximum. A minimum condition is that the number of parameters is less than the number of patterns, but this is not always sufficient (McHugh 1956). This prevents the fitting of even a 2 class latent class model with only two outcomes, as this requires 5 parameters with only 4 patterns. A 2 class model is just identifiable when fitted to data on 3 outcomes, as there are 7 parameters and 8 possible patterns, however if a random effect is included then 4 outcomes are required for identifiability. A check is made by randomLCA to determine if the model is identifiable based on the number of parameters. As this is not always sufficient, special cases, for example the requirement of at least 5 outcomes to fit a 3 class latent class model, are also flagged as non-identified. A further check is provided by determining that the rank of the Hessian is not less than the number of estimated parameters (pages 150-1, Skrondal and Rabe-Hesketh 2004), however this may also occur if a parameter is on the boundary of the parameter space.

2.6 Computational methods

The standard latent class methods are fitted using an EM algorithm, with multiple starting values, switching to a quasi-Newton method when near convergence. The multiple starting values increase the probability that the algorithm converges to the global maximum rather than a local maximum, with only the EM algorithm run to completion when selecting starting values to reduce execution time. For the latent class with random effects models it is necessary to integrate over the random effects for each subject, for which it is necessary to use an approximation. One of the methods that has been used for this is Gauss-Hermite quadrature, where the integral is approximated by a weighted sum of the function evaluated at defined points, which works well when the function is approximately standard normal. However, for random effects models the likelihood is not, but the standardised likelihood is. Therefore, Gauss-Hermite quadrature is applied to the standardised likelihood, known as adaptive Gauss-Hermite quadrature (Liu and Pierce 1994). This has the advantage over standard Gauss-Hermite quadrature that integration is only performed in the region of the mode, reducing the number of quadrature points required. For the two level random effects model, adaptive Gauss-Hermite quadrature is again used. An orthogonal transform is applied to reduce the integration to two one-dimensional integrations and the method of moments is used to determine the location of the modes, as described in Rabe-Hesketh et al. (2005).

The random effects latent class models are fitted using a generalized EM algorithm (GEM) (p173, Little and Rubin 2002) with maximisation using a quasi-Newton method. At each expectation step the location of the modes for the adaptive quadrature are recalculated, and the maximisation step performed based on these locations. The maximisation is not run to convergence at each step, but terminated after a number of quasi-Newton steps. For all models, the EM or GEM algorithm switches to a quasi-Newton for all parameters when close to convergence. For the random effects models starting values are obtained from the model without random effects, and a search performed over possible starting values for the random effect variance.

A difficulty with latent class models is the calculation of standard errors when the parameter estimates are near the boundary of the parameter space. There are several options, one of which is to use Bayesian maximum a posteriori (MAP) estimation (Galindo Garre and Vermunt 2006). This is a form of Bayesian estimation in that a prior probability is placed on the parameters, however the posterior distribution is then maximised similar to maximum likelihood estimation. This is equivalent to the penalised likelihood described by Firth (1993) in which a penalty is placed on extreme values of the logit or probit scale outcome probabilities. The posterior distribution will consequently have a mode, whereas the likelihood may not. In randomLCA the prior distribution used is the Dirichlet as in Galindo Garre and Vermunt (2006) but with a much smaller default penalty, where the penalty argument is equal to the number of extra observations that are effectively added to each cell. A penalty of 0.01 is used, which reduces the risk of numerical problems, without greatly affecting the estimated probabilities. A sensible upper limit for the penalty is 0.5, which has been found to perform well by a number of authors, for example Rubin and Schenker (1987) in application to binary proportions, and may be used when necessary. This will also produce results similar to the Latent GOLD software default settings (Vermunt and Magidson 2013). Setting the penalty to zero will produce results identical to maximum likelihood.

3 Examples

3.1 Myocardial infarction example

This example demonstrates the fitting of data from Rindskopf and Rindskopf (1986), where latent class analysis is used to determine diagnostic classifications based on medical tests. Although this example is for medical data, the model is simply standard latent class so the methods can be applied to data from other areas, for example psychology and sociology. The maximum number of classes that can be fitted is limited to 2 due to identifiability, so we fit models for 1 and 2 classes, assuming that the outcome probabilities are homogenous in each class.

The fitting function for randomLCA is the randomLCA function which fits both the standard and random effects models. The command, where only the patterns parameter is required, is:

randomLCA(patterns, freq = NULL, nclass = 2, calcSE = TRUE, notrials = 20,
          random = FALSE, byclass = FALSE, quadpoints = 21, constload = TRUE,
          blocksize = dim(patterns)[2], level2 = FALSE, probit = FALSE,
          level2size = blocksize, qniterations = 5, penalty = 0.01, EMtol = 1.0e-9,
          verbose = FALSE, seed = as.integer(runif(1, 0, .Machine$integer.max)),
          cores =  max(detectCores(logical = FALSE) \%/\% 2, 1))

For a standard latent class model the parameters of interest are

The remaining arguments will be considered when discussing the relevant models or may be obtained from the package documentation. Fitting the standard latent class models for one and two classes (note that a three class model is not identifiable) requires only the basic arguments:

data(myocardial)

myocardial.lca1 <- randomLCA(myocardial[, 1:4], freq = myocardial$freq,
                             nclass = 1)
myocardial.lca2 <- randomLCA(myocardial[, 1:4], freq = myocardial$freq,
                             nclass = 2)

The BIC values may be extracted from the fitted objects, and formed into a data frame:

myocardial.bic <- data.frame(classes = 1:2, bic = c(BIC(myocardial.lca1),
                                                    BIC(myocardial.lca2)))
print(myocardial.bic, row.names = FALSE)
#>  classes      bic
#>        1 524.7441
#>        2 402.2951

Using BIC as a selection method, this selects the 2 class model, indicating a breakdown into diseased and non-diseased, which it is assumed represent those with and without myocardial infarction, although the true nature of classes is always debatable. A characteristic of this data is that a single class random effect model has a lower BIC than the 2 class standard latent class, so we need to assume that there are at least two classes, or that the underlying latent variable is categorical rather than continuous.

An alternative is to use the parametric bootstrap (McLachlan 1987) to determine the number of classes, and this is easily performed using the simulate function to generate the samples under the null hypothesis and then refit to refit both models. The simulate function returns a list of data frames simulated from the specified model, and each of these can then be refitted using the specified null and alternative model, as shown in the following code.

nsims <- 999
obslrt <- 2 * (logLik(myocardial.lca2) - logLik(myocardial.lca1))
thesims <- simulate(myocardial.lca1, nsim = nsims)
simlrt <- as.vector(lapply(thesims, function(x) {
  submodel <- refit(myocardial.lca1, newpatterns = x)
  fullmodel <- refit(myocardial.lca2, newpatterns = x)
  return(2 * (logLik(fullmodel) - logLik(submodel)))
}))

A \(p\)-value is obtained by comparing the observed likelihood ratio test statistic to the simulated. This can be performed in a number of ways, of which the most commonly used is described by (page 148, Davison and Hinkley 1997).

print((sum(simlrt >= obslrt) + 1)/(nsims +  1))
#> [1] 0.001

Showing again the clear evidence in favour of the 2 class model.

Summary may be used to display the fitted results:

summary(myocardial.lca2)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik
#>         2 379.4054 402.2951 388.4054 -180.7027 -180.7829
#> Class probabilities 
#> Class  1 Class  2 
#>   0.5422   0.4578 
#> Outcome probabilities 
#>          Q.wave History    LDH    CPK
#> Class  1 0.0001  0.1951 0.0270 0.1956
#> Class  2 0.7668  0.7914 0.8279 0.9999

From this it is clear that Class 2 is the diseased class and Class 1 the non-diseased, where the disease is myocardial infarction, based on the higher outcome probabilities. Outcome probabilities are plotted using the plot function, and shown in Figure 1. Note that plot is based on xyplot so the additional graphical arguments must be the appropriate lattice ones.

plot(myocardial.lca2, type = "b", pch = 1:2, xlab = "Test",
       ylab = "Outcome Probability",
       scales = list(x = list(at = 1:4, labels = names(myocardial)[1:4])),
       key = list(corner = c(0.05, .95), border = TRUE, cex = 1.2,
                  text = list(c("Class 1", "Class 2")),
                  col = trellis.par.get()$superpose.symbol$col[1:2],
                  points = list(pch = 1:2)))
Outcome probabilities for 2 class latent class model for myocardial infarction data.

Figure 1: Outcome probabilities for 2 class latent class model for myocardial infarction data.

Individual results may be obtained from summary, for example the outcome probabilities but may also be obtained using the outcomeProbs function, which will also give the 95% confidence intervals.

outcomeProbs(myocardial.lca2)
#> Class  1 
#>            Outcome p        2.5 %    97.5 %
#> Q.wave  6.882258e-05 6.549913e-22 1.0000000
#> History 1.951154e-01 1.017696e-01 0.3415263
#> LDH     2.697838e-02 3.831996e-03 0.1665593
#> CPK     1.956158e-01 9.784094e-02 0.3528805
#> Class  2 
#>         Outcome p        2.5 %    97.5 %
#> Q.wave  0.7668256 5.919036e-01 0.8817497
#> History 0.7914045 6.388627e-01 0.8905522
#> LDH     0.8278976 6.552020e-01 0.9241148
#> CPK     0.9999367 1.464933e-13 1.0000000

For Q.wave in Class 1 there appears to be a problem with the standard errors, as assumptions about the normal approximation to the likelihood do not apply close to the boundary. Using the parametric bootstrap with boot = TRUE will produce improved results or alternatively, the value of the penalty argument could be increased.

outcomeProbs(myocardial.lca2, boot = TRUE)
#> Class  1 
#>            Outcome p        2.5 %       97.5 %
#> Q.wave  6.882867e-05 4.400879e-05 0.0001289999
#> History 1.951155e-01 8.419667e-02 0.3127504331
#> LDH     2.697840e-02 5.138272e-05 0.0877891115
#> CPK     1.956159e-01 7.967963e-02 0.3214772816
#> Class  2 
#>         Outcome p     2.5 %    97.5 %
#> Q.wave  0.7668256 0.6201619 0.9068112
#> History 0.7914045 0.6614429 0.9126307
#> LDH     0.8278976 0.7014661 0.9505552
#> CPK     0.9999367 0.9804981 0.9999489

The outcome probabilities give some interesting information. For example, in Class 1, those without myocardial infarction, will have absence of Q.wave but in those with myocardial infarction it will only be present in 76.7%. The class probabilities can be obtained as classProbs(myocardial.lca2) of 0.54 and 0.46 for Class 1 and 2 respectively.

One aspect of latent class is that no subject is uniquely allocated to a given class, although in some cases a subject may have an extremely high probability of being in a given class. The posterior class probabilities can be obtained as

print(postClassProbs(myocardial.lca2), row.names = FALSE)
#>  Q.wave History LDH CPK Freq      Class 1      Class 2
#>       1       1   1   1   24 1.670757e-07 9.999998e-01
#>       0       1   1   1    5 7.919814e-03 9.920802e-01
#>       1       0   1   1    4 2.614857e-06 9.999974e-01
#>       0       0   1   1    3 1.110642e-01 8.889358e-01
#>       1       1   0   1    3 2.898659e-05 9.999710e-01
#>       0       1   0   1    5 5.807210e-01 4.192790e-01
#>       1       0   0   1    2 4.534697e-04 9.995465e-01
#>       0       0   0   1    7 9.559025e-01 4.409746e-02
#>       0       0   1   0    1 9.998768e-01 1.232123e-04
#>       0       1   0   0    7 9.999889e-01 1.111584e-05
#>       0       0   0   0   33 9.999993e-01 7.102497e-07

This shows subjects with 3 or 4 positive tests to be strongly classified as having myocardial infarction, and even some with 2 positive tests are well classified. Having only one positive test makes it unlikely that it is myocardial infarction.

3.2 Dentistry example

An important area of application of latent class and random effects latent class is the development or comparison of diagnostic testing methods where there is no gold standard test. A gold standard test is one that is the best available and can usually be assumed to be close to perfect, but usually being more expensive or difficult to perform (Kraemer 1992). Given a gold standard it is easy to construct new tests or compare existing tests, as we know the true disease status of each subject. Latent class methods allow the construction of tests based on the assumption that subjects fall in either two or more classes, with diseased or non-diseased as a minimum, except that the classes can only be inferred from the observed test results. This has the consequence that the status of the subjects is not known exactly, which reduces the accuracy and relies upon the assumptions made about the test result distribution.

There are other methods that have been proposed, the major of which are discrepant resolution and composite reference (Chapter 7, Pepe 2003), both of which have disadvantages and advantages compared to latent class methods. The advantage of the latent class method is its statistical basis, however it has a disadvantage of dependency on assumptions about the diagnostic tests, especially the assumption of either conditional independence or normally distributed heterogeneity.

The further arguments to the randomLCA function required for a random effects model are:

random Specifies whether a random effect should be included. byclass Allow loadings for the random effect to vary by class. quadpoints Number of quadrature points for the adaptive quadrature. These specify how accurate the numerical approximation to the marginal likelihood is, and should be increased until there is negligible improvement in model fit. constload The same loading is used for all outcomes when using a random effects model. blocksize If the outcome loadings are broken into blocks what is the block size? This allows the number of \(b_cj\) parameters associated with the \(\lambda_i\) to be reduced when fitting models with a large number of outcomes and as tructure to the outcomes by placing constraints on the \(b_cj\). This will be demonstrated in the symptoms example. *normalfont Fit probit model rather than logistic for relationship between parameters and outcome probabilities. This is the relationship typically used in some disciplines.

This example shows the fitting of the dentistry data from Qu et al. (1996). The data consists of the results of five dentists evaluating x-rays for presence or absence of caries. For consistency with the original paper I have also set probit = TRUE to give the probit link. Fitting first the three possible models for one class:

dentistry.lca1 <- randomLCA(dentistry[, 1:5], freq = dentistry$freq, nclass = 1)
dentistry.lca1random <- randomLCA(dentistry[, 1:5], freq = dentistry$freq,
  nclass = 1, random = TRUE, probit = TRUE)
dentistry.lca1random2 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,
  nclass = 1, random = TRUE, probit = TRUE, constload = FALSE)

This can then be repeated for 2 to 4 classes, and using BIC the BIC extracted for each model, and then formed into a data frame to summarise. Note that we cannot use a parametric bootstrap based likelihood ratio test to compare the standard to random effects latent class as the models are not nested. Here the number of quadrature points will need to be increased for some models to allow convergence.

dentistry.lca2 <- randomLCA(dentistry[, 1:5], freq = dentistry$freq, nclass=2)
dentistry.lca3 <- randomLCA(dentistry[,1:5], freq= dentistry$freq, nclass = 3)
dentistry.lca4 <- randomLCA(dentistry[, 1:5], freq = dentistry$freq,nclass=4)
dentistry.lca2random <- randomLCA(dentistry[,1:5],freq= dentistry$freq,
  nclass = 2, random = TRUE, probit = TRUE)
dentistry.lca3random <- randomLCA(dentistry[, 1:5], freq =
  dentistry$freq,nclass=3,random=TRUE,probit=TRUE)
dentistry.lca4random <- randomLCA(dentistry[,1:5],freq= dentistry$freq,
  nclass = 4, random = TRUE, quadpoints = 31, probit = TRUE)
dentistry.lca2random2 <- randomLCA(dentistry[, 1:5], freq =
  dentistry$freq,nclass=2,random=TRUE,quadpoints=41,probit=TRUE,constload=FALSE)
  dentistry.lca3random2 <- randomLCA(dentistry[,1:5],freq= dentistry$freq,
  nclass = 3, random = TRUE, quadpoints = 41, probit = TRUE,
  constload = FALSE)

We can display the BIC values in a table, with the first column for standard latent class, second for random effects with a constant loading for each dentist, and the third with loading varying by dentist. Note that it is not possible to fit a 4 class random effects model with individual loadings for each dentist due to non-identifiability. As for the previous example we can form the BIC values into a table, where bic, bicrandom and bicrandom2 are the BIC values from the standard, random effects with constant loading and random effects with non-constant loading latent class models.

bic.data <- data.frame(classes = 1:4, bic = c(BIC(dentistry.lca1),
  BIC(dentistry.lca2), BIC(dentistry.lca3), BIC(dentistry.lca4)),
  bicrandom = c(BIC(dentistry.lca1random), BIC(dentistry.lca2random),
  BIC(dentistry.lca3random), BIC(dentistry.lca4random)),
  bicrandom2 = c(BIC(dentistry.lca1random2), BIC(dentistry.lca2random2),
  BIC(dentistry.lca3random2), NA))
kable(bic.data, format = "html")
Table 1: BIC for Dentistry Models
classes bic bicrandom bicrandom2
1 17531.13 14974.79 14938.27
2 15021.64 14944.69 14949.37
3 14962.89 14963.54 14992.33
4 15000.03 15007.19 NA

For the standard latent class models the minimum BIC of 1.49629^{4} is obtained for the 3 class model. With addition of the random effect with constant loading, minimum BIC is obtained with a 2 class model with a decrease from the latent class model to 1.49447^{4}. Allowing the loadings to vary by dentist (2LCR model obtained by Qu et al. (1996)) minimum BIC of 1.49383^{4} was obtained using a single class model, equivalent to a single factor Item Response Theory (IRT) model (Chapter 7, Bartholomew et al. 2002). This assumes that rather than subjects being grouped into classes they simply have different levels of an underlying latent variable, possibly severity. In the absence of any assumptions about the appropriate model this would be the model to be used, and we could conclude that severity of caries was on a continuous scales with each dentist having different thresholds for determining the presence or absence. In the paper by Qu et al. (1996) the assumption is made that the underlying latent variable is categorical, that is that there are two distinct types of subjects, and so the 2 class with random effect with constant loading will be used.

Summary may be used to display the fitted results:

summary(dentistry.lca2random)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik   Link
#>         2 14869.56 14944.69 14881.56 -7422.782 -7422.848 Probit
#> Class probabilities 
#> Class  1 Class  2 
#>    0.179    0.821 
#> Conditional outcome probabilities 
#>              V1     V2     V3     V4     V5
#> Class  1 0.3804 0.6770 0.6508 0.3999 0.9033
#> Class  2 0.0057 0.0836 0.0059 0.0258 0.2964
#> Marginal Outcome Probabilities 
#>              V1     V2     V3     V4     V5
#> Class  1 0.4017 0.6464 0.6243 0.4179 0.8561
#> Class  2 0.0192 0.1294 0.0199 0.0558 0.3310
#> Loadings 
#> 0.704417

Clearly, based on the lower outcome probabilities, Class 1 is the non-diseased and Class 2 the diseased.

For latent class models with random effects there are two additional arguments to plot

graphtype Type of graph, either “marginal” or “conditional”. For marginal the outcome probabilities integrated over the random effect are plotted, and for conditional they are plotted conditional on the random effect, with zero the default.

conditionalp For a conditional graph the percentile corresponding to the random effect at which the outcome probability is to be calculated. Fifty percent is the default, corresponding to a random effect value of zero.

The marginal outcome probabilities, obtained by integrating over the random effect can be plotted, as in Figure 2. The marginal outcome probabilities reflect the average probability for any subject in the class having a positive outcome. This differs from the conditional outcome probabilities which are for a subject with zero random effect, and thus represent a typical subject.

plot(dentistry.lca2random, graphtype = "marginal", type = "b", pch = 1:2,
  xlab = "Dentist", ylab = "Marginal Outcome Probability",
  key = list(corner = c(0.05, .95), border = TRUE, cex = 1.2,
  text = list(c("Class 1", "Class 2")),
  col = trellis.par.get()$superpose.symbol$col[1:2],
  points = list(pch = 1:2)))
Marginal outcome probabilities for 2 class latent class with random effect (2LCR) model for dentistry data.

Figure 2: Marginal outcome probabilities for 2 class latent class with random effect (2LCR) model for dentistry data.

Adding the boot = TRUE parameter to the outcomeProbs function will obtain bootstrap confidence intervals. Differences from the Qu et al paper result from their use of an individual loading for each dentist when calculating Table 6.

We can demonstrate the effect of the random effect by plotting the outcome probabilities by varying percentiles of the random effect, using the following code, which will place each percentile in a different panel.

plot(dentistry.lca2random, graphtype = "conditional", type = "b",
  pch = 1:2, conditionalp = c(0.025, 0.5, 0.975),
  scales = list(alternating = FALSE, x = list(cex = 0.8)),
  xlab = "Dentist", ylab = "Conditional Outcome Probability",
  key = list(corner = c(0.05, .9), border = TRUE, cex = 1.2,
  text = list(c("Class 1", "Class 2")),
  col = trellis.par.get()$superpose.symbol$col[1:2],
  points = TRUE))
Conditional outcome probabilities for 2 class latent class with random effect (2LCR) model for dentistry data.

Figure 3: Conditional outcome probabilities for 2 class latent class with random effect (2LCR) model for dentistry data.

If it is desired to have the plots on a single graph, it requires using the calcCondProb function which returns a data.frame containing the conditional probabilities conditional on the random effect for each class and outcome. We use the 2.5th and 97.5th percentiles, and superpose the plots on the same graph, as shown in Figure ??. As an alternative each class could be placed in a separate panel.

probs <- calcCondProb(dentistry.lca2random, conditionalp =c(0.025, 0.5, 0.975))

Two important concepts in diagnostic testing are sensitivity and specificity. Sensitivity is the probability of obtaining a positive result given that the true state is positive, and specificity is the probability of a negative result given that the true state is negative. Calculation of sensitivity and specificity is shown in the following code, where the number of quadrature points have been increased to ensure convergence for all simulated datasets but may also be obtained by increasing the penalty. The default for evaluation of outcome probabilities for random effect models is marginal, and this is appropriate for obtaining sensitivity and specificity.

dentistry.lca2random <- randomLCA(dentistry[, 1:5], freq = dentistry$freq,
  nclass = 2, random = TRUE, quadpoints = 71, probit = TRUE)
probs <- outcomeProbs(dentistry.lca2random, boot = TRUE)

It is necessary to determine which is the class with higher outcome probabilities, as it is the diseased class. The variable diseased gives the number of the diseased class, and notdiseased for the non-diseased, and is based on the outcome probabilities being lower for the non-diseased class.

diseased <- ifelse(probs[[1]]$Outcome[1] < probs[[2]]$Outcome[1], 2, 1)
notdiseased <- 3 - diseased
sens <- apply(probs[[diseased]], 1, function(x)
  sprintf("%3.2f (%3.2f, %3.2f)", x[1], x[2], x[3]))
spec <- apply(probs[[notdiseased]], 1, function(x)
  sprintf("%3.2f (%3.2f, %3.2f)", 1 - x[1], 1 - x[3], 1 - x[2]))

stable <- data.frame(sens, spec)
names(stable) <- c("Sensitivity", "Specificity")
kable(stable, format = "html")
Table 2: Sensitivity and Specificity for Dentistry data
Sensitivity Specificity
V1 0.40 (0.33, 0.47) 0.98 (0.97, 0.99)
V2 0.65 (0.57, 0.72) 0.87 (0.85, 0.89)
V3 0.62 (0.53, 0.78) 0.98 (0.96, 1.00)
V4 0.42 (0.35, 0.49) 0.94 (0.93, 0.96)
V5 0.86 (0.79, 0.91) 0.67 (0.64, 0.69)

The true and false positive rates can be calculated from the outcome probabilities, similar to sensitivity and sensitivity, and plotted for each dentist, and are shown in Figure 4.

rates <- data.frame(tpr=probs[[diseased]][, 1],
  fpr=probs[[notdiseased]][, 1])
plot(tpr~fpr, type= "p",
  xlab= "False Positive Rate\n(1-Specificity)",
  ylab= "True Positive Rate (Sensitivity)",
  xlim=c(0.0, 0.35), ylim=c(0.35, 0.9), data=rates)
text(rates$fpr, rates$tpr,labels=1:length(rates$fpr), pos=4)
True and false positive rates by dentist.

Figure 4: True and false positive rates by dentist.

This gives a better explanation. It appears that the difference between dentists is mainly related to the threshold for what they classify as diseased. Dentist 5 is more likely to correctly identify teeth as diseased but at the expense of being more likely to identify non-diseased teeth as diseased. Note that this is different from an ROC curve where the same data is used but the test threshold is adjusted. Here the dentists may choose different thresholds but may also have different levels of performance.

Posterior class probabilities may again be obtained with postClassProbs

print(postClassProbs(dentistry.lca2random), row.names = FALSE)
#>  V1 V2 V3 V4 V5 Freq    Class 1    Class 2
#>   0  0  0  0  0 1880 0.01510022 0.98489978
#>   0  0  0  0  1  789 0.06917821 0.93082179
#>   0  0  0  1  0   43 0.06932842 0.93067158
#>   0  0  0  1  1   75 0.20226869 0.79773131
#>   0  0  1  0  0   23 0.45787606 0.54212394
#>   0  0  1  0  1   63 0.70501594 0.29498406
#>   0  0  1  1  0    8 0.65707871 0.34292129
#>   0  0  1  1  1   22 0.84127347 0.15872653
#>   0  1  0  0  0  188 0.07570102 0.92429898
#>   0  1  0  0  1  191 0.22845928 0.77154072
#>   0  1  0  1  0   17 0.19730870 0.80269130
#>   0  1  0  1  1   67 0.44250810 0.55749190
#>   0  1  1  0  0   15 0.69765253 0.30234747
#>   0  1  1  0  1   85 0.86945157 0.13054843
#>   0  1  1  1  0    8 0.82628376 0.17371624
#>   0  1  1  1  1   56 0.93698413 0.06301587
#>   1  0  0  0  0   22 0.20830139 0.79169861
#>   1  0  0  0  1   26 0.43828314 0.56171686
#>   1  0  0  1  0    6 0.38511918 0.61488082
#>   1  0  0  1  1   14 0.63614568 0.36385432
#>   1  0  1  0  0    1 0.84992051 0.15007949
#>   1  0  1  0  1   20 0.93487342 0.06512658
#>   1  0  1  1  0    2 0.91123326 0.08876674
#>   1  0  1  1  1   17 0.96597650 0.03402350
#>   1  1  0  0  0    2 0.42926621 0.57073379
#>   1  1  0  0  1   20 0.68678171 0.31321829
#>   1  1  0  1  0    6 0.61133970 0.38866030
#>   1  1  0  1  1   27 0.82621016 0.17378984
#>   1  1  1  0  0    3 0.92818913 0.07181087
#>   1  1  1  0  1   72 0.97420229 0.02579771
#>   1  1  1  1  0    1 0.95985001 0.04014999
#>   1  1  1  1  1  100 0.98854606 0.01145394

Clearly, as the number of dentists identifying the subject as diseased increases, the posterior probability of being diseased increases, until it is almost one when all dentists identify the subject as diseased. The observed and fitted values may be obtained using the fitted method which returns a data frame containing them. Again, differences from the Qu et al paper result from a model with different loading for each dentist. We can obtain the fitted values for the two models as follows:

dentistry.lca2.fitted <- fitted(dentistry.lca2)
dentistry.lca2random.fitted <- fitted(dentistry.lca2random)
dentistry.fitted <- merge(dentistry.lca2.fitted,
dentistry.lca2random.fitted, by = names(dentistry.lca2.fitted)[1:6])
names(dentistry.fitted)[6:8] <- c("Obs", "Exp 2LC", "Exp 2LCR")
print(dentistry.fitted, row.names = FALSE)
#>  V1 V2 V3 V4 V5  Obs     Exp 2LC    Exp 2LCR
#>   0  0  0  0  0 1880 1836.269026 1882.192314
#>   0  0  0  0  1  789  830.354551  779.797946
#>   0  0  0  1  0   43   61.935085   56.084218
#>   0  0  0  1  1   75   49.638377   72.333935
#>   0  0  1  0  0   23   28.631092   25.843402
#>   0  0  1  0  1   63   47.477140   60.365790
#>   0  0  1  1  0    8    4.035499    4.693958
#>   0  0  1  1  1   22   35.146035   25.089527
#>   0  1  0  0  0  188  213.894932  176.134369
#>   0  1  0  0  1  191  152.205162  209.645061
#>   0  1  0  1  0   17   12.146925   17.583716
#>   0  1  0  1  1   67   61.010487   53.823638
#>   0  1  1  0  0   15   11.208779   14.481080
#>   0  1  1  0  1   85   91.568856   78.991523
#>   0  1  1  1  0    8    8.068104    5.582298
#>   0  1  1  1  1   56   86.407020   67.144795
#>   1  0  0  0  0   22   21.211304   16.998029
#>   1  0  0  0  1   26   25.170898   30.581732
#>   1  0  0  1  0    6    2.100454    2.526789
#>   1  0  0  1  1   14   16.081237   10.598156
#>   1  0  1  0  0    1    2.541628    3.323290
#>   1  0  1  0  1   20   24.707588   20.750499
#>   1  0  1  1  0    2    2.180103    1.437973
#>   1  0  1  1  1   17   23.519057   17.796229
#>   1  1  0  0  0    2    6.023203    7.401449
#>   1  1  0  0  1   20   42.001102   31.873369
#>   1  1  0  1  0    6    3.694996    2.415691
#>   1  1  0  1  1   27   39.260319   23.634112
#>   1  1  1  0  0    3    5.667911    4.566133
#>   1  1  1  0  1   72   61.064430   59.378207
#>   1  1  1  1  0    1    5.392063    3.466836
#>   1  1  1  1  1  100   58.386637  102.463936

It can be seen how the random effects model more accurately models the data, with fitted values closer to the observed data.

3.3 Symptoms example

This comprises data on the presence or absence of respiratory and allergy symptoms in the Childhood Asthma Prevention Study (CAPS) (Mihrshahi et al. 2001) and was used as the example in Beath and Heller (2009). The symptoms of night cough, wheeze, itchy rash and flexural dermatitis since the previous visit were recorded at one month, then quarterly for the first year and then twice yearly until age two years. For analysis these are aggregated for each six month period to avoid numerical problems associated with very small probabilities. The aim of the analysis is to identify the number of classes of subjects based on their respiratory and allergy symptoms combined. As well as allowing for the classes to be defined by different levels of the symptoms it will also allow for changes over time.

For randomLCA the data is required in wide format with the four outcomes repeated in order for the number of periods. While it is not necessary, each outcome is suffixed with a period and the corresponding time identifier. This makes interpretation of the results easier and also will be used in labelling of the graphs. The structure of the data is:

names(symptoms)
#>  [1] "Nightcough.13" "Wheeze.13"     "Itchyrash.13"  "FlexDerma.13" 
#>  [5] "Nightcough.45" "Wheeze.45"     "Itchyrash.45"  "FlexDerma.45" 
#>  [9] "Nightcough.6"  "Wheeze.6"      "Itchyrash.6"   "FlexDerma.6"  
#> [13] "Nightcough.7"  "Wheeze.7"      "Itchyrash.7"   "FlexDerma.7"  
#> [17] "Freq"

The following two additional arguments to randomLCA are used to define the two level random effects model.

level2 Fit 2 level random effects model.

level2size Size of level 2 blocks for fitting 2 level models.

The first model fitted is a standard latent class, to allow for no subject or period effect.

symptoms.lca2 <- randomLCA(symptoms[,1:16],
  freq=symptoms$Freq,nclass=2)

A variation of the random effects latent class model can be fitted allowing the loadings (\(b_cj\)parameters$) for the outcomes to be repeated, that is wheeze at the different time points, for example, will always have the same loading, using the blocksize argument. This is equivalent to the 2 level model with the time-dependent random variable having a variance of zero, that is no time-dependent effect. The outcomes have been set to have non-constant loading, although in practice models for constant loading would also be fitted.

symptoms.lca2random <- randomLCA(symptoms[,1:16], freq=symptoms$Freq,
  random=TRUE, nclass=2, blocksize=4, constload=FALSE)

The two level models are specified through the level2 argument and the number of outcomes at each time through the level2size argument. For these models the penalty is increased to 0.1 to reduce the execution time, but for a two class model will be about an hour and for a three class model about two hours.

symptoms.lca2random2 <- randomLCA(symptoms[,1:16], freq=symptoms$Freq,
  random=TRUE, level2=TRUE, nclass=2, level2size=4, constload=FALSE,
  penalty=0.1)
symptoms.lca1 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,nclass=1)
summary(symptoms.lca1)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik
#>         1 10774.49 10844.19 10790.49 -5371.244 -5371.372
#> Class probabilities 
#> Class  1 
#>        1 
#> Outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5536    0.4007        0.406       0.2953
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.6252    0.3784        0.382       0.2739
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.4557   0.2459      0.2713       0.179       0.3952
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1893      0.2574      0.1452

symptoms.lca2 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,nclass=2)
summary(symptoms.lca2)
#>   Classes      AIC      BIC     AIC3   logLik penlogLik
#>         2 9980.059 10123.81 10013.06 -4957.03 -4957.172
#> Class probabilities 
#> Class  1 Class  2 
#>   0.6534   0.3466 
#> Outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5272    0.3763       0.1893       0.1008
#> Class  2        0.6032    0.4466       0.8135       0.6609
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.6196    0.3477       0.1598       0.0753
#> Class  2        0.6359    0.4366       0.8028       0.6501
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.4207   0.2160      0.0706      0.0492       0.3284
#> Class  2       0.5208   0.3018      0.6451      0.4210       0.5206
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1316      0.1034      0.0493
#> Class  2   0.2977      0.5462      0.3251

symptoms.lca3 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,nclass=3)
summary(symptoms.lca3)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik
#>         3 9743.034 9960.839 9793.034 -4821.517 -4821.686
#> Class probabilities 
#> Class  1 Class  2 Class  3 
#>   0.2900   0.3874   0.3225 
#> Outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.7744    0.5915       0.2295       0.1293
#> Class  2        0.3471    0.2238       0.1769       0.0968
#> Class  3        0.6038    0.4422       0.8393       0.6824
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.9014    0.6353       0.2137       0.0748
#> Class  2        0.4122    0.1493       0.1508       0.0968
#> Class  3        0.6348    0.4246       0.8124       0.6666
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.6380   0.4260      0.0388      0.0196       0.5993
#> Class  2       0.2703   0.0705      0.1124      0.0803       0.1407
#> Class  3       0.5178   0.2980      0.6662      0.4376       0.5209
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.3278      0.0777      0.0403
#> Class  2   0.0000      0.1303      0.0565
#> Class  3   0.2947      0.5702      0.3454

symptoms.lca4 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,nclass=4)
summary(symptoms.lca4)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik
#>         4 9618.883 9910.742 9685.883 -4742.441 -4742.619
#> Class probabilities 
#> Class  1 Class  2 Class  3 Class  4 
#>   0.1620   0.2115   0.3559   0.2706 
#> Outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.4443    0.2579       0.9501       0.8971
#> Class  2        0.6858    0.5412       0.7834       0.5721
#> Class  3        0.3674    0.2387       0.0818       0.0000
#> Class  4        0.7623    0.5913       0.2088       0.1034
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.4875    0.2188       0.6063       0.4802
#> Class  2        0.6884    0.5113       0.7970       0.6580
#> Class  3        0.4362    0.1621       0.1492       0.0941
#> Class  4        0.9085    0.6564       0.2304       0.0872
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.3524   0.1093      0.2834      0.1339       0.4054
#> Class  2       0.5977   0.4081      0.8189      0.5650       0.5780
#> Class  3       0.2787   0.0686      0.1244      0.0949       0.1178
#> Class  4       0.6429   0.4377      0.0232      0.0112       0.6135
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1036      0.1523      0.0555
#> Class  2   0.3843      0.7685      0.5003
#> Class  3   0.0000      0.1413      0.0562
#> Class  4   0.3396      0.0720      0.0379

symptoms.lca5 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,nclass=5)
summary(symptoms.lca5)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik
#>         5 9553.011 9918.924 9637.011 -4692.506  -4692.72
#> Class probabilities 
#> Class  1 Class  2 Class  3 Class  4 Class  5 
#>   0.1085   0.3377   0.2557   0.1030   0.1951 
#> Outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.4019    0.2164       1.0000       0.9999
#> Class  2        0.3829    0.2454       0.0897       0.0064
#> Class  3        0.7831    0.6160       0.2508       0.1410
#> Class  4        0.4478    0.2633       0.4622       0.3177
#> Class  5        0.6913    0.5657       0.7943       0.5904
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.4494    0.0990       0.5370       0.3378
#> Class  2        0.4465    0.1726       0.0758       0.0149
#> Class  3        0.8929    0.6436       0.1687       0.0001
#> Class  4        0.6178    0.4597       0.9526       1.0000
#> Class  5        0.6872    0.5008       0.8041       0.6608
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.3131   0.0667      0.3482      0.2066       0.3660
#> Class  2       0.2868   0.0734      0.1064      0.0857       0.1288
#> Class  3       0.6232   0.4461      0.0458      0.0212       0.6311
#> Class  4       0.4772   0.1942      0.2563      0.0896       0.3730
#> Class  5       0.6004   0.4149      0.8096      0.5749       0.5777
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.0897      0.2303      0.0786
#> Class  2   0.0000      0.1489      0.0635
#> Class  3   0.3630      0.0785      0.0398
#> Class  4   0.0968      0.0000      0.0000
#> Class  5   0.3962      0.8267      0.5365

symptoms.lca1random <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,nclass=1,blocksize=4,constload=FALSE)
summary(symptoms.lca1random)
#>   Classes      AIC     BIC     AIC3    logLik penlogLik  Link
#>         1 9900.388 9987.51 9920.388 -4930.194 -4930.336 Logit
#> Class probabilities 
#> Class  1 
#>        1 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5547    0.3977       0.3456       0.2022
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.6274    0.3747       0.3097       0.1776
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.4546   0.2393      0.1636      0.0865       0.3932
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1828      0.1496      0.0626
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5537    0.4009       0.4046       0.2947
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.6254    0.3785       0.3809       0.2734
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.4554   0.2456       0.269      0.1781        0.395
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1892       0.256      0.1453
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.2654 0.3717    2.0099    1.8752

symptoms.lca2random <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,nclass=2,blocksize=4,constload=FALSE)
summary(symptoms.lca2random)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik  Link
#>         2 9608.567 9769.743 9645.567 -4767.284 -4767.437 Logit
#> Class probabilities 
#> Class  1 Class  2 
#>   0.4431   0.5569 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.7997    0.6151       0.2763       0.1250
#> Class  2        0.3598    0.2267       0.4035       0.2703
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.8788    0.6631       0.2908       0.1294
#> Class  2        0.4266    0.1508       0.3245       0.2158
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.6435   0.4398      0.1350      0.0573       0.5842
#> Class  2       0.3065   0.0850      0.1861      0.1093       0.2455
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.3394      0.1349      0.0718
#> Class  2   0.0625      0.1598      0.0527
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.7962    0.6103       0.3586       0.2252
#> Class  2        0.3621    0.2353       0.4415       0.3496
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.8758    0.6567       0.3687       0.2298
#> Class  2        0.4279    0.1591       0.3913       0.3080
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.6411   0.4424      0.2428      0.1399       0.5827
#> Class  2       0.3096   0.0911      0.2896      0.2080       0.2490
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.3458      0.2428      0.1613
#> Class  2   0.0674      0.2665      0.1325
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.2742  0.436    2.0264     1.914

symptoms.lca3random <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,nclass=3,blocksize=4,constload=FALSE)
summary(symptoms.lca3random)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik  Link
#>         3 9465.032 9700.262 9519.032 -4678.516 -4678.701 Logit
#> Class probabilities 
#> Class  1 Class  2 Class  3 
#>   0.3452   0.2835   0.3713 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.3715    0.2409       0.0106       0.0000
#> Class  2        0.4550    0.2859       0.9812       0.8942
#> Class  3        0.7993    0.6355       0.2502       0.0891
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.4424    0.1550       0.1033       0.0597
#> Class  2        0.4917    0.2343       0.5985       0.4366
#> Class  3        0.8998    0.6970       0.3552       0.1549
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.2784   0.0580      0.0785      0.0605       0.1218
#> Class  2       0.3864   0.1694      0.3152      0.1565       0.4230
#> Class  3       0.6757   0.4776      0.1578      0.0626       0.6315
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.0107      0.0994      0.0377
#> Class  2   0.1543      0.2091      0.0619
#> Class  3   0.3814      0.1550      0.0851
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.3717    0.2438       0.0475       0.0001
#> Class  2        0.4551    0.2885       0.9284       0.7962
#> Class  3        0.7989    0.6336       0.3401       0.1839
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.4425    0.1577       0.2093       0.1432
#> Class  2        0.4917    0.2372       0.5597       0.4604
#> Class  3        0.8995    0.6945       0.4113       0.2548
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.2787   0.0595      0.1785      0.1445       0.1221
#> Class  2       0.3866   0.1722      0.3853      0.2563       0.4232
#> Class  3       0.6754   0.4779      0.2650      0.1477       0.6313
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.0110      0.2047      0.1060
#> Class  2   0.1570      0.3087      0.1467
#> Class  3   0.3831      0.2624      0.1787
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.0808  0.248    2.0313    1.9101

symptoms.lca4random <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,nclass=4,blocksize=4,constload=FALSE)
summary(symptoms.lca4random)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik  Link
#>         4 9417.299 9726.582 9488.299 -4637.649 -4637.852 Logit
#> Class probabilities 
#> Class  1 Class  2 Class  3 Class  4 
#>   0.4017   0.1672   0.1929   0.2382 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.7002    0.4968       0.0726       0.0272
#> Class  2        0.5221    0.3135       0.9905       0.8801
#> Class  3        0.4729    0.2867       0.8562       0.7516
#> Class  4        0.4147    0.3047       0.0735       0.0000
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.8252    0.5147       0.1510       0.0424
#> Class  2        0.5200    0.2514       0.2391       0.0524
#> Class  3        0.5615    0.3238       0.9998       1.0000
#> Class  4        0.4520    0.1872       0.1036       0.0688
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.5392   0.2891      0.0078      0.0025       0.4335
#> Class  2       0.3264   0.1112      0.3102      0.1510       0.4491
#> Class  3       0.4838   0.2327      0.6224      0.3631       0.4265
#> Class  4       0.3574   0.1139      0.2979      0.1946       0.2113
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1740      0.0260      0.0128
#> Class  2   0.1895      0.1979      0.1004
#> Class  3   0.1119      0.3756      0.1656
#> Class  4   0.0993      0.4168      0.1858
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.6791    0.4975       0.1405       0.0598
#> Class  2        0.5194    0.3470       0.9718       0.8126
#> Class  3        0.4762    0.3238       0.7763       0.6902
#> Class  4        0.4247    0.3394       0.1418       0.0000
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.7994    0.5118       0.2310       0.0856
#> Class  2        0.5176    0.2925       0.3104       0.1013
#> Class  3        0.5542    0.3558       0.9991       1.0000
#> Class  4        0.4577    0.2322       0.1798       0.1247
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.5345   0.3259      0.0238      0.0068       0.4414
#> Class  2       0.3454   0.1526      0.3662      0.2204       0.4552
#> Class  3       0.4858   0.2754      0.5848      0.4003       0.4352
#> Class  4       0.3735   0.1557      0.3569      0.2625       0.2369
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.2191      0.0652      0.0313
#> Class  2   0.2344      0.2751      0.1650
#> Class  3   0.1533      0.4138      0.2349
#> Class  4   0.1389      0.4428      0.2543
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.7844 1.0926    1.6055    1.4514

symptoms.lca5random <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,nclass=5,blocksize=4,constload=FALSE)
summary(symptoms.lca5random)
#>   Classes      AIC     BIC     AIC3    logLik penlogLik  Link
#>         5 9347.853 9731.19 9435.853 -4585.926 -4586.138 Logit
#> Class probabilities 
#> Class  1 Class  2 Class  3 Class  4 Class  5 
#>   0.1524   0.1163   0.1337   0.4878   0.1099 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5744    0.3744       1.0000       0.7889
#> Class  2        0.5520    0.2952       0.5568       0.3817
#> Class  3        0.5840    0.3923       0.3122       0.0850
#> Class  4        0.5895    0.3889       0.0597       0.0167
#> Class  5        0.4422    0.3235       1.0000       0.9999
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.5686    0.3409       0.2182       0.0777
#> Class  2        0.6887    0.5063       0.9999       1.0000
#> Class  3        0.5496    0.2286       0.4126       0.3216
#> Class  4        0.7125    0.3538       0.1115       0.0106
#> Class  5        0.5593    0.2493       1.0000       0.8699
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.3848   0.1811      0.2501      0.1144       0.4727
#> Class  2       0.4946   0.1364      0.2682      0.0965       0.3977
#> Class  3       0.4640   0.1563      0.6840      0.4644       0.3097
#> Class  4       0.4356   0.1810      0.0253      0.0226       0.3274
#> Class  5       0.5156   0.2992      0.8452      0.6611       0.5141
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.2144      0.1598      0.0535
#> Class  2   0.0816      0.0288      0.0000
#> Class  3   0.1765      0.8251      0.5918
#> Class  4   0.1006      0.0603      0.0214
#> Class  5   0.1498      0.7536      0.4011
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5615    0.4061       0.9999       0.7763
#> Class  2        0.5429    0.3439       0.5511       0.3887
#> Class  3        0.5695    0.4197       0.3292       0.0941
#> Class  4        0.5740    0.4172       0.0725       0.0191
#> Class  5        0.4523    0.3665       0.9999       0.9999
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.5566    0.3802       0.2395       0.0862
#> Class  2        0.6586    0.5047       0.9999       1.0000
#> Class  3        0.5409    0.2875       0.4213       0.3315
#> Class  4        0.6795    0.3903       0.1303       0.0121
#> Class  5        0.5490    0.3056       0.9999       0.8585
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.4044   0.2438      0.2705      0.1251       0.4775
#> Class  2       0.4955   0.1988      0.2879      0.1063       0.4152
#> Class  3       0.4703   0.2195      0.6672      0.4666       0.3401
#> Class  4       0.4468   0.2438      0.0317      0.0257       0.3555
#> Class  5       0.5129   0.3471      0.8241      0.6520       0.5116
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.2749      0.1811      0.0600
#> Class  2   0.1353      0.0359      0.0000
#> Class  3   0.2394      0.8036      0.5863
#> Class  4   0.1587      0.0732      0.0243
#> Class  5   0.2129      0.7331      0.4070
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      1.0108 1.3599      0.71    0.5295

symptoms.lca1random2 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,level2=TRUE,nclass=1, level2size=4,constload=FALSE,
  penalty=0.1)
summary(symptoms.lca1random2)
#>   Classes      AIC    BIC     AIC3    logLik penlogLik  Link
#>         1 9504.222 9595.7 9525.222 -4731.111 -4733.267 Logit
#> Class probabilities 
#> Class  1 
#>        1 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5543    0.3988       0.1054        0.083
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.6265     0.376        0.067       0.0643
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.4553   0.2416      0.0037      0.0155       0.3943
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1851      0.0026      0.0092
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.5538     0.401       0.4078       0.2959
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.6254    0.3787       0.3867        0.275
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.4557   0.2459      0.2698      0.1772       0.3952
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1895      0.2567      0.1484
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.1433  0.224      6.52    2.9941
#> Tau 
#> 0.9417

symptoms.lca2random2 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,level2=TRUE,nclass=2, level2size=4,constload=FALSE, penalty=0.1)
summary(symptoms.lca2random2)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik  Link
#>         2 9198.261 9363.793 9236.261 -4561.131 -4563.424 Logit
#> Class probabilities 
#> Class  1 Class  2 
#>   0.4296   0.5704 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.8031    0.6242       0.1255       0.0611
#> Class  2        0.3674    0.2319       0.0834       0.0931
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.8731    0.6722       0.1586       0.0542
#> Class  2        0.4407    0.1580       0.0297       0.0631
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.6525   0.4576      0.0046      0.0100       0.6216
#> Class  2       0.3089   0.0853      0.0023      0.0169       0.2276
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.3771      0.0030      0.0157
#> Class  2   0.0476      0.0017      0.0045
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.8021    0.6220       0.4185       0.2766
#> Class  2        0.3680    0.2356       0.3997       0.3107
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.8723    0.6693       0.4297       0.2675
#> Class  2        0.4410    0.1616       0.3547       0.2792
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.6518   0.4584      0.2839      0.1593       0.6210
#> Class  2       0.3097   0.0879      0.2590      0.1890       0.2285
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.3793      0.2679      0.1846
#> Class  2   0.0493      0.2474      0.1210
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.1006 0.2039    6.6664    3.0751
#> Tau 
#> 0.9544

# these can take considerable execution time

symptoms.lca3random2 <- randomLCA(symptoms[,1:16],freq=symptoms$Freq,
  random=TRUE,level2=TRUE,nclass=3, level2size=4,constload=FALSE,
  penalty=0.1)
summary(symptoms.lca3random2)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik  Link
#>         3 9159.706 9399.292 9214.706 -4524.853 -4527.507 Logit
#> Class probabilities 
#> Class  1 Class  2 Class  3 
#>   0.3128   0.2706   0.4166 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.3241    0.1426       0.9717       0.6494
#> Class  2        0.4414    0.3383       0.0003       0.0029
#> Class  3        0.8009    0.6285       0.0566       0.0434
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.4312    0.1144       0.7348       0.4045
#> Class  2        0.4596    0.2047       0.0008       0.0062
#> Class  3        0.8826    0.6864       0.0728       0.0408
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.3531   0.1199      0.0406      0.0518       0.4369
#> Class  2       0.2516   0.0396      0.0006      0.0121       0.0025
#> Class  3       0.6659   0.4620      0.0027      0.0081       0.6220
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1174      0.0044      0.0075
#> Class  2   0.0003      0.0062      0.0091
#> Class  3   0.3508      0.0015      0.0132
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.3274    0.1554       0.6566       0.5564
#> Class  2        0.4426    0.3480       0.1729       0.0900
#> Class  3        0.7968    0.6204       0.3744       0.2383
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.4326    0.1260       0.5464       0.4645
#> Class  2        0.4605    0.2183       0.2084       0.1212
#> Class  3        0.8792    0.6755       0.3862       0.2337
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.3559   0.1317      0.3592      0.2517       0.4382
#> Class  2       0.2556   0.0451      0.1992      0.1558       0.0026
#> Class  3       0.6628   0.4645      0.2507      0.1342       0.6196
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.1292      0.2689      0.1306
#> Class  2   0.0003      0.2822      0.1401
#> Class  3   0.3599      0.2297      0.1604
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.2114 0.3951    6.1728    2.8442
#> Tau 
#> 0.9707

Repeating for up to five classes or when the BIC increases gives the following results. It should be noted that the two level models can take considerable time to fit. This is due to the relatively large number of quadrature points required for this data, as a consequence of a large number of patterns consisting entirely of zeroes. As for the previous examples we can form the BIC values into a table, where bic, bic.random and bic.random2 are the BIC values from the standard, random effects and 2 level random effects latent class models.

symptoms.bic <- data.frame(class=1:5,
  bic=c(BIC(symptoms.lca1),BIC(symptoms.lca2),BIC(symptoms.lca3), BIC(symptoms.lca4),
    BIC(symptoms.lca5)),
  bic.random=c(BIC(symptoms.lca1random),BIC(symptoms.lca2random),
    BIC(symptoms.lca3random),BIC(symptoms.lca4random),BIC(symptoms.lca5random)),
  bic.random2=c(BIC(symptoms.lca1random2),BIC(symptoms.lca2random2),
    BIC(symptoms.lca3random2),NA,NA))
kable(symptoms.bic, format = "html")
Table 3: BIC for Symptoms models
class bic bic.random bic.random2
1 10844.187 9987.510 9595.700
2 10123.811 9769.743 9363.793
3 9960.839 9700.262 9399.292
4 9910.742 9726.582 NA
5 9918.924 9731.190 NA

This shows the optimal model is the 2 Class model with random effects for both subject and period.

summary(symptoms.lca2random2)
#>   Classes      AIC      BIC     AIC3    logLik penlogLik  Link
#>         2 9198.261 9363.793 9236.261 -4561.131 -4563.424 Logit
#> Class probabilities 
#> Class  1 Class  2 
#>   0.4296   0.5704 
#> Conditional outcome probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.8031    0.6242       0.1255       0.0611
#> Class  2        0.3674    0.2319       0.0834       0.0931
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.8731    0.6722       0.1586       0.0542
#> Class  2        0.4407    0.1580       0.0297       0.0631
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.6525   0.4576      0.0046      0.0100       0.6216
#> Class  2       0.3089   0.0853      0.0023      0.0169       0.2276
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.3771      0.0030      0.0157
#> Class  2   0.0476      0.0017      0.0045
#> Marginal Outcome Probabilities 
#>          Nightcough.13 Wheeze.13 Itchyrash.13 FlexDerma.13
#> Class  1        0.8021    0.6220       0.4185       0.2766
#> Class  2        0.3680    0.2356       0.3997       0.3107
#>          Nightcough.45 Wheeze.45 Itchyrash.45 FlexDerma.45
#> Class  1        0.8723    0.6693       0.4297       0.2675
#> Class  2        0.4410    0.1616       0.3547       0.2792
#>          Nightcough.6 Wheeze.6 Itchyrash.6 FlexDerma.6 Nightcough.7
#> Class  1       0.6518   0.4584      0.2839      0.1593       0.6210
#> Class  2       0.3097   0.0879      0.2590      0.1890       0.2285
#>          Wheeze.7 Itchyrash.7 FlexDerma.7
#> Class  1   0.3793      0.2679      0.1846
#> Class  2   0.0493      0.2474      0.1210
#> Loadings 
#>  Nightcough Wheeze Itchyrash FlexDerma
#>      0.1006 0.2039    6.6664    3.0751
#> Tau 
#> 0.9544

The marginal outcome probabilities are plotted in Figure 5 as follows.

plot(symptoms.lca2random2, type="b",
  scales=list(x=list(at=1:4, labels=c(6, 12, 18, 24))), pch = 1:4,
  xlab="Period",
  key = list(corner = c(0.75, .85),
    text = list(c("Night Cough", "Wheeze", "Itch Rash", "Flex. Derma.")),
    points = list(pch= 1:4), 
    col = trellis.par.get()$superpose.symbol$col[1:4], border = TRUE))
Marginal outcome probabilities for 2 class latent class for symptoms data.

Figure 5: Marginal outcome probabilities for 2 class latent class for symptoms data.

It can be seen that in Class 1 the outcome probabilities are greater than for Class 2, with the difference greatest for the two respiratory symptoms night cough and wheeze. Also over time the outcome probabilities decrease. Similarly to the dentistry example the conditional probabilities can be obtained using calcCond2Prob and plotted.

4 Summary

It has been shown how the randomLCA package may be used to determine classes of subjects based on observed binary data, and how this may be used to determine sensitivity and specificity for diagnostic tests. In the dentistry and symptoms examples an assumption of heterogeneous classes, where the outcome probabilities are allowed to vary within a class was shown to be an improvement over traditional latent class analysis. The randomLCA package also produces a range of plots for describing the classes, allows for bootstrapped standard errors, calculation of marginal outcome probabilities when using random effects models and the use of penalised likelihood. A further extension to randomLCA would be extension to latent class regression models (Dayton and Macready 1988), where the class probabilities are determined by the covariates. Another extension is to allow for ordinal data using a Graded Response model (Samejima 1970). A possible extension to the package is to allow for other data types. However for the random effects models it is difficult to extend the models except for ordinal data.

References

Bartholomew, D J, F Steele, I Moustaki, and J I Galbraith. 2002. The Analysis and Interpretation of Multivariate Data for Social Scientists. Chapman & Hall/CRC.
Beath, K. J, and G. Z Heller. 2009. “Latent Trajectory Modelling of Multivariate Binary Data.” Statistical Modelling 9 (3): 199–213. https://doi.org/10.1177/1471082X0800900302.
Davison, A. C., and D. V. Hinkley. 1997. Bootstrap Methods and Their Application. Cambridge University Press.
Dayton, C M, and G B Macready. 1988. “Concomitant-Variable Latent-Class Models.” Journal of the American Statistical Association 83 (401): 173–78. https://doi.org/10.2307/2288938.
Dias, J G. 2006. “Latent Class Analysis and Model Selection.” In From Data and Information Analysis to Knowledge, edited by M Spiliopoulou, R Kruse, C Borgelt, A Nürnberger, and W Gaul. Springer-Verlag.
Firth, D. 1993. “Bias Reduction of Maximum Likelihood Estimates.” Biometrika 80 (1): 27–38. https://doi.org/10.1093/biomet/80.1.27.
Galindo Garre, Francisca, and Jeroen K. Vermunt. 2006. “Avoiding Boundary Estimates in Latent Class Analysis by Bayesian Posterior Mode Estimation.” Behaviormetrika 33 (1): 43–59. https://doi.org/10.2333/bhmk.33.43.
Golden, R. R. 1982. “A Taxometric Model for the Detection of a Conjectured Latent Taxon.” Multivariate Behavioral Research 17: 389–416. https://doi.org/10.1207/s15327906mbr1703_6.
Kraemer, Helena Chmura. 1992. Evaluating Medical Tests: Objective and Quantitative Guidelines. SAGE Publications.
Langeheine, R., and F. van de Pol. 1990. “A Unifying Framework for Markov Modelling in Discrete Space and Discrete Time.” Sociological Methods and Research 18 (4): 416–41. https://doi.org/10.1177/0049124190018004002.
Lazarsfeld, P. F., and N. W. Henry. 1968. Latent Structure Analysis. Houghton Mifflin.
Lin, Ting Hsiang, and C Mitchell Dayton. 1997. “Model Selection Information Criteria for Non-Nested Latent Class Models.” Journal of Educational and Behavioral Statistics 22 (3): 249–64. https://doi.org/10.3102/10769986022003249.
Linzer, Drew A, and Jeffrey B Lewis. 2011. “poLCA: An r Package for Polytomous Variable Latent Class Analysis.” Journal of Statistical Software 42 (10): 1–29. https://doi.org/10.18637/jss.v042.i10.
Little, R. J. A., and D. B. Rubin. 2002. Statistical Analysis with Missing Data. 2nd ed. John Wiley & Sons.
Liu, Qing, and Donald A Pierce. 1994. “A Note on Gauss-Hermite Quadrature.” Biometrika 81 (3): 624–29. https://doi.org/10.2307/2337136.
McHugh, R. B. 1956. “Efficient Estimation and Local Identification in Latent Class Analysis.” Psychometrika 21 (4): 331–47. https://doi.org/10.1007/BF02296300.
McLachlan, G. J. 1987. “On Bootstrapping the Likelihood Ratio Test Statistic for the Number of Components in a Normal Mixture.” Applied Statistics 36 (3): 318–24. https://doi.org/10.2307/2347790.
McLachlan, G. J., and D. Peel. 2000. Finite Mixture Models. John Wiley & Sons.
Mihrshahi, Seema, Jennifer K Peat, Karen Webb, et al. 2001. “The Childhood Asthma Prevention Study (CAPS): Design and Research Protocol of a Randomized Trial for the Primary Prevention of Asthma.” Controlled Clinical Trials 22 (3): 333–54. https://doi.org/10.1016/S0197-2456(01)00112-X.
Muthén, Bengt. 2006. “Should Substance Use Disorders Be Considered as Categorical or Dimensional?” Addiction 101 (Supplement 1): 6–16. https://doi.org/10.1111/j.1360-0443.2006.01583.x.
Muthén, B, and K Shedden. 1999. “Finite Mixture Modeling with Mixture Outcomes Using the EM Algorithm.” Biometrics 55 (2): 463–69. https://doi.org/10.1111/j.0006-341X.1999.00463.x.
Muthén, L. K, and B. O. Muthén. 2015. Mplus User’s Guide. Seventh. Muthén & Muthén.
Nyholt, D. R., N. G. Gillespie, A. C. Heath, K. R. Merikangas, D. L. Duffy, and N. G. Martin. 2004. “Latent Class and Genetic Analysis Does Not Support Migraine with Aura and Migraine Without Aura as Separate Entities.” Genetic Epidemiology 26: 231–44. https://doi.org/10.1002/gepi.10311.
Nylund, Karen L, and Bengt O Muthén. 2007. “Deciding on the Number of Classes in Latent Class Analysis and Growth Mixture Modeling: A Monte Carlo Simulation Study.” Structural Equation Modeling 14 (4): 535–69. https://doi.org/10.1080/10705510701575396.
Pepe, Margaret Sullivan. 2003. The Statistical Evaluation of Medical Tests for Classification and Prediction. Oxford University Press.
Pickles, A., and A. Angold. 2003. “Natural Categories or Fundamental Dimensions: On Carving Nature at the Joints and the Rearticulation of Psychopathology.” Development and Psychopathology 15: 529–51. https://doi.org/10.1017/S0954579403000282.
Qu, Y, M Tan, and M H Kutner. 1996. “Random Effects Models in Latent Class Analysis for Evaluating Accuracy of Diagnostic Tests.” Biometrics 52 (3): 797–810. https://doi.org/10.2307/2533043.
R Core Team. 2016. R: A Language and Environment for Statistical Computing. R Founda- tion for Statistical Computing.
Rabe-Hesketh, Sophia, Anders Skrondal, and Andrew Pickles. 2005. “Maximum Likelihood Estimation of Limited and Discrete Dependent Variable Models with Nested Random Effects.” Journal of Econometrics 128 (2): 301–23. https://doi.org/10.1016/j.jeconom.2004.08.017.
Rindskopf, D., and W. Rindskopf. 1986. “The Value of Latent Class Analysis in Medical Diagnosis.” Statistics in Medicine 5: 21–27. https://doi.org/10.1002/sim.4780050105.
Rubin, Donald B., and N. Schenker. 1987. “Logit-Based Interval Estimation for Binomial Data Using the Jeffreys Prior.” Sociological Methodology 17: 131–44. https://doi.org/10.2307/271031.
Samejima, Fumiko. 1970. “Estimation of Latent Ability Using a Response Pattern of Graded Scores.” Psychometrika 35 (1): 139–39. https://doi.org/10.1007/BF02290599.
Skrondal, Anders, and Sophia Rabe-Hesketh. 2004. Generalized Latent Variable Modeling: Multilevel, Longitudinal and Structural Equation Models. Chapman & Hall/CRC.
Uebersax, J. S. 1999. “Probit Latent Class Analysis with Dichotomous or Ordered Category Measures: Conditional Independence/Dependence Models.” Applied Psychological Measurement 23 (4): 283–97. https://doi.org/10.1177/01466219922031400.
Vacek, P. M. 1985. “The Effect of Conditional Dependence on the Evaluation of Diagnostic Tests.” Biometrics 41: 959–68. https://doi.org/10.2307/2530967.
Vermunt, Jeroen K., and Jay Magidson. 2013. LG-Syntax User’s Guide: Manual for Latent GOLD 5.0 Syntax Module. Statistical Innovations Inc.
White, Arthur, and Thomas Brendan Murphy. 2014. “BayesLCA: An r Package for Bayesian Latent Class.” Journal of Statistical Software 61 (13): 1–28. https://doi.org/10.18637/jss.v061.i13.
Young, M. A. 1982. “Evaluating Diagnostic Criteria: A Latent Class Paradigm.” Journal of Psychiatric Research 17 (3): 285–96. https://doi.org/10.1016/0022-3956(82)90007-3.