Group-Based Trajectory Models and Latent Class Analysis

In plain language

Group-based trajectory models and latent class analysis are statistical tools that look at a group of patients and automatically discover hidden subgroups — clusters of people who follow the same pattern over time or share the same combination of health conditions. Instead of summarizing everyone's medication use as one average number, these models reveal that the average may mask very different groups: for example, half the patients taking their statin consistently, 30% tapering off month by month, and 20% stopping within weeks of starting. The number of clusters is chosen by comparing how well models with two, three, or four groups fit the data, using a scoring rule called BIC. The key honesty caveat is that the groups are mathematical constructs produced by the model, not biological disease subtypes discovered in nature — two analysts with slightly different choices may identify different groups from the same patients.

What GBTM and LCA actually are

Group-based trajectory modeling (GBTM) and latent class analysis (LCA) are finite mixture methods that represent a heterogeneous population as a blend of k latent subpopulations, each with its own distinct pattern. GBTM, developed by Nagin, is specifically designed for repeated-measures longitudinal outcomes: it estimates a small number of trajectory groups — each described by a polynomial or logistic curve over time — and assigns each patient to the group whose trajectory best matches their observed sequence. LCA is the cross-sectional analogue: it identifies clusters of patients based on profiles of binary or categorical indicators measured at a single point in time (for example, presence or absence of ten comorbidity flags from claims). Both models belong to the same family — finite mixture models — and share the same core idea: the population is not homogeneous, and summarizing it with a single mean discards the heterogeneity that may be the most policy-relevant feature of the data.

The finite mixture idea

A finite mixture model with k components posits that each observation is drawn from one of k distributions, each weighted by a class-membership probability (the mixing proportion pi_k, with the sum of all pi_k equal to 1). Neither the class membership nor the class-specific trajectories are observed — both are estimated from the data simultaneously by maximum likelihood or Bayesian inference. The likelihood function sums the probability of the observed data under each class, weighted by its mixing proportion. Estimation typically uses the EM algorithm: the E-step assigns soft class memberships (posterior probabilities), and the M-step updates the class-specific parameters and mixing proportions. GBTM applies this framework to longitudinal binary or count outcomes — the canonical software, PROC TRAJ in SAS, uses a censored-normal, Poisson, zero-inflated Poisson, or logistic link depending on the outcome type. LCA applies the same framework to cross-sectional multivariate categorical data.

GBTM for medication adherence: richer than a single PDC cutoff

The Franklin et al. (2013) argument for GBTM over a single PDC scalar is direct: a PDC threshold of 0.80 collapses the entire longitudinal adherence trajectory to a pass/fail, discarding when the patient was adherent, whether they restarted after a gap, and how steeply adherence declined over time. GBTM instead identifies trajectory shapes — for example, four groups in a statin cohort: consistent adherers (high flat PDC throughout the year), slow decliners (high early but falling month by month), rapid discontinuers (high at the index date then abrupt drop at month two or three), and intermittent users (cycling on and off with repeated gaps). These groups differ not only in their 12-month PDC totals but in the pattern of exposure, which can drive very different clinical outcomes — a pattern that is invisible to the scalar. For study design, this means that conditioning a downstream outcome analysis on trajectory-group membership can recover treatment-effect heterogeneity across adherence phenotypes that a single PDC-adjusted regression cannot.

LCA for cross-sectional latent subtypes

LCA applies the same finite mixture framework to a patient's profile of binary claim flags measured at baseline — for example, ten comorbidity indicators from the Charlson index, drug-class exposure flags, and care-utilization categories. The model asks: what is the smallest number of latent patient types that can explain the observed co-occurrence patterns among these flags? A four-class solution might identify a metabolic-syndrome cluster (diabetes, hypertension, dyslipidemia co-occurring), a cardiovascular disease cluster (prior MI, stroke, heart failure), a low-burden cluster (few comorbidities), and an oncology cluster (cancer, anemia, opioid use). This LCA-derived phenotype may be a more powerful confounder-adjustment variable than any single summary score, because it captures the joint distributional structure of comorbidities rather than a weighted sum.

