Latent Transition Analysis¶

Latent classes that move — LCA in time, via forward-filter backward-sample¶

Everything so far has been a snapshot: each subject sits in one latent class. Latent transition analysis (LTA; Collins & Lanza 2010) puts the classes in motion — each subject occupies a class at every wave, and the class evolves as a Markov chain. It is a per-subject hidden Markov model with parameters shared across the sample:

$$S_{i,1}\sim\text{Categorical}(\pi),\qquad S_{i,t}\mid S_{i,t-1}=a\sim\text{Categorical}(\tau_{a,\cdot}),\qquad x_{i,t,j}\mid S_{i,t}=c\sim\text{Categorical}(\delta_{c,j,\cdot}).$$

$\pi$ is the class mix at wave 1; $\tau$ is the $C\times C$ transition matrix ($\tau_{ab}=\Pr(\text{class }a\to\text{class }b)$); $\delta$ is the familiar LCA measurement model, held invariant over time so a "class" means the same thing at every wave. Ordinary LCA is the case $T=1$; the new object is $\tau$, which tells us how people move — who stays, who escalates, who recovers.

From scratch the sampler is a data-augmentation Gibbs whose engine is forward-filter backward-sample (FFBS) — draw each subject's whole state sequence $S_{i,1:T}$ jointly from its exact conditional, exactly the latent-state sampler behind Markov-switching models. Given the sequences, $\pi$, $\tau$ and $\delta$ are all conjugate Dirichlet. We validate on a simulated cohort, then fit the classic National Youth Survey marijuana data (237 adolescents, 5 annual waves) and read off the escalation dynamics. A PyMC marginalised-HMM fit confirms it.

