Multilevel latent profile analysis

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)

Multilevel latent profile analysis (MLPA) estimates latent profiles at the observation level and, at a higher level, identifies latent group classes that differ in how probable each profile is for their observations. Each profile is characterized by a mean and variance for each indicator and is defined identically across all group classes. The group classes therefore differ only in the probabilities associated with the profiles.

Data

course_engagement has one row per course enrolment: 1,422 enrolments from 106 students. multilpa() needs two kinds of variables.

The data are simulated, and two further variables record the structure that generated them: engagement for each enrolment and student_type for each student.

library(latents)

Clustering within groups

To check how much each indicator varies between groups, we call descriptives() with vars for the indicators and id for the group. icc is the share of an indicator's variance that lies between groups.

descriptives(course_engagement,
             vars = c("browse", "lectures", "forum_read", "forum_post",
                      "attendance"),
             id = "student")

Fit

To fit the model, we call multilpa() with vars for the indicators, id for the group, n_profiles for the number of profiles, n_group_classes for the number of group classes, n_starts for the number of random EM starts, and seed for reproducible starts. Printing the fit gives the profile means.

fit <- multilpa(course_engagement,
                vars = c("browse", "lectures", "forum_read", "forum_post",
                         "attendance"),
                id = "student", n_profiles = 2, n_group_classes = 2,
                n_starts = 10, seed = 1)
fit

To check whether the starts agree, we call get_results() with "starts". Starts that reach the same log likelihood support the solution.

get_results(fit, "starts")

Profiles

To obtain the profile means and variances, we call get_results() on the fit. There is one row per profile and indicator.

get_results(fit)

To see the profiles, we call plot() on the fit. Each line is a profile's mean on each indicator, and its label gives the profile's share of the cases.

plot(fit)

Group classes

To obtain the profile probabilities of each group class, we call get_results() with "profile_probabilities". probability is the chance that an observation of a group in that class falls in each profile, and group_class_probability is the share of groups in the class.

get_results(fit, "profile_probabilities")

To compare the group classes, we call plot() with what = "probabilities": one line per group class across the profiles.

plot(fit, what = "probabilities")

Class sizes

To obtain class sizes, we call get_results() with "counts". effective_count sums the posterior probabilities, so it need not be a whole number.

get_results(fit, "counts")

Classification

To obtain classification certainty, we call get_results() with "entropy". Relative entropy is reported for observations and for groups: 1 means every case is assigned with certainty, 0 means no separation.

get_results(fit, "entropy")

To add each observation's assigned profile and group class to the data, we call get_results() with "assignments" and data for the data frame. It returns the data in its original row order with five added columns: profile and group_class hold the most probable classes, uncertainty is one minus the posterior probability of the assigned profile, and the posterior_profile_ columns hold the posterior probability of each profile. Assignment takes the most probable class, so profile and group_class alone ignore classification uncertainty.

get_results(fit, "assignments", data = course_engagement) |> head()

Recovery of the known structure

To compare the fitted classes with known labels, we call get_results() with "assignments", data for the data frame, and truth for the columns that hold the labels. engagement is compared with the profiles and student_type with the group classes. Class numbers are arbitrary, so the table is read by its largest cells.

get_results(fit, "assignments", data = course_engagement,
            truth = c("engagement", "student_type"))

Reference

The calls below list the other tables, functions and views available for a fitted model. Each comment states what the call returns.

# Tables: get_results(fit, what)
get_results(fit, "all")                          # every table, as a named list
get_results(fit, "covariances")                  # within-profile residual covariances
get_results(fit, "posteriors")                   # one row per observation and profile
get_results(fit, "posteriors", format = "wide")  # one row per observation
get_results(fit, "group_posteriors")             # one row per group and group class
get_results(fit, "classification")               # class sizes and average posteriors
get_results(fit, "average_posteriors")           # mean posterior by assigned class
get_results(fit, "classification_errors")        # P(assigned class | true class)
get_results(fit, "bch_weights")                  # weights for three-step analysis
get_results(fit, "residuals")                    # residual correlations within profiles
get_results(fit, "information_criteria")         # AIC, BIC and the other criteria
get_results(fit, "model")                        # one row describing the fit
get_results(fit, "data")                         # the variables used in the fit
as.data.frame(fit)                               # the profile table

# Summaries and diagnostics
descriptives(fit)                    # indicator summaries with the ICC
descriptives(fit, by = "profile")    # the same, by assigned profile
diagnostics(fit)                     # every classification diagnostic
report(fit)                          # summary, diagnostics and plots together

# Inference and likelihood
parameter_inference(fit)             # standard errors and Wald intervals
confint(fit)                         # confidence intervals
vcov(fit)                            # covariance matrix of the estimates
logLik(fit)                          # log likelihood
AIC(fit)                             # Akaike information criterion
BIC(fit)                             # BIC with the number of groups
nobs(fit)                            # number of groups

# Plots: plot(fit, what)
plot_views()                # every view and what it shows
plot(fit, what = "bars")             # profile means with 95% intervals
plot(fit, what = "heatmap")          # means in standard deviations
plot(fit, what = "sizes")            # effective size of each profile
plot(fit, what = "entropy")          # per-case classification uncertainty
plot(fit, what = "posteriors")       # posterior of each assigned profile
plot(fit, what = "avepp")            # average posteriors, assigned by posterior class
plot(fit, what = "all")              # every view this fit supports

# Other functions
sensitivity(fit, seeds = 1:10)       # refits under several seeds
starting_values(fit)                 # parameters to start another fit

The evaluation vignette covers enumerate_classes(), candidate_fit() and bootstrap_lrt(); the covariates vignette covers three_step(), r3step() and fit_staged(); the transitions vignette covers lta(), get_tna() and get_group_tna().

Sampling weights

When students were sampled unequally, for example wavering students at half the rate of committed ones, each student carries a sampling weight, and the fit should estimate the population rather than the sample. We pass the weight column as weights. The weight belongs to the independent unit, here the student, so it must be constant within id. The fit maximizes the weighted (pseudo) log likelihood, with the weights scaled to sum to the number of students, and reports Kish's effective sample size.

weighted_data <- transform(course_engagement,
                           survey_weight = ifelse(student_type == "wavering", 2, 1))
weighted_fit <- multilpa(weighted_data,
                         vars = c("browse", "lectures", "forum_read",
                                  "forum_post", "attendance"),
                         id = "student", n_profiles = 2, n_group_classes = 2,
                         weights = "survey_weight", n_starts = 10, seed = 1)
weighted_fit

Under weights the inverse information is not the variance of the estimates, so parameter_inference() reports sandwich standard errors by default and refuses vcov_type = "observed".

head(parameter_inference(weighted_fit), 4)

lta(), mixture_regression() and every family of multilpa() take weights the same way. The bootstrap likelihood ratio test does not: it simulates from the model, not from the sampling design.

All tables

To see every table at once, we call summary() on the fit, with rows for the number of rows shown per table.

print(summary(fit), rows = 3)

Assumptions

options(old_options)


Try the latents package in your browser

Any scripts or data that you put into this service are public.

latents documentation built on Oct. 16, 2026, 5:09 p.m.