Two-level latent class analysis
In latents: Multilevel Latent Profile, Class and Transition Models

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.5, fig.align = "center", dpi = 96)
old_options <- options(digits = 3, width = 100)

Two-level latent class analysis (LCA) is used when the indicators are categorical: yes/no answers, ratings on a few ordered options, or unordered choices. Like multilevel latent profile analysis, it estimates latent classes of observations, called profiles in the output, and latent group classes that differ in how probable each profile is for their observations. The difference is in how a profile is described: by the probability of each answer to each item, instead of a mean and a variance. multilca() fits LCA, treating every indicator as categorical; multilpa() with some indicators named in categorical fits a mixed model with both kinds of indicator.

What the package supports

Item types. Binary, ordinal and nominal items share one parameterization: a free probability for every category in every profile. Numeric, integer, logical, character and factor columns are accepted; factor levels keep their order, other types are sorted. Ordered items can instead be modelled as ordinal (an adjacent-category logit with far fewer parameters), and counts as count (Poisson); see below.
Mixed measurement. Continuous and categorical indicators in one model.
Missing answers. missing = "fiml" uses the observed answers of each row.
Standard errors. Observed-information and robust (vcov_type = "robust") standard errors through parameter_inference().
Model comparison. enumerate_classes() and, for complete data, bootstrap_lrt().
External variables. three_step() and r3step() work on a categorical fit.
Membership covariates. profile_covariates and group_covariates work with categorical indicators, and parameter_inference() reports their standard errors.
Transitions. lta() accepts categorical indicators.
Bounds. min_probability keeps every response probability above a lower bound.

Data

student_esm holds experience-sampling data: 2,582 prompts from 100 university students, answered at times when the student had not been studying. student identifies the group. Eight yes/no items record the leisure activities done since the previous prompt, and four items rate affect from 1 to 7. day is the study day, from 0 to 13. ?student_esm gives the source.

library(latents)

To keep the example short, we call subset() with day <= 6, which keeps the first week and all 100 students.

first_week <- subset(student_esm, day <= 6)

The eight activity items are used in several fits, so we store their names once.

activities <- c("time_with_friends", "on_social_media", "tv_video_games",
                "listened_music", "sports", "walking", "reading",
                "part_time_job")

Fitting the model

To fit the model, we call multilca() with vars for the categorical indicators, id for the group, n_profiles, n_group_classes, n_starts, seed, and tol for a tighter convergence tolerance, which standard errors need.

lca <- multilca(first_week, vars = activities, id = "student",
                n_profiles = 2, n_group_classes = 2, n_starts = 3,
                tol = 1e-10, seed = 1)
lca

To obtain the response probabilities, we call get_results() with "responses". There is one row per profile, item and category; probability is the chance of that answer in that profile.

get_results(lca, "responses")

To see them, we call plot() with what = "responses": one line per profile across the items.

plot(lca, what = "responses")

To obtain the profile probabilities of each group class, we call get_results() with "profile_probabilities", and to compare the classes, plot() with what = "probabilities".

get_results(lca, "profile_probabilities")

plot(lca, what = "probabilities")

To obtain classification certainty, we call get_results() with "entropy".

get_results(lca, "entropy")

Mixed indicators

To combine both kinds, we call multilpa() with all indicators in vars and only the categorical ones in categorical. Here three affect ratings are continuous and the activities categorical. worried is left out: about half its ratings are 1, which can push a profile's variance to its lower bound.

mixed <- multilpa(first_week,
                  vars = c("happy", "relaxed", "exhausted", activities),
                  categorical = activities, id = "student", n_profiles = 2,
                  n_group_classes = 2, n_starts = 3, seed = 1)

To obtain the means and variances of the continuous indicators, we call get_results() on the fit; for the answer probabilities of the categorical ones, get_results() with "responses".

get_results(mixed)
get_results(mixed, "responses")

Ordinal and count indicators

The affect ratings run from 1 to 7. As categorical items each would have six free probabilities in every profile; as ordinal items each has six category intercepts shared by every profile and one location per profile, which says how far up the scale that profile sits. We pass their names as ordinal.

affect <- c("happy", "relaxed", "worried", "exhausted")
ordinal <- multilpa(first_week, vars = c(affect, activities), id = "student",
                    ordinal = affect, categorical = activities, n_profiles = 2,
                    n_group_classes = 2, n_starts = 3, seed = 1)
ordinal

get_results() with "ordinal" gives each profile's probability of every rating and its location.

head(get_results(ordinal, "ordinal"), 7)

Indicators holding counts (number of messages sent, of absences) go in count: each is Poisson with one mean per profile, read with get_results(fit, "count_means"). Counts more variable than a Poisson within a profile take count_model = "negative_binomial", which adds a dispersion per profile (or one shared, count_dispersion = "equal"). Both types work with membership covariates and in lta() as well.

Reference

The calls below list other functions that work with categorical and mixed models. Each comment states what the call returns.

get_results(lca, "classification")               # class sizes and average posteriors
get_results(lca, "residuals")                    # residual association within profiles
parameter_inference(lca)                         # standard errors and intervals
enumerate_classes(first_week, vars = activities, categorical = activities,
                  id = "student", n_profiles = 2:3,
                  n_group_classes = 1:2, seed = 1)  # compare numbers of classes (slow)
plot(mixed)                                      # means of the continuous indicators
plot(mixed, what = "responses")                  # answers to the categorical ones
multilca(first_week, vars = activities, id = "student", n_profiles = 2,
         n_group_classes = 2, missing = "fiml")  # incomplete answers, if any

Assumptions

Items are independent within a profile.
Every category of every item has its own probability in each profile, so rare answers can push probabilities to their lower bound (min_probability).
Categories are taken from the observed values; factor levels keep their order.