Latent Transition Analysis (R)¶
Forward-filter backward-sample from scratch in base R¶
The R counterpart to lta_python.ipynb. Latent transition analysis is a shared hidden Markov model over latent classes: an initial distribution $\pi$, a $C\times C$ transition matrix $\tau$, and the LCA measurement model $\delta$ held invariant across waves. We port the data-augmentation Gibbs to base R — the key step is forward-filter backward-sample (FFBS), drawing each subject's whole latent trajectory jointly — validate it on a simulated cohort, then fit the National Youth Survey marijuana data (237 adolescents, 5 waves).
A note on packages.
LMestis the standard R package for latent Markov / LTA models, but on this Windows build loading its data and fitting routines segfaults the R session, so we keep the notebook fully self-contained in base R (as we do elsewhere when a compiled R package is unstable here). The from-scratch FFBS Gibbs is the same algorithmLMestruns internally, and it cross-checks the Python engine and the PyMC forward-marginalised HMM exactly.
options(repr.plot.width=12, repr.plot.height=4.5)
GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; RED<-"#c53030"; BLUE<-"#2b6cb0"; GREY<-"#718096"
rdir <- function(a){ g<-rgamma(length(a),a,1); g/sum(g) }
cat_rows <- function(P){ u<-runif(nrow(P)); max.col((t(apply(P,1,cumsum))>=u), "first") }
emis_loglik <- function(X, delta, C){ # X: array(N,T,J) in 1..L ; delta: array(C,J,L)
N<-dim(X)[1]; Ti<-dim(X)[2]; J<-dim(X)[3]; logB<-array(0,c(N,Ti,C))
for(c in 1:C){ acc<-matrix(0,N,Ti)
for(j in 1:J){ ld<-log(delta[c,j,]); acc<-acc+matrix(ld[X[,,j]],N,Ti) }
logB[,,c]<-acc }
logB }
ffbs <- function(logB, pi, tau, C){ # sample state sequences S:(N,T)
N<-dim(logB)[1]; Ti<-dim(logB)[2]; mx<-apply(logB,c(1,2),max); B<-exp(logB-array(mx,c(N,Ti,C)))
alpha<-array(0,c(N,Ti,C))
a<-B[,1,]*matrix(pi,N,C,byrow=TRUE); sc<-rowSums(a); alpha[,1,]<-a/sc
for(t in 2:Ti){ pred<-alpha[,t-1,]%*%tau; a<-pred*B[,t,]; sc<-rowSums(a); alpha[,t,]<-a/sc }
S<-matrix(0L,N,Ti); S[,Ti]<-cat_rows(alpha[,Ti,])
for(t in (Ti-1):1){ w<-alpha[,t,]*t(tau[,S[,t+1]]); w<-w/rowSums(w); S[,t]<-cat_rows(w) }
S }
lta_gibbs <- function(X, C, draws=2000, burn=2000){
N<-dim(X)[1]; Ti<-dim(X)[2]; J<-dim(X)[3]; L<-max(X); lev<-1:L
S<-matrix(sample(1:C,N*Ti,replace=TRUE),N,Ti)
PI<-matrix(0,draws,C); TAU<-array(0,c(draws,C,C)); DE<-array(0,c(draws,C,J,L))
for(it in 1:(draws+burn)){
delta<-array(0,c(C,J,L))
for(c in 1:C) for(j in 1:J){ vals<-X[,,j][S==c]; delta[c,j,]<-rdir(1+tabulate(vals,nbins=L)) }
pi<-rdir(1+tabulate(S[,1],nbins=C))
trans<-matrix(0,C,C); for(t in 2:Ti){ tb<-table(factor(S[,t-1],1:C),factor(S[,t],1:C)); trans<-trans+matrix(tb,C,C) }
tau<-t(sapply(1:C, function(a) rdir(1+trans[a,])))
S<-ffbs(emis_loglik(X,delta,C), pi, tau, C)
score<-sapply(1:C, function(c) sum(apply(array(delta[c,,],c(J,L)),2,mean)*lev)) # mean emitted level
o<-order(score); inv<-order(o); delta<-delta[o,,,drop=FALSE]; pi<-pi[o]; tau<-tau[o,o]; S<-matrix(inv[S],N,Ti)
if(it>burn){ PI[it-burn,]<-pi; TAU[it-burn,,]<-tau; DE[it-burn,,,]<-delta } }
list(pi=colMeans(PI), tau=apply(TAU,c(2,3),mean), delta=apply(DE,c(2,3,4),mean)) }
sim_lta <- function(N,Ti,pi,tau,delta){ C<-dim(delta)[1]; J<-dim(delta)[2]; L<-dim(delta)[3]
S<-matrix(0,N,Ti); S[,1]<-sample(1:C,N,replace=TRUE,prob=pi)
for(t in 2:Ti) for(i in 1:N) S[i,t]<-sample(1:C,1,prob=tau[S[i,t-1],])
X<-array(0,c(N,Ti,J)); for(i in 1:N) for(t in 1:Ti) for(j in 1:J) X[i,t,j]<-sample(1:L,1,prob=delta[S[i,t],j,])
list(X=X,S=S) }
cat("from-scratch base-R FFBS Gibbs for latent transition analysis ready\n")
from-scratch base-R FFBS Gibbs for latent transition analysis ready
1. Recovering a known transition matrix¶
800 subjects, $T=3$ waves, two classes (low-/high-risk) each measured by 4 binary items; 20% of low-risk escalate per wave, high-risk is 90% sticky. FFBS should recover $\pi$, $\tau$ and $\delta$.
set.seed(1)
C<-2; J<-4; L<-2
delta_t<-array(0,c(C,J,L)); for(j in 1:J){ delta_t[1,j,]<-c(.85,.15); delta_t[2,j,]<-c(.20,.80) }
pi_t<-c(.7,.3); tau_t<-matrix(c(.8,.2,.1,.9),2,2,byrow=TRUE)
sim<-sim_lta(800,3,pi_t,tau_t,delta_t); fit<-lta_gibbs(sim$X, 2, draws=1500, burn=1500)
cat("pi true", pi_t, " est", round(fit$pi,3), "\n")
cat("tau true:", paste(round(as.vector(t(tau_t)),2),collapse=" "), "\n")
cat("tau est :", paste(round(as.vector(t(fit$tau)),2),collapse=" "), "\n")
cat("tau max abs error", round(max(abs(fit$tau-tau_t)),3), "\n")
par(mfrow=c(1,2), mar=c(4,4,3,1))
for(k in 1:2){ Mx<-if(k==1) tau_t else fit$tau; ttl<-if(k==1)"true tau" else "estimated tau"
image(1:2,1:2,t(Mx[2:1,]),col=hcl.colors(20,"Blues",rev=TRUE),zlim=c(0,1),axes=FALSE,xlab="to class",ylab="from class",main=ttl)
axis(1,1:2,c("low","high")); axis(2,1:2,c("high","low"))
for(r in 1:2) for(cc in 1:2) text(cc, 3-r, sprintf("%.2f",Mx[r,cc])) }
par(mfrow=c(1,1))
cat("FFBS reconstructs the dynamics, not just the classes: the estimated transition matrix matches the truth.\n")
pi true 0.7 0.3 est 0.698 0.302
tau true: 0.8 0.2 0.1 0.9
tau est : 0.79 0.21 0.12 0.88
tau max abs error 0.017
FFBS reconstructs the dynamics, not just the classes: the estimated transition matrix matches the truth.
2. The marijuana data — adolescent use dynamics¶
237 adolescents, 5 annual waves, marijuana use coded 1 = none, 2 = occasional, 3 = frequent. A single indicator per wave, so each latent state is a de-noised use-status with a misclassification distribution; $\tau$ gives the year-to-year dynamics. We fit $C=3$ states, ordered low→high.
seqs <- as.matrix(read.csv("marijuana.csv")) # already coded 1..3
X <- array(seqs, c(nrow(seqs), ncol(seqs), 1))
cat(nrow(seqs), "adolescents, T=", ncol(seqs), "waves\n")
set.seed(2); fit <- lta_gibbs(X, 3, draws=2500, burn=2500)
lab<-c("none","occasional","frequent")
cat("\ninitial mix pi:", paste(sprintf("%s %.2f",lab,fit$pi),collapse=" "), "\n")
cat("transition matrix tau (row=from, col=to):\n"); m<-fit$tau; dimnames(m)<-list(lab,lab); print(round(m,3))
cat("emission delta (state -> level):\n"); e<-fit$delta[,1,]; dimnames(e)<-list(lab,c("none","occ","freq")); print(round(e,3))
prev<-matrix(0,5,3); prev[1,]<-fit$pi; for(t in 2:5) prev[t,]<-prev[t-1,]%*%fit$tau
par(mfrow=c(1,2), mar=c(4,4,3,1))
image(1:3,1:3,t(m[3:1,]),col=hcl.colors(20,"Oranges",rev=TRUE),zlim=c(0,1),axes=FALSE,xlab="to state",ylab="from state",main="Year-to-year transition tau")
axis(1,1:3,lab,cex.axis=.8); axis(2,1:3,rev(lab),cex.axis=.8)
for(r in 1:3) for(cc in 1:3) text(cc,4-r,sprintf("%.2f",m[r,cc]))
cols<-c(GREEN,ORANGE,RED)
matplot(1:5, prev, type="b", pch=19, lwd=2.2, col=cols, ylim=c(0,1), xlab="wave (annual)", ylab="share of cohort", main="Latent use-state prevalence")
legend("topright", lab, col=cols, lwd=2, pch=19, bty="n")
par(mfrow=c(1,1))
cat(sprintf("\nPersistence + escalation: non-users stay %.0f%% but %.0f%% drift up; occasional is the volatile middle\n", 100*m[1,1],100*m[1,2]))
cat(sprintf("(%.0f%% escalate to frequent); frequent use is the stickiest state (%.0f%% stay).\n", 100*m[2,3],100*m[3,3]))
237 adolescents, T= 5 waves
initial mix pi: none 0.91 occasional 0.07 frequent 0.02
transition matrix tau (row=from, col=to):
none occasional frequent none 0.843 0.135 0.022 occasional 0.109 0.629 0.262 frequent 0.021 0.105 0.874
emission delta (state -> level):
none occ freq none 0.984 0.011 0.005 occasional 0.234 0.720 0.046 frequent 0.032 0.085 0.883
Persistence + escalation: non-users stay 84% but 13% drift up; occasional is the volatile middle
(26% escalate to frequent); frequent use is the stickiest state (87% stay).
3. Summary¶
The from-scratch base-R FFBS Gibbs recovers a known transition matrix on simulated data and, on the National Youth Survey marijuana panel, reproduces the Python and PyMC result: three latent use-states with non-use leaking, occasional use as the volatile middle, and frequent use the stickiest state, and a cohort drifting steadily toward use over five years. The whole story lives in the transition matrix $\tau$ — invisible to a single-wave LCA.
LMest is the standard R package for these latent Markov models (and allows covariates on the transitions, à la latent-transition regression), but it is unstable on this Windows build, so we ran the algorithm directly. This closes the latent-class arc and links it to the time-series arc: latent transition analysis is a hidden Markov model, and FFBS is the same forward–backward sampler used for Markov-switching regime models.