multilpa: Fit a two-level latent profile model

In latents: Multilevel Latent Profile, Class and Transition Models

Fit a two-level latent profile model

Description

Fits individual Gaussian profiles nested within observed groups. A discrete latent group class determines the profile proportions. Profile means and variances are shared across group classes (measurement invariance), and indicators are independent conditional on individual profile membership by default. Full residual covariance and observed-data maximum likelihood for missing indicators are available.

Usage

multilpa(
  data,
  vars,
  id,
  n_profiles,
  n_group_classes = 2L,
  profile_covariates = character(),
  group_covariates = character(),
  variance_model = c("varying", "equal"),
  n_starts = 10L,
  max_iter = 1000L,
  tol = 1e-08,
  min_variance = 1e-06,
  seed = NULL,
  start = NULL,
  missing = c("error", "fiml"),
  profile_slopes = c("shared", "group_class"),
  covariance_model = c("diagonal", "full"),
  categorical = character(),
  min_probability = 1e-10,
  time = NULL,
  fixed = character(),
  centering = c("none", "person", "grand"),
  volume = NULL,
  shape = NULL,
  orientation = NULL,
  model = NULL,
  select_start = c("likelihood", "converged"),
  prior = NULL,
  noise = FALSE,
  acceleration = c("squarem", "none"),
  family = c("profiles", "additive", "dispersion", "additive_dispersion",
    "restricted_cross_level", "full_cross_level"),
  between_variance = c("varying", "equal"),
  weights = NULL,
  ordinal = character(),
  count = character(),
  count_model = c("poisson", "negative_binomial"),
  count_dispersion = c("varying", "equal")
)

Arguments