Model selection: BIC, entropy, minimum class size, and the reification trap

Choosing k is the most consequential analytic decision in any mixture model. The standard toolkit: (1) Bayesian Information Criterion (BIC): fit models with k = 2, 3, 4, and more classes and pick the k where BIC stops improving (lower is better; look for an elbow, not a global minimum that may be k = n). (2) Entropy: a summary of how cleanly the model separates classes (ranges 0 to 1; values at or above 0.80 indicate distinct, well- separated groups). Low entropy means posterior probabilities are diffuse — every patient has moderate probability of belonging to multiple classes — which undermines the interpretability of any class-specific estimate. (3) Minimum class size: classes containing fewer than five percent of the sample are usually not scientifically credible and are artifacts of overfitting; in a cohort of 1000, a class of 15 patients warrants skepticism. (4) Interpretability: does the extra class tell a new clinical story, or does it merely split a prior class into near-identical twins? The best-fitting model by BIC may not be the most scientifically useful.

The reification trap is the central intellectual hazard of mixture modeling. Because PROC TRAJ or poLCA will always produce a clean k-class solution with interpretable labels, it is tempting to conclude that the classes are real, discovered entities — biological subtypes, clinical phenotypes, or adherence archetypes that exist independent of the model. They are not. The classes are the model's best finite approximation to a continuous distribution of heterogeneity, chosen jointly by the EM algorithm and the analyst's choice of k, link function, and indicator set. A different dataset, a different gap definition, or a different analyst's choice of k will produce different classes. Explicit acknowledgment of this uncertainty is not optional in any rigorous manuscript: the classes are model constructs, not discovered diseases.

Posterior class probabilities and classify-analyze bias

After the EM algorithm converges, each patient has a vector of posterior probabilities — the estimated probability of belonging to each class given their observed trajectory. The tempting shortcut is modal assignment: assign each patient to their most probable class and then treat class membership as a known, perfectly measured variable in a downstream regression. This introduces classify-analyze bias: it ignores the uncertainty in class membership and typically underestimates standard errors and produces biased class-specific estimates. The theoretically correct alternative is the three-step approach (Vermunt 2010; Bakk and Vermunt 2016): first, fit the mixture model; second, assign classes by modal assignment; third, correct the downstream regression for classification error using the average posterior probabilities as a misclassification matrix. In practice, the three-step correction materially changes estimates when entropy is low (below 0.70); when entropy exceeds 0.90, modal assignment is approximately unbiased.

Growth mixture models vs GBTM

The growth mixture model (GMM) adds random effects within each class — allowing individual- level deviations around the class-average trajectory — whereas GBTM constrains each class to a fixed (deterministic) trajectory shape with no within-class residual variation beyond measurement error. GMMs are far more difficult to identify (non-convergence and class- collapse are common) and require substantially larger samples; GBTM's simpler fixed- trajectory-per-class assumption often produces more stable and reproducible results at typical RWE sample sizes of 500 to 5000 patients.

Downstream use: trajectory group as exposure and the immortal-time hazard

A common downstream use is to treat the trajectory class as an exposure variable in a Cox or logistic model — for example, asking whether rapid discontinuers have higher hospitalization risk than consistent adherers. This creates a latent immortal-time problem: the trajectory class is defined using outcome-free follow-up during the trajectory window (for example, months 1 through 12 post-index). If a patient dies or is hospitalized in month 6, they cannot complete a 12-month trajectory; survivors who complete it are inherently selected to be alive and event-free for the duration. Conditioning on trajectory- class membership therefore implicitly conditions on the future, biasing any hazard estimate for events that occur during the trajectory window. The only clean solution is to end the trajectory window before the outcome window begins — use months 1 through 6 for trajectory definition, then months 7 through 24 for outcome follow-up — ensuring that trajectory-group membership is fully determined before the outcome period starts.

