--- title: "Two-level latent class analysis" author: "Written by [Mohammed Saqr](https://saqr.me) and [Sonsoles López-Pernas](https://sonsoles.me)" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Two-level latent class analysis} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r options, include = FALSE} 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) represents heterogeneity in categorical response patterns among observations nested within groups. Observation-level classes, labelled profiles in the output, are defined by item-response probabilities. Group classes differ in their conditional probabilities of observation-level class membership. By default, `multilca()` treats all selected indicators as categorical. `multilpa()` supports mixed continuous and categorical indicators through the `categorical` argument. ## What the package supports - **Item types.** Binary, ordinal and nominal items share one parameterization: unrestricted category probabilities in each profile, subject to summing to one. 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 shared category intercepts), 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 `bootstrap_lrt()`, which keeps each row's pattern of missing answers in a FIML fit. - **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. ```{r setup} library(latents) ``` The example uses the first week of observations (`day <= 6`), retaining all 100 students while limiting computation. ```{r first-week} first_week <- subset(student_esm, day <= 6) ``` The eight activity items are used in several fits, so we store their names once. ```{r activities} activities <- c("time_with_friends", "on_social_media", "tv_video_games", "listened_music", "sports", "walking", "reading", "part_time_job") ``` ## Fitting the model The model specifies two observation-level classes and two group classes. The `tol` argument sets a stricter convergence tolerance for this example; accurate numerical standard errors require a sufficiently well-converged solution. ```{r lca} lca <- multilca(first_week, vars = activities, id = "student", n_profiles = 2, n_group_classes = 2, n_starts = 3, tol = 1e-10, seed = 1) lca ``` The `"responses"` table contains one row per profile, item and category. The column `probability` is the estimated conditional probability of that response within the profile. ```{r responses} get_results(lca, "responses") ``` The response plot compares the fitted item-response probabilities across profiles. ```{r responses-plot} plot_responses(lca) ``` For numeric categories, such as ratings from 1 to 5, a profile plot can summarize each item's fitted distribution on the original score scale: ```{r categorical-summaries, eval = FALSE} plot_profiles(lca) # probability-weighted means plot_profiles(lca, statistic = "median") # fitted medians plot_profiles(lca, statistic = "mode") # most probable scores ``` Means assume meaningful spacing between category scores. The median is the smallest score whose cumulative probability reaches 0.5; tied modes use the smallest score. These summaries use fitted probabilities, including in FIML models, rather than averages of cases assigned to each class. They require numeric category labels and use the raw scale. Confidence intervals are not computed, so these plots avoid the parameter-inference calculation used by the response-probability plot. In mixed models, `plot_profiles()` continues to show only continuous indicators. The `"profile_probabilities"` table and the prevalence plot summarize how the distribution of observation-level classes differs across group classes. ```{r profile-probabilities} get_results(lca, "profile_probabilities") ``` ```{r probability-plot} plot_prevalence(lca) ``` The `"entropy"` table summarizes posterior classification separation at both levels. ```{r entropy} get_results(lca, "entropy") ``` ## Mixed indicators Mixed models specify all indicators in `vars` and the categorical subset in `categorical`. This example treats three affect ratings as continuous and the activity indicators as categorical. Treating ratings as continuous is a measurement assumption that should be evaluated against an ordinal specification. The `worried` item is omitted because its concentration at the lowest category can produce a component with variance near the lower bound under a Gaussian specification. ```{r mixed} 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) ``` The default results table reports the continuous indicators' means and variances. The `"responses"` table reports categorical response probabilities. ```{r mixed-results} get_results(mixed) get_results(mixed, "responses") ``` ## Ordinal and count indicators The affect ratings range from 1 to 7. An unrestricted categorical specification requires six free probability parameters per item and profile. The ordinal specification instead uses shared category intercepts and profile locations under an adjacent-category logit model, with identifying constraints described in `?multilpa`. The `ordinal` argument identifies the items to which this specification applies. ```{r 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 ``` The `"ordinal"` table reports fitted category probabilities and profile locations for the ordinal indicators. ```{r ordinal-results} head(get_results(ordinal, "ordinal"), 7) ``` Count indicators, such as numbers of messages or absences, are specified in `count`. The default Poisson model estimates one mean per profile, available through `get_results(fit, "count_means")`. The negative-binomial option, `count_model = "negative_binomial"`, permits within-profile overdispersion through profile-specific dispersion parameters, or a common parameter with `count_dispersion = "equal"`. Both count specifications support membership covariates and latent transition models. ## Reference The calls below list other functions that work with categorical and mixed models. Each comment states what the call returns. ```{r reference, eval = FALSE} 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_profiles(mixed) # means of the continuous indicators plot_responses(mixed) # 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. - Under unrestricted categorical measurement, each item has a separate response-probability distribution in every profile. Rare responses can produce estimates at the lower bound (`min_probability`). - Categories are taken from the observed values; factor levels keep their order. ```{r restore-options, include = FALSE} options(old_options) ```