data	A data frame containing indicators and a group identifier.
vars	Unique character vector of continuous indicator column names.
id	Name of the observed group identifier column. Character, factor, or numeric identifiers are supported; missing identifiers are not. id has no default: omitting it raises latents_bad_argument, because a forgotten grouping would otherwise be fitted as a different model without saying so. Passing id = NULL explicitly fits a single-level model: the observations are treated as independent, each row is its own unit, and n_group_classes becomes one. That fit raises a latents_single_level message saying so. This is the ordinary Gaussian or latent-class mixture that the two-level model reduces to, and every verb of this package works on it. The unit column is fabricated internally as .observation; it is not returned by get_results(x, "data"), and a data that already has a column of that name raises latents_bad_data. Asking for more than one group class without an id raises latents_bad_argument, because one observation per unit leaves no composition for a second-level class to differ in.
n_profiles	Positive integer number of individual profiles.
n_group_classes	Positive integer number of latent group classes. With id = NULL it is one, and naming anything else is an error.
profile_covariates, group_covariates	Names of numeric columns of data predicting individual-profile and group-class membership through multinomial logits, with the final class as reference. Naming either one fits the one-step covariate model and returns a multilpa_covariates object: profile intercepts differ by group class, profile slopes are shared across group classes unless profile_slopes = "group_class", and a group_covariate must be constant within each group. Covariates enter in their supplied units, so centre or scale them beforehand if that is what you want. This is one-step maximum likelihood, not a regression on assigned classes; three_step() and r3step() are the staged alternatives. The covariate path supports neither start nor fixed, and refuses them by name rather than ignoring them; it supports missing = "fiml" for the indicators, but covariates must be complete. Use one for an ordinary, pooled LPA.
variance_model	Either "varying" (profile-specific variances or covariance matrices) or "equal" (shared across profiles).
n_starts	Positive integer number of EM starts. The first is deterministic: Ward's hierarchical clustering of the standardized indicators (on at most 2000 evenly spaced rows), each cluster described by its own means and spread, which is the start mclust's name refers to and finds solutions whose profiles differ in shape. The rest are k-means from random centres. When start is supplied, it supplies the first start instead; remaining starts are random initializations. It is ignored when max_iter = 0 and start is supplied: that call evaluates the supplied parameters and nothing else, so exactly one start is run whatever n_starts says.
max_iter	Nonnegative integer maximum number of EM updates per start. max_iter = 0 performs no update. With start, it returns the model evaluated at exactly those values, whatever n_starts is, so that logLik() scores a parameter set supplied from elsewhere rather than the best of some random initializations that were never asked for. Without start there is nothing to evaluate at, so the random initializations are scored and the highest is returned.
tol	Positive relative log-likelihood tolerance. Convergence requires absolute change no greater than ⁠tol * (1 + abs(previous log likelihood))⁠.
min_variance	Positive lower bound on each variance, or each covariance eigenvalue for full covariance, in squared input units. This defines a constrained maximum-likelihood problem. Bound-active estimates are explicitly reported and generate a warning.
seed	Optional random seed: any whole number set.seed() accepts, negative ones included. With a supplied seed, the caller's random-number state is restored on exit.
start	Optional list of means and variances (profiles by indicators), profile_probabilities (group classes by profiles), and group_probabilities (vector). Starting probabilities must be positive. For full covariance, supply covariances (indicators by indicators by profiles); variances may be omitted or must match their diagonals. A categorical model also takes response_probabilities, a list of one profiles-by-categories matrix per indicator named in categorical. That list is read by position unless it is labelled: name its elements after the indicators, or its columns after the categories, and each block is matched to the indicator and category its labels name, so a differently ordered categorical or a differently ordered set of factor levels cannot attach a distribution to the wrong item. A label naming an indicator or a category this fit does not have raises latents_bad_start rather than being aligned by position.
missing	"error" rejects missing indicators; "fiml" maximizes the observed-data likelihood under an ignorable missingness mechanism (MAR). Missing indicators are integrated out, not filled in for likelihood fitting. This holds with membership covariates too; the covariates themselves must be complete (latents_bad_data).
profile_slopes	"shared" (the default) gives each profile covariate one slope per profile logit, common to every group class. "group_class" lets each group class have its own slopes, so a covariate can predict profile membership differently in different kinds of group (a cross-level interaction); its coefficients are reported with terms such as z:group_class_1. It adds (n_group_classes - 1) times as many slope parameters again, each estimated from the groups in its class, so it needs enough groups per class. Naming it without profile_covariates raises latents_bad_argument.
covariance_model	"diagonal" assumes conditional independence; "full" estimates within-profile residual covariances.
categorical	Character vector naming indicators to treat as categorical. Each is modelled by unrestricted, profile-specific response probabilities over its observed categories, which is the latent class measurement model. Binary, ordinal and unordered indicators are all handled by the same unrestricted parameterization; numeric, integer, logical, character and factor columns are accepted. Indicators not named here stay Gaussian, so naming a subset fits a mixed-mode model.
min_probability	Positive lower bound on every categorical response probability, defining a constrained maximum-likelihood problem in the same way min_variance does for Gaussian indicators.
time	Optional name of a column giving each observation's position within its group, such as a wave, occasion or course number. The model does not use it; it is stored so that get_results(x, "sequences"), get_results(x, "sequence_summary") and plot(what = "sequences") can read the assignments back in order. Values must be complete and unique within each group.
fixed	Character vector naming measurement blocks to hold at the values start supplies, instead of estimating them: any of "means", "variances" and "response_probabilities", or "measurement" for every block the model has. Under covariance_model = "full", "variances" holds the residual covariance matrices. A held block stays exactly as supplied, in every restart, and stops counting towards n_parameters, so this is a different model rather than a different starting point for the same one. start must carry the named blocks; starting_values(fit, what = "measurement") produces them, and fit_staged() wraps the whole two-stage workflow in one call. The mixing parameters cannot be held: they are what a fixed-measurement fit is for.
centering	How to centre the continuous indicators before fitting. "none", the default, fits them as supplied. "person" subtracts each group's own mean from its rows, so a value reads as a deviation from that unit's average and the profiles become profiles of change rather than of level: this is the within-person, person-mean-centred design (Quintana, 2021; Voelkle, Brose, Schmiedek, & Lindenberger, 2014). "grand" subtracts one mean per indicator, which moves the origin without touching the within-group structure. The offsets are kept on the fit, so get_results(x, "data") still returns the columns you supplied and every verb that checks row alignment still checks it. Centring removes exactly the between-unit variation, so "person" refuses with latents_bad_data when it leaves an indicator constant — which is what happens when a unit has one observation of it. With "person" the group classes become types of change pattern, not types of unit.
volume, shape, orientation	The covariance structure, in the three pieces it is made of. Each profile's covariance decomposes as ⁠Sigma_k = lambda_k * D_k * A_k * D_k'⁠: a volume ⁠lambda_k = |Sigma_k|^(1/d)⁠, an orientation D_k of eigenvectors, and a shape A_k, diagonal with determinant one. Constraining the three across profiles gives the fourteen models mclust names with three letters, and this package fits all of them. volume is "equal" or "varying". shape is "equal", "varying" or "spherical", the last making every indicator's spread equal within a profile, which leaves no orientation to constrain. orientation is "axis" (axis-parallel, a diagonal covariance), "equal" (one orientation shared by every profile) or "varying". Each is NULL by default, which follows variance_model and covariance_model, so a call that names none of them fits exactly what it always did. volume shape orientation model equal spherical --- EII varying spherical --- VII equal/varying equal/varying axis EEI, VEI, EVI, VVI equal/varying equal/varying equal EEE, VEE, EVE, VVE equal/varying equal/varying varying EEV, VEV, EVV, VVV The estimates are those of Celeux and Govaert (1995); the two models with a shared orientation and a free shape, EVE and VVE, have no closed form and use the minorize-maximize step of Browne and McNicholas (2014). Parameter counts match mclust's own for all fourteen. Anything other than EEI, VVI, EEE or VVV is maximized across every profile at once, so it cannot be combined with a held variances block. parameter_inference() reports Wald standard errors for all fourteen: a constrained structure is differentiated in its own free coordinates (log volumes, determinant-one log shapes and orientations), which have exactly as many dimensions as the structure has parameters, and carried to the reported variances and covariances by the delta method. method = "bootstrap" is available for all fourteen as well.
model	The covariance structure named by its mclust code, one of "EII", "VII", "EEI", "VEI", "EVI", "VVI", "EEE", "VEE", "EVE", "VVE", "EEV", "VEV", "EVV", "VVV": the first letter is the volume, the second the shape and the third the orientation, each equal (E) or varying (V) across profiles, with I for axis-aligned (no covariances). The same codes enumerate_classes() takes. Use it instead of variance_model/covariance_model or volume/shape/orientation, not together with them.
select_start	Which start the fit reports. "likelihood", the default, takes the highest log likelihood across starts, preferring a converged start only among those tied with it to within a relative 1e-10. "converged" takes the highest log likelihood among the converged starts whenever at least one converged, and falls back to every start when none did. Use it when an unconverged start edges out a converged one by a negligible amount, which happens when a membership logit creeps along a flat ridge of the likelihood: the converged start is then a maximum that parameter_inference() can work with, and the other is not. summary() and get_results(fit, "starts") show every start either way.
prior	NULL, the default, for maximum likelihood, or prior_control() for maximum a posteriori estimation of the Gaussian means and covariances under the conjugate prior of Fraley and Raftery (2007) — mclust's priorControl(). The hyperparameters default to mclust's defaultPrior() computed from the indicators being fitted, and each M-step is mclust's own, so a fit from the same start reproduces mclust::Mclust(..., prior = priorControl()). It is defined for EII, VII, EEI, VEI, EVI, VVI, EEE, EEV, VEV and VVV (mclust has none for VEE, EVE, VVE and EVV), for complete continuous indicators without fixed or membership covariates; anything else raises latents_unsupported_prior. The mixing proportions are not given a prior. Following mclust, the fit's log_likelihood is the unpenalized log likelihood evaluated at the posterior mode, and aic, bic and their variants are computed from it with the usual parameter count; EM convergence is judged on its relative change. A prior's MAP step can lower that likelihood, so the decrease check made under maximum likelihood is not applied. parameter_inference(method = "wald") is refused for such a fit (latents_unsupported_inference), because the likelihood's score is not zero at a posterior mode; method = "bootstrap" refits every resample with the same prior. The fit records the request as prior and the resolved hyperparameters, in the indicators' units, as prior_parameters.
noise	FALSE, the default, or TRUE to add a noise component: one more mixture component with a constant density 1 / V over the hypervolume V of the data, which absorbs observations no Gaussian profile explains — outliers and scatter — instead of letting them distort a profile (Banfield & Raftery, 1993; mclust's initialization = list(noise = ) and Vinv). V is mclust's hypvol(): the smaller of the volumes of the data's axis-aligned and principal-component-aligned bounding boxes. The component has no measurement parameters; its mixing proportion is estimated with the others, and, as mclust counts it, the model has two more parameters (the proportion and the hypervolume), so its criteria equal mclust's. A random start puts a tenth of the mass on the noise; a start supplies noise_probability beside profile_probabilities, the two summing to one. Available for continuous, complete indicators with one group class (id = NULL, or n_group_classes = 1) and any of the fourteen covariance structures, alone or with prior; otherwise latents_unsupported_noise. In the fit, the noise component is profile 0, as in mclust's classification: subject_profiles is 0 for an observation assigned to it, get_results(fit, "assignments") adds a posterior_noise column, get_results(fit, "posteriors") and get_results(fit, "profile_probabilities") have profile-0 rows, and the classification diagnostics count it as a class. means, variances, covariances, profile_probabilities and subject_posteriors describe the Gaussian profiles only, so the latter two sum to one minus the noise share; the fit adds noise_probability, noise_posteriors, n_noise and hypervolume. Verbs that do not yet account for the component — parameter_inference(), three_step(), r3step(), bootstrap_lrt(), starting_values(), bivariate residuals, fit_staged() and the posterior plots — refuse such a fit with latents_unsupported_noise.
acceleration	"squarem", the default, accelerates EM with SQUAREM (Varadhan & Roland, 2008): each cycle extrapolates along the last two EM steps and finishes with an ordinary EM step, and falls back to plain EM whenever the extrapolated point would lower the likelihood, so the iteration stays monotone and converges to the same kind of maximum. The convergence test is read over a whole cycle, so it is never looser than plain EM's. "none" runs plain EM, which follows the same path as other EM implementations from the same start (mclust, Mplus) and is the choice for step-by-step reproduction. iterations and max_iter count EM steps either way. Fits with prior, and models with membership covariates, always use plain EM. For EVE and VVE, whose shared orientation has no closed form, each EM iteration takes a fixed number of warm-started orientation steps rather than solving it (a generalized EM step), under either setting; this keeps every iteration monotone and makes the EM map a fixed function, which SQUAREM requires.
family	The model family. "profiles", the default, is the model described above: individual profiles whose prevalences differ across group classes. The other three are group-class-only families with no individual profiles (Houle, Morin & Harvey, 2026): each group has a Gaussian intercept per indicator, drawn from its group class's distribution, and its members vary around it. "additive": group classes differ in their means (and, optionally, in between-group variances); one within-group variance per indicator is shared by every class. "dispersion": group classes differ in their within-group variances; means and between-group variances are shared. "additive_dispersion": group classes differ in means and in within-group variances. These families take only data, vars, id, n_group_classes, between_variance, n_starts, max_iter, tol, min_variance and seed; any other argument is refused with latents_bad_argument. Indicators must be complete and continuous, and at least one group must have more than one row (latents_unidentified otherwise). Between variances estimated at zero are reached by a boundary maximization with a Karush-Kuhn-Tucker check and flagged (latents_boundary); interior fits are finished by Newton steps. The result has class multilpa_additive; bic uses the number of groups and bic_individual the number of rows. Read it with get_results.multilpa_additive(); summary(), plot(), coef(), vcov(), confint() and parameter_inference() work on it. These families are experimental: enumeration over families, bootstrap inference, missing data and covariates are not yet available. See the group-class and cross-level models article. Two cross-level families estimate individual profiles and group classes from the same indicators, following the manifest-aggregation specification of Houle et al. (2026, supplement): each group's means of the indicators are between-level indicators of its group class, with class-specific means and variances. "restricted_cross_level": profile prevalences are the same in every group class, so profiles and group classes are estimated separately; the profile composition of each class is reported descriptively. "full_cross_level": profile prevalences differ across group classes (as in "profiles"), and the group means also inform the classes. These take data, vars, id, n_profiles, n_group_classes, variance_model (individual profiles), between_variance (group classes), n_starts, max_iter, tol, min_variance and seed. The group means are computed from the same ratings, so the likelihood is the specification's working likelihood: comparable between cross-level fits of the same data, not with the other families. Standard errors are not available (latents_unsupported_inference). The result has class multilpa_cross_level, read with get_results.multilpa_cross_level().
between_variance	For the additive and additive-dispersion families: "varying" (the default) estimates between-group variances per group class, "equal" one set shared by all classes. The dispersion family holds them equal; "varying" is refused there. For the cross-level families, the variances of the group means within each group class.
weights	NULL, or the name of a numeric column of data holding a sampling weight for each independent unit: each id group of a two-level fit, each row of a single-level one. A two-level weight must be constant within its group (latents_bad_weights otherwise); within-unit weights are not supported. The fit maximizes the pseudo log likelihood ⁠sum_j w_j log L_j⁠ (Skinner, 1989) with the weights scaled to sum to the number of units, as Mplus does, so integer weights give the fit to the data with each unit repeated that many times and the criteria stay on the sample's scale. Standard errors are the sandwich: parameter_inference() defaults to vcov_type = "robust" and refuses "observed" and "opg" (latents_unsupported_weights), and method = "bootstrap" resamples the units with their weights. Classification tables report each unit's own posterior; effective counts and proportions are weighted. bootstrap_lrt(), three_step(), r3step(), prior and noise refuse a weighted fit. Works with every family, membership covariates, categorical indicators and missing = "fiml".
ordinal	Names of indicators in vars to model as ordinal: an (ordered) factor, or whole numbers, whose categories are ordered by factor level or by value. Each follows an adjacent-category logit with category intercepts shared by every profile and one location per profile, ⁠log P(k | c) / P(k - 1 | c) = a_k - a_(k-1) + eta_c⁠, with the last profile's location fixed at zero: (K - 1) + (C - 1) parameters per indicator against C (K - 1) for a categorical one. This is Latent GOLD's default ordinal model, against which it is checked. Read the category probabilities and locations with get_results(fit, "ordinal").
count	Names of indicators in vars holding non-negative whole numbers, modelled as Poisson with one mean per profile (or negative binomial, see count_model); read the means with get_results(fit, "count_means"). Ordinal and count indicators take missing = "fiml", weights, two-level fits, membership covariates, Wald and bootstrap inference, bootstrap_lrt(), predict() and lta(); start, fixed, prior, noise and the group-class families do not take them yet (latents_unsupported_indicator).
count_model	"poisson" (the default) or "negative_binomial": NB2, with mean mu and variance ⁠mu + alpha mu^2⁠, for counts more variable than a Poisson within a profile (Latent GOLD's ⁠poisson overdispersed⁠). A dispersion estimated at zero is the Poisson limit, a boundary fit (latents_boundary warning; Wald inference is refused there).
count_dispersion	For count_model = "negative_binomial": one dispersion per profile ("varying", the default) or one shared by every profile ("equal").

Details

Infinite and constant observed indicators are rejected. Each indicator must have at least two distinct observed values. No rows are silently dropped. In FIML mode, fully missing individuals contribute no direct measurement likelihood but receive posterior probabilities from their group's information. Observation counts retain these individuals; the individual-level BIC excludes them through n_informative, while the group-level BIC uses all observed groups. Missing patterns can prevent parameter identification; no general identification guarantee is made. Initialization alone uses indicator-mean filling. EM uses conditional Gaussian sufficient statistics and observed marginal densities. More than one group class requires more than one profile and at least one group with multiple individuals; these checks are necessary but do not establish identification. By default the highest finite likelihood across starts is returned, even if that start did not converge; select_start = "converged" prefers a converged start, and converged and starts show which was chosen. Profile and group-class labels are arbitrary.

Value

An multilpa object containing means, variances, optional covariances (indicators by indicators by profiles), profile_probabilities, group_probabilities, posterior matrices, classifications, log likelihood, information criteria, restart diagnostics, and convergence history. Individual rows retain their input order; groups retain first-occurrence order. bic and bic_groups use the observed group count; bic_individual uses n_informative, the number of rows carrying at least one observed indicator. These are alternative conventions, not interchangeable criteria. boundary identifies variance bounds; small_classes flags effective memberships below one. The fit itself carries no standard errors: parameter_inference() computes them from the fit and the data it was fitted to, and bootstrap_lrt() tests nested models. No guarantee of global optimality is given, whatever n_starts is used. Read the tidy form with as.data.frame(); what selects which table.

References

Vermunt, J. K. (2003). Multilevel latent class models. Sociological Methodology, 33, 213–239. doi:10.1111/j.0081-1750.2003.t01-1-00131.x.

Banfield, J. D., & Raftery, A. E. (1993). Model-based Gaussian and non-Gaussian clustering. Biometrics, 49, 803–821. doi:10.2307/2532201.

Houle, S. A., Morin, A. J. S., & Harvey, J.-F. (2026). Multilevel latent profile analyses: A comprehensive guide. Organizational Research Methods. doi:10.1177/10944281261469432.

Skinner, C. J. (1989). Domain means, regression and multivariate analysis. In C. J. Skinner, D. Holt, & T. M. F. Smith (Eds.), Analysis of complex surveys (pp. 59–87). Wiley.

Asparouhov, T. (2005). Sampling weights in latent variable modeling. Structural Equation Modeling, 12, 411–434. doi:10.1207/s15328007sem1203_4.

Fraley, C., & Raftery, A. E. (2007). Bayesian regularization for normal mixture estimation and model-based clustering. Journal of Classification, 24, 155–181. doi:10.1007/s00357-007-0004-5.

See Also

multilca() for a model in which every indicator is categorical.

Examples

# Course sessions nested in students: engagement profiles of sessions, and
# classes of students that differ in how often they are engaged.
activity <- c("browse", "lectures", "forum_read", "forum_post", "attendance")
fit <- multilpa(course_engagement, activity, "student",
                n_profiles = 2, n_group_classes = 2, n_starts = 4,
                seed = 42)
summary(fit)
as.data.frame(fit)
get_results(fit, what = "profile_probabilities")

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