Pros, cons, and trade-offs

GBTM and LCA pros: they reveal heterogeneity invisible to single-summary statistics; class labels are interpretable to clinical and health-economics audiences; they handle binary, count, and censored-normal longitudinal outcomes (GBTM) or multivariate binary cross- sectional profiles (LCA); and they integrate naturally with downstream effect-modification analyses when trajectory group is an effect modifier.

Cons: classes are model constructs — results depend on k, link function, indicator selection, and sample; the reification trap is endemic; modal assignment ignores classification uncertainty; Python lacks a dedicated GBTM implementation (GaussianMixture approximates LCA-style clustering on cross-sectional data but does not model longitudinal trajectories); convergence failures are common at high k values.

Trade-offs vs a single PDC scalar: GBTM recovers trajectory shape at the cost of model dependence; a PDC threshold of 0.80 is reproducible, regulation-endorsed, and transparent but collapses heterogeneity. Both should be reported when the research question spans how much and what pattern.

Trade-offs vs mixed-effects models: mixed-effects models treat adherence heterogeneity as a continuous random effect (each patient has their own random intercept and slope), whereas GBTM imposes a discrete class structure. Mixed-effects models are the appropriate primary model when the goal is a population-average or subject-specific trajectory mean; GBTM is the appropriate model when the goal is to identify and label discrete subpopulations with distinct behavioral phenotypes.

When to use

Use GBTM when the research question concerns the distribution of longitudinal adherence patterns (rather than a single mean), when defining adherence phenotypes as effect- modifiers or confounders in a downstream outcome regression provided the immortal-time hazard is avoided by pre-specifying non-overlapping trajectory and outcome windows, and when payer or manufacturer audiences require patient segmentation with interpretable labels. Use LCA when the goal is to identify multimorbidity phenotypes from claims or EHR comorbidity flags as a supplement to Charlson or Elixhauser summary scores, or when binary-indicator co-occurrence patterns are the target of description rather than a weighted sum.

When NOT to use

Do not use GBTM as the primary adherence endpoint in a regulatory submission or payer dossier where a reproducible PDC rate is required. Do not allow k to be chosen post-hoc to make the result look clinically meaningful; the number of classes must be pre-specified or determined by objective criteria (BIC, entropy, and minimum class size) with the final choice locked before examining downstream outcomes. Do not report low-entropy solutions (entropy below 0.70) as if classes were cleanly distinct without applying three-step corrections. Do not use GBTM as the only analysis of heterogeneity when a continuous treatment-effect modifier (for example, continuous PDC as an interaction term) would answer the same question with fewer assumptions. Critically, do not define trajectory groups over a window that overlaps the outcome window; classes defined by future data create implicit immortal time and biased hazard estimates in any downstream survival analysis.

Interpreting the output

In the worked example, a 3-class GBTM fit to 1000 statin users returns three trajectory groups: consistent adherers (n = 500, mean 12-month PDC = 0.90), moderate adherers (n = 300, mean PDC = 0.60), and rapid discontinuers (n = 200, mean PDC = 0.20).

(1) Formal interpretation. The model-implied overall population mean PDC is the class- share-weighted average: (0.5 0.9) + (0.3 0.6) + (0.2 * 0.2) = 0.45 + 0.18 + 0.04 = 0.67. This matches a naive cohort-level PDC of 0.67, but the GBTM reveals that this single number masks a trimodal distribution. The consistent-adherer class (50% of patients) drives a PDC cluster centered near 0.90, well above the conventional 0.80 adherence threshold, while the rapid-discontinuer class (20%) drives a cluster near 0.20. The moderate-adherer class (30%) sits below the 0.80 threshold at a mean of 0.60. Class membership uncertainty is quantified by the entropy statistic: if entropy is at or above 0.80, modal class assignments are approximately valid for downstream analyses; if entropy is below 0.70, three-step corrections should be applied before class membership is used as a covariate.

