---
title: "Multilevel latent profile analysis"
author: "Written by [Mohammed Saqr](https://saqr.me) and [Sonsoles López-Pernas](https://sonsoles.me)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Multilevel latent profile 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)
```

Multilevel latent profile analysis captures heterogeneity in continuous
responses by modelling variation within and across groups. Observation-level
profiles describe patterns across indicators, while latent group classes
capture differences in the prevalence of these profiles across groups.

In this vignette, course enrolments are nested within students. Engagement
profiles describe variation in activity patterns across enrolments, while
student classes capture differences in how frequently students exhibit those
patterns. The profiles share the same indicator means and covariance structure
across student classes, allowing their prevalence to be compared on a common
basis.

## Data

`course_engagement` has one row per course enrolment: 1,422 enrolments from
106 students. `multilpa()` requires a group identifier and indicator variables.

- **Group:** `student` identifies which enrolments belong to the same
  student.
- **Indicators:** five activity measures (`browse`, `lectures`, `forum_read`,
  `forum_post`, `attendance`), each standardized within course, so 0 is the
  course average.

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

```{r setup}
library(latents)
```

## Clustering within groups

The `descriptives()` function summarizes the indicators specified in `vars`
and their clustering within the groups specified by `id`. The intraclass
correlation (`icc`) estimates the proportion of an indicator's variance
attributable to between-group differences.

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

## Fit

The arguments `n_profiles` and `n_group_classes` specify the numbers of
observation-level profiles and group classes. Estimation uses multiple
expectation-maximization (EM) starts, controlled by `n_starts`; `seed`
controls reproducibility. The `cores` argument permits parallel estimation
of starts. The printed model summary includes the estimated profile means.

```{r fit}
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
```

The starting-value summary reports the log likelihood reached by each start.
Replication of the selected likelihood provides evidence of numerical
stability, although it does not establish that the global maximum has been
found.

```{r starts}
get_results(fit, "starts")
```

## Profiles

By default, `get_results()` returns profile means and variances, with one
row per profile and continuous indicator.

```{r profiles}
get_results(fit)
```

The profile plot displays the estimated indicator means for each profile.
The labels also report the profiles' estimated shares of observations.

```{r profile-plot}
plot_profiles(fit)
```

## Group classes

The `"profile_probabilities"` table reports the conditional probability of
each observation-level profile within each group class (`probability`).
The column `group_class_probability` reports the estimated proportion of
groups in that class.

```{r profile-probabilities}
get_results(fit, "profile_probabilities")
```

The prevalence plot compares these conditional profile probabilities, with
one line per group class.

```{r probability-plot}
plot_prevalence(fit)
```

## Class sizes

The `"counts"` table summarizes class sizes. The column `effective_count`
sums posterior membership probabilities and can therefore be noninteger.

```{r counts}
get_results(fit, "counts")
```

## Classification

The `"entropy"` table summarizes classification uncertainty at both levels.
For a model with more than one class, relative entropy is 1 when posterior
membership is certain and 0 when every case has equal posterior probability
of membership in every class. It describes separation under the fitted model,
not overall model fit.

```{r entropy}
get_results(fit, "entropy")
```

The `"assignments"` table appends classification results to `data` in its
original row order. The columns `profile` and `group_class` contain the most
probable classes. The column `uncertainty` is one minus the posterior
probability of the assigned profile, and the `posterior_profile_` columns
contain membership probabilities for each profile. The number of posterior
columns depends on `n_profiles`. Analyses based only on the assigned labels
omit this classification uncertainty.

```{r assignments}
get_results(fit, "assignments", data = course_engagement) |> head()
```

## Recovery of the known structure

For these simulated data, the `truth` argument permits comparison with the
generating labels: `engagement` for profiles and `student_type` for group
classes. Class numbers are arbitrary, so correspondence is assessed from the
cross-tabulation rather than by matching numeric labels directly.

```{r recovery}
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.

```{r reference, eval = FALSE}
# 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: the plot_ functions, each with what = and "all"
plot_views()                                   # every view and what it shows
plot_profiles(fit, what = "bars")              # profile means with 95% intervals
plot_profiles(fit, what = "heatmap")           # means in standard deviations
plot_profiles(fit, what = "all")               # every profile view
plot_diagnostics(fit, what = "sizes")          # effective size of each profile
plot_diagnostics(fit, what = "entropy")        # per-case classification uncertainty
plot_diagnostics(fit, what = "posteriors")     # posterior of each assigned profile
plot_diagnostics(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
```

The [class enumeration](https://pak.dynasite.org/latents/articles/enumeration.html)
article covers `enumerate_classes()` and `get_candidate()`; the
[latent transition analysis](https://pak.dynasite.org/latents/articles/lta.html)
article covers `lta()`, `get_tna()` and `get_group_tna()`. The evaluation article on the
[package website](https://pak.dynasite.org/latents/articles/) covers
`bootstrap_lrt()`, and the covariates article covers `three_step()`,
`r3step()` and `fit_staged()`.

## Sampling weights

Sampling weights can represent unequal selection probabilities. The example
below assigns twice the weight to the simulated `"wavering"` student class,
as would be appropriate if its sampling probability were half that of the
other class. The weight column is supplied through `weights` and must be
constant within the independent sampling unit, here `student`. Estimation
maximizes a weighted pseudo-log-likelihood after normalizing weights to sum
to the number of students. The output also reports Kish's effective sample
size.

```{r weights}
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
```

For weighted fits, `parameter_inference()` uses sandwich standard errors by
default. The package does not support `vcov_type = "observed"` for these
fits because inverse-information standard errors do not generally represent
the sampling variance under unequal weighting.

```{r weights-inference}
head(parameter_inference(weighted_fit), 4)
```

The `weights` argument is also available in `lta()`, `mixture_regression()`
and the model families fitted through `multilpa()`. The parametric bootstrap
likelihood ratio test does not support weighted fits; its simulations do not
reproduce a sampling design.

## All tables

The model summary collects the result tables. The `rows` argument limits
the number of rows printed from each table.

```{r summary}
print(summary(fit), rows = 3)
```

## Assumptions

- Profiles have the same means and variances in every group class.
- Observations are independent given their group's class, and groups are
  independent of one another.
- Indicators are independent within a profile (the default diagonal
  covariance; `covariance_model = "full"` relaxes it).
- The numbers of profiles and group classes are chosen by the analyst; see
  the [class enumeration](https://pak.dynasite.org/latents/articles/enumeration.html)
  article for comparing them.

```{r restore-options, include = FALSE}
options(old_options)
```