Data: marijuana.csv — 237 respondents, 5 waves, each self-reporting marijuana use as 1 = none, 2 = occasional, 3 = frequent (from the LMest data_drug set; Elliott's National Youth Survey).

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import lta as M
rng = np.random.default_rng(3)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
np.set_printoptions(precision=3, suppress=True)
print("Latent transition analysis = a shared hidden Markov model over latent classes.")

Latent transition analysis = a shared hidden Markov model over latent classes.

1. Does it work? — recovering a known transition matrix¶

A simulated cohort of 800 subjects over $T=3$ waves, two classes (low-risk / high-risk) each measured by 4 binary items. The truth has 20% of low-risk subjects escalating each wave and high-risk being "sticky" (90% stay). We check that FFBS recovers $\pi$, the transition matrix $\tau$, and the item profiles $\delta$.

In [2]:
C,J,L,T = 2,4,2,3
delta_t = np.zeros((C,J,L))
for j in range(J):
    delta_t[0,j]=[0.85,0.15]        # low-risk: mostly the '0' response
    delta_t[1,j]=[0.20,0.80]        # high-risk: mostly the '1' response
pi_t  = np.array([0.70,0.30])
tau_t = np.array([[0.80,0.20],[0.10,0.90]])
Xs,Ss = M.simulate_lta(800,T,pi_t,tau_t,delta_t,rng)
sim = M.lta_gibbs(Xs, C, rng, draws=2000, burn=2000)
pi_e, tau_e = sim['pi'].mean(0), sim['tau'].mean(0)
print("pi   true", pi_t, " est", pi_e)
print("tau  true\n", tau_t, "\ntau  est\n", tau_e, f"\n(max abs error {np.abs(tau_e-tau_t).max():.3f})")

fig,ax=plt.subplots(1,2,figsize=(11,4.2))
for a,(Mx,ttl) in zip(ax,[(tau_t,"true $\\tau$"),(tau_e,"estimated $\\tau$")]):
    im=a.imshow(Mx,cmap="Blues",vmin=0,vmax=1)
    for r in range(C):
        for cc in range(C): a.text(cc,r,f"{Mx[r,cc]:.2f}",ha="center",va="center",color="k")
    a.set_xticks(range(C)); a.set_yticks(range(C)); a.set_xticklabels(["low","high"]); a.set_yticklabels(["low","high"])
    a.set_xlabel("to class"); a.set_ylabel("from class"); a.set_title(ttl)
plt.tight_layout(); plt.show()
print("FFBS recovers the transition matrix (and pi, delta): the sampler reconstructs not just the classes")
print("but the DYNAMICS -- who stays low-risk, who escalates, how sticky the high-risk state is.")

pi   true [0.7 0.3]  est [0.703 0.297]
tau  true
 [[0.8 0.2]
 [0.1 0.9]] 
tau  est
 [[0.801 0.199]
 [0.128 0.872]] 
(max abs error 0.028)
No description has been provided for this image 
FFBS recovers the transition matrix (and pi, delta): the sampler reconstructs not just the classes
but the DYNAMICS -- who stays low-risk, who escalates, how sticky the high-risk state is.

2. The marijuana data — how adolescent drug use unfolds¶

237 adolescents report marijuana use over 5 annual waves (none / occasional / frequent). A single indicator per wave, so the latent state is a de-noised use-status with a misclassification (measurement) distribution, and $\tau$ tells us how status changes year to year. We fit $C=3$ latent states, ordered low→high use.

In [3]:
seqs = pd.read_csv("marijuana.csv").to_numpy() - 1     # 0=none,1=occasional,2=frequent
X = seqs[:, :, None]                                    # (237,5,1)
print(f"{X.shape[0]} adolescents, T={X.shape[1]} waves, single ordinal indicator (3 levels)")
fit = M.lta_gibbs(X, 3, rng, draws=3000, burn=3000)
pi, tau, delta = fit['pi'].mean(0), fit['tau'].mean(0), fit['delta'].mean(0)
lab=["none","occasional","frequent"]
print("\ninitial mix pi (wave 1):", dict(zip(lab, pi.round(2))))
print("transition matrix tau (row=from, col=to):\n", tau)
print("emission delta (latent state -> observed level):\n", delta[:,0,:])

# implied class prevalence over the 5 waves (pi propagated through tau)
prev = np.zeros((5,3)); prev[0]=pi
for t in range(1,5): prev[t]=prev[t-1]@tau
fig,ax=plt.subplots(1,2,figsize=(12.5,4.4))
im=ax[0].imshow(tau,cmap="Oranges",vmin=0,vmax=1)
for r in range(3):
    for c in range(3): ax[0].text(c,r,f"{tau[r,c]:.2f}",ha="center",va="center")
ax[0].set_xticks(range(3)); ax[0].set_yticks(range(3)); ax[0].set_xticklabels(lab); ax[0].set_yticklabels(lab)
ax[0].set_xlabel("to state next wave"); ax[0].set_ylabel("from state"); ax[0].set_title("Year-to-year transition matrix $\\tau$")
cols=[GREEN,ORANGE,RED]
for c in range(3): ax[1].plot(range(1,6), prev[:,c], 'o-', color=cols[c], lw=2.2, label=lab[c])
ax[1].set_xlabel("wave (annual)"); ax[1].set_ylabel("share of cohort"); ax[1].set_ylim(0,1)
ax[1].set_title("Latent use-state prevalence across waves"); ax[1].legend(frameon=False)
plt.tight_layout(); plt.show()
print(f"\nPersistence + escalation: non-users stay non-users {tau[0,0]:.0%} of the time but {tau[0,1]:.0%} drift to occasional;")
print(f"occasional users are the volatile middle ({tau[1,2]:.0%} escalate to frequent); frequent use is sticky ({tau[2,2]:.0%} stay).")
print("The cohort steadily shifts out of non-use across the five years -- the hallmark adolescent escalation curve.")

237 adolescents, T=5 waves, single ordinal indicator (3 levels)

initial mix pi (wave 1): {'none': np.float64(0.91), 'occasional': np.float64(0.07), 'frequent': np.float64(0.02)}
transition matrix tau (row=from, col=to):
 [[0.842 0.137 0.021]
 [0.099 0.64  0.261]
 [0.022 0.098 0.88 ]]
emission delta (latent state -> observed level):
 [[0.983 0.012 0.005]
 [0.25  0.703 0.048]
 [0.028 0.092 0.88 ]]
No description has been provided for this image 
Persistence + escalation: non-users stay non-users 84% of the time but 14% drift to occasional;
occasional users are the volatile middle (26% escalate to frequent); frequent use is sticky (88% stay).
The cohort steadily shifts out of non-use across the five years -- the hallmark adolescent escalation curve.

The transition matrix is the payoff that a snapshot LCA cannot give. Non-use is fairly stable but leaks (about one non-user in six begins using each year); occasional use is the unstable middle state, feeding escalation to frequent use; and frequent use is the stickiest state — once there, most stay. Propagating the wave-1 mix through $\tau$ shows the cohort steadily moving out of non-use over the five years. The emission matrix confirms the three latent states are cleanly identified by the observed levels, with modest misclassification (some genuine occasional users report "none").

3. The same model in PyMC — marginalising the state sequence with the forward algorithm¶

NUTS cannot sample the discrete states $S_{i,t}$, so we marginalise them with the forward algorithm: the observed-data likelihood is a product of matrix–vector recursions over the $T=5$ waves, which we simply unroll (no scan needed for a short chain — a bonus on Windows). Dirichlet priors on $\pi$, each row of $\tau$, and each emission $\delta$; the forward log-likelihood goes in with pm.Potential. States are ordered post-hoc by mean emitted level so the transition matrix lines up with the from-scratch fit.

In [4]:
import pymc as pm, pytensor.tensor as pt
Xoh = np.eye(3)[seqs].astype(float)         # (237,5,3) one-hot of the observed level
with pm.Model() as mod:
    pi_   = pm.Dirichlet("pi", np.ones(3))
    tau_  = pm.Dirichlet("tau", np.ones(3), shape=(3,3))
    delta_= pm.Dirichlet("delta", np.ones(3), shape=(3,3))       # C x L (single item)
    B = pt.tensordot(Xoh, delta_, axes=[[2],[1]])                # (N,5,C) emission probs
    a = pi_[None,:]*B[:,0,:]; sc=a.sum(1); a=a/sc[:,None]; ll=pt.log(sc)
    for t in range(1,5):
        a = pt.dot(a, tau_)*B[:,t,:]; sc=a.sum(1); a=a/sc[:,None]; ll=ll+pt.log(sc)
    pm.Potential("like", ll.sum())
    idata = pm.sample(1000, tune=1500, chains=4, target_accept=0.9, random_seed=5, progressbar=False)

# order each posterior draw's states by mean emitted level, then average tau
dd = idata.posterior["delta"].stack(s=("chain","draw")).transpose("s",...).values   # (S,3,3)
tt = idata.posterior["tau"].stack(s=("chain","draw")).transpose("s",...).values      # (S,3,3)
lev=np.arange(3); tau_pm=np.zeros((3,3))
for s in range(dd.shape[0]):
    o=np.argsort((dd[s]*lev).sum(1)); tau_pm += tt[s][o][:,o]
tau_pm/=dd.shape[0]
print("PyMC transition matrix (ordered):\n", tau_pm)
print("from-scratch Gibbs tau:\n", tau)
print(f"\nmax abs difference between the two samplers: {np.abs(tau_pm-tau).max():.3f}")

fig,ax=plt.subplots(1,2,figsize=(11,4.2))
for a_,(Mx,ttl) in zip(ax,[(tau,"from-scratch Gibbs"),(tau_pm,"PyMC (marginalised)")]):
    a_.imshow(Mx,cmap="Oranges",vmin=0,vmax=1)
    for r in range(3):
        for c in range(3): a_.text(c,r,f"{Mx[r,c]:.2f}",ha="center",va="center")
    a_.set_xticks(range(3)); a_.set_yticks(range(3)); a_.set_xticklabels(lab); a_.set_yticklabels(lab)
    a_.set_title("$\\tau$: "+ttl); a_.set_xlabel("to"); a_.set_ylabel("from")
plt.tight_layout(); plt.show()
print("Two samplers, one transition matrix: the from-scratch FFBS Gibbs and the PyMC forward-marginalised HMM agree.")
print("\n(States are relabelled WITHIN each draw before averaging — sorted by mean emitted level — so the")
print("comparison is label-invariant. Any per-component R-hat/ESS warnings above reflect that same")
print("unidentified labelling, not a disagreement about the dynamics.)")

g++ not available, if using conda: `conda install gxx`

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [pi, tau, delta]

Sampling 4 chains for 1_500 tune and 1_000 draw iterations (6_000 + 4_000 draws total) took 9 seconds.

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

PyMC transition matrix (ordered):
 [[0.843 0.135 0.022]
 [0.102 0.64  0.258]
 [0.022 0.099 0.879]]
from-scratch Gibbs tau:
 [[0.842 0.137 0.021]
 [0.099 0.64  0.261]
 [0.022 0.098 0.88 ]]

max abs difference between the two samplers: 0.003
No description has been provided for this image 
Two samplers, one transition matrix: the from-scratch FFBS Gibbs and the PyMC forward-marginalised HMM agree.

(States are relabelled WITHIN each draw before averaging — sorted by mean emitted level — so the
comparison is label-invariant. Any per-component R-hat/ESS warnings above reflect that same
unidentified labelling, not a disagreement about the dynamics.)

4. Summary¶

Latent transition analysis is latent class analysis with a clock: the same measurement model $\delta$ at every wave, plus an initial distribution $\pi$ and a transition matrix $\tau$ that governs how subjects move between classes. On the National Youth Survey marijuana data the three latent use-states reproduce the familiar adolescent pattern — non-use leaks, occasional use is the volatile middle, frequent use is sticky — and the cohort drifts steadily toward use across five years. None of this is visible to a single-wave LCA; it lives entirely in $\tau$.

From scratch the model needed one powerful new ingredient — forward-filter backward-sample to draw each subject's latent trajectory jointly — after which $\pi$, $\tau$ and $\delta$ are all conjugate Dirichlet updates. We validated recovery of a known transition matrix on simulated data (max absolute error ≈0.03) and reproduced the real-data $\tau$ with a PyMC forward-marginalised HMM (the discrete states integrated out, the recursion unrolled over the five waves).

This closes the latent-class arc by connecting it to time series: LTA is a hidden Markov model, and the FFBS sampler here is the same one behind the Markov-switching models in the volatility arc — there the states switch a variance regime, here they switch a response profile. The two arcs meet in the forward–backward algorithm.