(2) Practical interpretation. A single average PDC of 0.67 for the cohort would suggest the population is modestly non-adherent as a group. The GBTM tells a more nuanced story for a payer or formulary decision-maker: half the cohort is highly adherent (mean PDC near 0.90) and would likely benefit from continued access to the drug; one in five patients stopped within the first few months and represents an opportunity for targeted adherence- support intervention rather than access restriction. Reporting only the 0.67 average obscures both the high-adherers who should not face access barriers and the rapid discontinuers who need a different kind of support. The trajectory classes, despite being model constructs rather than biological discoveries, provide a clinically actionable segmentation that a scalar cannot.

Worked example

Scenario

A pharmacoepidemiologist fits a 3-class GBTM to 1000 new statin users followed for 12 months after their first fill. Each patient's 12 monthly binary adherence indicators (yes/no for having statin supply on hand that month) are the input. The model returns three trajectory groups with estimated class proportions of 50%, 30%, and 20% and class-specific mean 12-month PDC values of 0.90, 0.60, and 0.20 respectively. The analyst wants to verify that the class-specific means reproduce the cohort-level mean PDC observed in the claims data.

Dataset

Class-level summary from the fitted 3-class GBTM (1000 statin users, first 12 months post-index). Each row is one trajectory class; patient counts and mean PDC are model-derived.

Steps

• Convert class shares to patient counts. Consistent adherers: 0.5 1000 = 500 patients. Moderate adherers: 0.3 1000 = 300 patients. Rapid discontinuers: 0.2 * 1000 = 200 patients. Check: 500 + 300 + 200 = 1000 (total cohort size confirmed).

• Compute each class's contribution to the overall weighted mean PDC. Consistent-adherer contribution: 0.5 0.9 = 0.45. Moderate-adherer contribution: 0.3 0.6 = 0.18. Rapid-discontinuer contribution: 0.2 * 0.2 = 0.04.

• Sum the three contributions to get the overall weighted mean PDC: 0.45 + 0.18 + 0.04 = 0.67.

• Interpretation: a naive cohort-level PDC of 0.67 is below the conventional 0.80 adherence threshold, which would suggest the cohort is non-adherent as a whole. The GBTM reveals that this single number hides a trimodal reality: 500 patients (50%) have a mean PDC near 0.90 (above the threshold), 300 patients (30%) have a mean PDC of 0.60 (below it), and 200 patients (20%) have a mean PDC of 0.20 (far below it). These three groups face very different clinical trajectories and warrant different intervention strategies.

Result

Class sizes: 0.5 1000 = 500 consistent adherers, 0.3 1000 = 300 moderate adherers, 0.2 1000 = 200 rapid discontinuers. Weighted contributions: 0.5 0.9 = 0.45, 0.3 0.6 = 0.18, 0.2 0.2 = 0.04. Overall weighted mean PDC: 0.45 + 0.18 + 0.04 = 0.67. The cohort-level mean of 0.67 obscures a trimodal adherence distribution that the three trajectory classes make explicit.

Runnable example

# ── LCA-style clustering with sklearn GaussianMixture ──
# NOTE: Approximates cross-sectional LCA on binary indicators.
# For longitudinal GBTM, use R (lcmm) or SAS (PROC TRAJ).

import numpy as np
from sklearn.mixture import GaussianMixture

# ── Worked example: verify weighted mean PDC arithmetic ──
# 1000 patients, 3 classes with shares 0.50 / 0.30 / 0.20
class_shares   = np.array([0.50, 0.30, 0.20])
class_pdc_mean = np.array([0.90, 0.60, 0.20])

# Weighted overall mean PDC
# 0.5 * 0.9 = 0.45; 0.3 * 0.6 = 0.18; 0.2 * 0.2 = 0.04
contributions = class_shares * class_pdc_mean
# contributions = [0.45, 0.18, 0.04]
overall_mean_pdc = contributions.sum()
# 0.45 + 0.18 + 0.04 = 0.67
print(f"Contributions: {contributions}")
print(f"Overall weighted mean PDC: {overall_mean_pdc:.2f}")  # 0.67

