How can you implement latent class analysis with distal outcomes in R? I am looking to fit a fairly straightforward latent class analysis (LCA) model to derive phenotypes / clusters of a disease (in R). My dataset contains the manifest variables used to derive the clusters (as in any other LCA model), which are categorical.
Found packages that do the trick in deriving the classes and that can also incorporate covariates in a regression to see whether they are related to the classes (i.e., poLCA).
My problem is that I haven't found a way to incorporate distal outcomes in any R packages that deal with LCA. Does anyone have any experience in implementing distal outcomes in any existing R packages or any way to get around to this?
 A: I don't know of many statistical software packages in general that implement latent class analysis with distal outcomes.
Background
Just to clarify for readers: in latent class analysis, we're saying that we have a number of indicators $X_j$. We assume there's a latent, unordered categorical variable that causes responses to those indicators, and that it has $k$ classes. We iteratively increase $k$, and for each value of $k$, we estimate the means of all the $X_j$s. Formally, for the usual application of latent class analysis with binary indicators, we'd estimate:
\begin{align}
P(X_j = 1 | C = k) &= {\rm logit}^{-1}(\alpha_{jk})  \\[5pt]
P(C = k) &= \frac{\exp(\gamma_k)}{\sum^C_{j=1}\gamma_j}
\end{align}
$\gamma_1 = 0$ because you need a base class. The above is a multinomial regression model with intercepts only. That is, the likelihood function for LCA models contains a multinomial model. For a model with binary indicators
$\sum^K_{k=1}P(Class = k)\prod^M_{m=1}P(Indicator_m=1|Class=k)[1-P(Indicator_m=1|Class=k)^{1-I(Indicator_m=1|Class=k)]}$
Variables related to the latent class
An obvious step after you fit an LCA model is asking if some variable(s) is/are related to latent class membership. For example, how many percent of each class are male. Formally, say we're interested in a vector of some variables $Y$, and we want to see something like $E(Y | Class = k)$.
If we had some observed categorical variable, we would just tabulate $Y$ by that variable. One of the harder things to grasp about LCA is that you can't do this. You do not know which class each observation belongs to. You do know the vector of latent classes to which each observation belongs to, e.g.




Obs.
P(Class = 1)
P(Class = 2)
P(Class = 3)




Mrs. Smith
0.1
0.2
0.7


Mrs. Wang
0.3
0.3
0.4




You can assume modal class assignment, i.e. you assume that each observation is in the class with the highest posterior probability of membership. However, if a lot of people look like Mrs. Wang, we can see the problem - you aren't actually very sure which class Mrs. Wang belongs to. Whatever relationships exist between class membership and $Y$ may be distorted. If the model's entropy is high, i.e. you're relatively sure which classes people belong to, it may be justifiable to just do modal class assignment. If it were me, I would note this in limitations. If your entropy is relatively low, I'd prefer some sort of alternate method.
Jeroen Vermunt and colleagues have written extensively about this problem. One example paper is here (ungated). He references Bolck, Croon, and Hagenaars (2004, gated), who showed that modal class assignment should bias the relationships between latent class and $Y$ downwards. Two other simple solutions would seem to be random assignment (we could multiply impute the class memberships and use Rubin's rules) or proportional assignment (use the probability of membership as a weight in regular tabulation), but they showed that these also are problematic.
BCH and Vermunt et al have both proposed ways to correct for the bias in the results. One of these is the three-step procedure referenced in the question. Because Vermunt's math is too complex for me to follow, I won't discuss it and I commend it to your reading.
Latent class regression
In latent class regression, we're asking how does the level of some covariate $Z$ change the probability of being in each latent class. That is, we can add $X$s and $\beta$s to the multinomial model.
$$P(C = k | Z) = \frac{\exp(\gamma_k + \gamma Z)}{\sum^C_{j=1}\exp(\gamma_j + \gamma Z)}$$
This is one way to get at the association between latent class and some vector of variables. Covariates that predict latent class membership are usually denoted as $Z$ in the literature. Distal covariates are not used at all in estimating the latent class model, and I'm denoting them $Y$ above. Substantively, though, there isn't a difference in that you just want the relationship between the latent class and some variables.
We have $P(Class|Z)$. We can present that on its own terms, although sometimes multinomial logit models can be hard to interpret for those not used to them. If we would really rather present $E(Z|Class)$, though, then for binary $Z$s, we do already have the marginal probabilities of $Z$ and the marginal probability of being in each latent class. Stephanie Lanza and colleagues (2013, ungated) showed that we can use Bayes theorem on this problem. This doesn't require additional software and can be done in Excel or by hand.
They also showed an extension for continuous $Z$s. The thing is that if $Z$ is continuous, you need some assumption about the functional form of $f(Z)$. The authors wrote add-on packages for SAS (PROC LCA) and Stata that use kernel density estimation for $f(Z)$. I'm not sure how to implement that in R, and I'm not sure which (if any) latent class packages implement this.
poLCA and Stata can do latent class regression. There are probably more R packages that can do it.
Why would you not want to use latent class regression?
My understanding from people who regularly use multinomial logit is that it can be a fussy model in terms of convergence. For example, if some covariates or covariate combinations are not present in one latent class, the model can refuse to converge. In Stata at least, there's no clear error message, and you are thus left to infer that this may have happened.
Also, you probably started with a plain LCA model, then added covariates. (You can save the parameters from your chosen solution and supply them to the latent class regression model as start values, which can reduce estimation time.) What do you do if the regression model changes the substantive nature of the classes identified? I'm not sure that there's a universally accepted solution.
Software discussion
I don't have direct experience with MPlus or Latent Gold, as they are both specialist packages. I believe that both may implement the three-step procedure.
The R packages poLCA and flexmix don't implement the three-step procedure. Neither does Stata. I don't know about SAS. The relative rarity of software that does this procedure could be given as a reason that you didn't implement, if a reviewer really presses you.
A: Have you tried Structural Equation Modeling (SEM) packages? lavaan, sem, and OpenMX come to mind. I'm not sure about LCA in particular, but as I understand it, LCA is a subset of SEM, and if memory serves me the SEM packages that I've looked at in the past support both directions of outcomes and latent variables.