# Class patient counts from shares (n=1000)
n_total = 1000
n_per_class = class_shares * n_total  # [500, 300, 200]
print(f"Patients per class: {n_per_class.astype(int)}")

# ── LCA-style GaussianMixture on simulated comorbidity indicators ──
# 5 binary flags: diabetes, hypertension, dyslipidemia, heart_failure, prior_MI
rng = np.random.default_rng(42)
true_class = rng.choice([0, 1, 2], p=[0.50, 0.30, 0.20], size=1000)
class_probs = {
    0: [0.70, 0.80, 0.85, 0.10, 0.08],  # metabolic cluster
    1: [0.30, 0.50, 0.35, 0.65, 0.55],  # cardiovascular cluster
    2: [0.10, 0.15, 0.12, 0.05, 0.03],  # low-burden cluster
}
X = np.array([
    [rng.binomial(1, class_probs[c][j]) for j in range(5)]
    for c in true_class
], dtype=float)

# ── BIC comparison across k = 2, 3, 4 ──
print("\nBIC by k (lower = better; look for elbow):")
bic_results = {}
fitted_models = {}
for k in range(2, 5):
    gm = GaussianMixture(
        n_components=k, covariance_type="full",
        n_init=5, random_state=42
    )
    gm.fit(X)
    bic_results[k] = gm.bic(X)
    fitted_models[k] = gm
    print(f"  k={k}: BIC={bic_results[k]:.1f}")

# ── Fit k=3 and extract posterior probabilities ──
gm3 = fitted_models[3]
post_probs  = gm3.predict_proba(X)   # shape (1000, 3), rows sum to 1
modal_class = gm3.predict(X)         # modal assignment (ignores uncertainty)

# ── Entropy: classification quality ──
# Entropy = 1 + (1 / (n * ln(k))) * sum_ij [p_ij * ln(p_ij)]
eps = 1e-10
n, k = post_probs.shape
entropy = 1 + np.sum(post_probs * np.log(post_probs + eps)) / (n * np.log(k))
print(f"\nEntropy (k=3): {entropy:.3f}  (>=0.80 = good separation)")
print("If entropy < 0.70, apply 3-step correction (Bakk & Vermunt 2016) for inference.")

# ── Class mixing proportions ──
print("\nEstimated class proportions (k=3):")
for i, w in enumerate(gm3.weights_):
    n_cls = np.sum(modal_class == i)
    print(f"  Class {i}: pi={w:.3f}  ({n_cls} patients by modal assignment)")

# ── Mean simulated PDC by modal class ──
# Simulate 12-month PDC from known class-level means
pdc_true_means = {0: 0.90, 1: 0.60, 2: 0.20}
pdc = np.array([
    rng.beta(a=pdc_true_means[c] * 10, b=(1 - pdc_true_means[c]) * 10)
    for c in true_class
])
print("\nMean 12-month PDC by modal class assignment:")
for i in range(3):
    mask = modal_class == i
    print(f"  Class {i}: mean PDC = {pdc[mask].mean():.3f}")

# ── Group-Based Trajectory Modeling in R ──
# install.packages(c("lcmm", "poLCA"))
library(lcmm)
library(poLCA)

# ── 1. Simulate longitudinal adherence data (1000 patients, 12 months) ──
set.seed(42)
n <- 1000
true_cls <- sample(1:3, size = n, prob = c(0.50, 0.30, 0.20), replace = TRUE)
months    <- 1:12

sim_data <- do.call(rbind, lapply(seq_len(n), function(pid) {
  cls <- true_cls[pid]
  mu  <- c(0.90, 0.60, 0.20)[cls]
  pdc_vals <- pmin(pmax(rnorm(12, mean = mu, sd = 0.07), 0), 1)
  data.frame(id = pid, month = months, pdc = pdc_vals)
}))

# ── 2. Fit GBTM with k = 1, 2, 3 classes (hlme: hidden latent classes) ──
# hlme fits a linear mixed model with ng latent classes.
# mixture = ~month: class-specific intercept and slope on time.
# B = gridsearch(): multi-start to avoid local maxima — always use with ng >= 2.

m1 <- hlme(pdc ~ month, subject = "id", ng = 1, data = sim_data)

m2 <- hlme(pdc ~ month, subject = "id", ng = 2,
           mixture   = ~month,
           B         = gridsearch(m1, rep = 20, maxiter = 10, minit = m1),
           data      = sim_data)

m3 <- hlme(pdc ~ month, subject = "id", ng = 3,
           mixture   = ~month,
           B         = gridsearch(m1, rep = 30, maxiter = 10, minit = m1),
           data      = sim_data)

# ── 3. Compare models by BIC (lower = better) and entropy ──
cat("Model comparison — BIC and entropy:\n")
summarytable(m1, m2, m3, which = c("G", "loglik", "AIC", "BIC", "entropy"))
# Entropy >= 0.80: good class separation
# Entropy < 0.70: apply 3-step correction before downstream regression

# ── 4. Extract class assignments and posterior probabilities (k=3) ──
# m3$pprob: data.frame with columns id, class (modal), prob1, prob2, prob3
pp3 <- m3$pprob

cat("\nModal class proportions (k=3):\n")
print(prop.table(table(pp3$class)))
# Expected: ~0.50 / 0.30 / 0.20

cat(sprintf("Entropy (k=3): %.3f\n", m3$entropy))

# ── 5. Mean 12-month PDC by assigned class ──
pdcs <- aggregate(pdc ~ id, data = sim_data, FUN = mean)
pdcs$class <- pp3$class[match(pdcs$id, pp3$id)]

cat("\nMean 12-month PDC by modal trajectory class:\n")
class_means <- aggregate(pdc ~ class, data = pdcs, FUN = mean)
print(class_means)
# Expected: class means near 0.90, 0.60, 0.20 (in class-sorted order)

# ── 6. Verify worked-example weighted mean arithmetic ──
# Class shares: 0.50, 0.30, 0.20  |  PDC means: 0.90, 0.60, 0.20
shares <- c(0.50, 0.30, 0.20)
pdcs_e <- c(0.90, 0.60, 0.20)
contrib <- shares * pdcs_e          # 0.45, 0.18, 0.04
overall <- sum(contrib)             # 0.67
cat(sprintf("\nWeighted mean PDC: %.2f  (expected 0.67)\n", overall))

# ── 7. Cross-sectional LCA with poLCA ──
# poLCA requires indicators coded as 1, 2, ... (not 0/1)
set.seed(42)
sim_lca <- data.frame(
  diabetes     = sample(1:2, 500, replace = TRUE, prob = c(0.40, 0.60)),
  hypertension = sample(1:2, 500, replace = TRUE, prob = c(0.50, 0.50)),
  dyslipidemia = sample(1:2, 500, replace = TRUE, prob = c(0.45, 0.55)),
  hf           = sample(1:2, 500, replace = TRUE, prob = c(0.15, 0.85)),
  prior_mi     = sample(1:2, 500, replace = TRUE, prob = c(0.12, 0.88))
)
f_lca <- cbind(diabetes, hypertension, dyslipidemia, hf, prior_mi) ~ 1

lca_bic <- sapply(2:4, function(k) {
  fit <- poLCA(f_lca, data = sim_lca, nclass = k, verbose = FALSE, nrep = 5)
  fit$bic
})
names(lca_bic) <- paste0("k=", 2:4)
cat("\nLCA BIC by k (lower = better):\n")
print(lca_bic)

best_k <- which.min(lca_bic) + 1
lca_final <- poLCA(f_lca, data = sim_lca,
                   nclass = best_k, verbose = FALSE, nrep = 10)
cat(sprintf("\nSelected LCA k=%d. Class mixing proportions:\n", best_k))
print(lca_final$P)
# Class-conditional item probabilities reveal the comorbidity profile of each class