Latent Class Regression

Python · R · Download LCA-regression module

Model

Ordinary latent class analysis hands every subject the same prior class weights. But usually there are covariates — age, education, party identification — that ought to shift who lands in which class. Latent class regression, the concomitant-variable model, replaces the fixed weights with a multinomial logit of the covariates. The measurement half is untouched: the classes are still defined by response patterns. What changes is that the mix of classes now bends with who the respondent is, and setting every coefficient to zero recovers ordinary LCA exactly.

\Pr(T_i=c\mid w_i)=\frac{\exp(w_i^{\top}\gamma_c)}{\sum_{k=1}^{C}\exp(w_i^{\top}\gamma_k)}\qquad\text{(membership)}

x_{ij}\mid T_i=c\ \sim\ \text{Categorical}(\delta_{c,j,\cdot})\qquad\text{(measurement)}

$\gamma_c$ — membership log-odds coefficients; all zero recovers ordinary LCA
$\delta_{c,j,\cdot}$ — class $c$ 's response probabilities for item $j$
$w_i$ — the covariate vector; here a 7-point party-identification scale

One new ingredient

From scratch it needs exactly one new ingredient over ordinary LCA. The data-augmentation Gibbs sampler keeps its conjugate blocks — $\delta\mid T$ is Dirichlet, $T\mid\delta,\gamma$ is Categorical — but $\gamma\mid T$ is now a multinomial-logit regression of the imputed labels on the covariates, which is not conjugate and is updated by a small random-walk Metropolis step. That single non-conjugate block is the whole cost of admitting covariates. A PyMC fit with the class marginalised and R's poLCA formula interface corroborate it.

Block	Update	Conjugate?
$\delta\mid T$	Dirichlet	yes
$T\mid\delta,\gamma$	Categorical (the augmentation step)	yes
$\gamma\mid T$	multinomial logit on the imputed labels — random-walk Metropolis	no

What is actually identified

One subtlety the notebooks handle carefully. Multinomial-logit membership coefficients are only identified as contrasts — the raw per-class slopes can all shift by a constant without changing a single fitted probability. So the reported quantity is the contrast: +1.38 log-odds per party point for the pro-Bush class relative to pro-Gore, and $+0.57$ for the ambivalent class. All three engines agree on that number — the from-scratch Gibbs, PyMC's marginalised fit, and the R sampler — and the identified quantities that carry the substance are the membership probabilities $\Pr(T=c\mid w)$ themselves, which is what the crossing curves plot.

Notebooks

Validation comes first, on data with a known answer: the sampler recovers both halves of the model, with membership-curve error under $0.03$ , item-profile error under $0.06$ , and class recovery at an adjusted Rand index of $0.94$ . That matters here more than usual, because it is checking something ordinary LCA cannot even express — not just the measurement model, but how the covariate re-weights the classes.

The application is the classic 2000 American National Election Study data: 1,294 complete respondents rating Gore and Bush on six traits each (12 items, four levels), with 7-point party identification as the covariate. Three latent vote-classes come out exactly as candidate-trait ratings should produce them — one rates Gore well and Bush poorly, one the reverse, and one is lukewarm on both. Then party identification sharply re-weights them: the membership curves cross as you move from strong Democrat to strong Republican, with the ambivalent class peaking in the middle. A fixed-weight LCA reports a single average mix (30% / 28% / 42%) for everyone and misses this entirely.

Python Latent class regression from scratch — data-augmentation Gibbs with a Metropolis step for the membership coefficients, validated on simulated data, applied to the 2000 election study, plus a PyMC marginalised-class fit (Python) Open → R R companion — poLCA’s native formula interface for the concomitant model, a from-scratch base-R Gibbs sampler, and the two engines compared (R) Open →

Downloads

lcareg.py Concomitant-variable LCA — data-augmentation Gibbs with a Metropolis block for the membership coefficients, membership curves, relabelling and the adjusted Rand index (NumPy only) election.csv 2000 American National Election Study — 1,294 respondents, 12 candidate-trait items, party identification and age

LCA-Regression Module — Source Code

"""
lcareg.py -- LATENT CLASS REGRESSION / concomitant-variable LCA (from scratch).

Backs the notebooks in  "Latent Class Regression".

Ordinary LCA has FIXED class weights lambda_c. Latent class regression makes the
weights depend on covariates w_i (age, party, ...) through a MULTINOMIAL LOGIT
(the 'concomitant-variable' model, Dayton & Macready 1988; this is what poLCA fits
when you give it a formula):

    P(T_i = c | w_i) = exp(w_i' gamma_c) / sum_k exp(w_i' gamma_k)   (membership model)
    x_ij | T_i = c  ~ Categorical(delta_{c,j,.})                     (measurement model)

gamma_c are the membership log-odds coefficients (gamma is identified by a ridge prior
gamma_c ~ N(0, sd_gamma^2 I) for ALL classes -- the softmax is invariant to a constant
shift across classes, which the prior removes). delta_{c,j,.} is class c's response-
probability vector for item j over its L levels, with a Dirichlet(1) prior.

The from-scratch sampler is a DATA-AUGMENTATION Gibbs:
  1. delta_{c,j,.} | T   ~ Dirichlet(1 + level counts among class-c members)   (conjugate)
  2. T_i | delta, gamma  ~ Categorical( softmax(w_i'gamma) x prod_j delta ... )  (conjugate)
  3. gamma | T           ~ multinomial-logit regression of the imputed labels on W,
                           updated by a per-class random-walk METROPOLIS step
                           (non-conjugate -- the one Metropolis block in the sampler).

Class labels are permutation-invariant, so posterior draws are RELABELLED post-hoc by
matching each draw's item profiles to a reference draw (Hungarian is overkill for small
C -- we just minimise over the C! permutations).
"""

import numpy as np
from itertools import permutations
from scipy.special import logsumexp


# --------------------------------------------------------------------------- #
#  helpers                                                                      #
# --------------------------------------------------------------------------- #

def _softmax_rows(A):
    return np.exp(A - logsumexp(A, axis=1, keepdims=True))

def _item_loglik(X, delta):
    """log P(x_ij | class) summed over items j, for every class. X:(N,J) int in 0..L-1,
    delta:(C,J,L) -> returns (N,C)."""
    N, J = X.shape; C = delta.shape[0]
    out = np.empty((N, C))
    ar = np.arange(J)
    for c in range(C):
        ld = np.log(delta[c])                       # (J,L)
        out[:, c] = ld[ar[None, :], X].sum(1)       # gather ld[j, X[i,j]] then sum over j
    return out

def _mnl_loglik(gamma, T, W):
    """multinomial-logit log-likelihood of labels T given covariates W and coeffs gamma."""
    eta = W @ gamma.T                               # (N,C)
    return float((eta[np.arange(len(T)), T] - logsumexp(eta, axis=1)).sum())

def membership_probs(gamma, W):
    """predicted class-membership probabilities P(T=c | w) = softmax(W gamma^T)."""
    return _softmax_rows(W @ gamma.T)


# --------------------------------------------------------------------------- #
#  EM-style init: a few covariate-free LCA sweeps to separate the classes       #
# --------------------------------------------------------------------------- #

def _init_labels(X, C, L, rng, iters=25):
    """Quick covariate-free LCA (EM on delta only) to get well-separated starting labels."""
    N, J = X.shape
    delta = rng.dirichlet(np.ones(L), size=(C, J))
    lam = np.full(C, 1.0 / C)
    for _ in range(iters):
        r = np.log(lam)[None, :] + _item_loglik(X, delta)
        r = _softmax_rows(r)                        # (N,C) responsibilities
        lam = r.mean(0) + 1e-8; lam /= lam.sum()
        for c in range(C):
            w = r[:, c]
            for j in range(J):
                cnt = np.array([ (w * (X[:, j] == l)).sum() for l in range(L) ])
                delta[c, j] = (cnt + 1e-3) / (cnt.sum() + L * 1e-3)
    return r.argmax(1)


# --------------------------------------------------------------------------- #
#  Gibbs sampler for latent class regression                                    #
# --------------------------------------------------------------------------- #

def lcareg_gibbs(X, W, C, rng, L=None, draws=4000, burn=2000, sd_gamma=5.0,
                 step=0.08, adapt=True, relabel=True):
    """Data-augmentation Gibbs for concomitant-variable LCA.
    X:(N,J) integer responses coded 0..L-1;  W:(N,p) covariate matrix (include an intercept
    column of ones).  C classes.  Returns posterior draws of the membership coefficients
    gamma:(draws,C,p) and item profiles delta:(draws,C,J,L), relabelled for label switching."""
    X = np.asarray(X); W = np.asarray(W, float)
    N, J = X.shape; p = W.shape[1]
    if L is None:
        L = int(X.max()) + 1
    T = _init_labels(X, C, L, rng)
    gamma = np.zeros((C, p))
    steps = np.full(C, float(step)); acc = np.zeros(C); att = np.zeros(C)
    Bp = 1.0 / sd_gamma ** 2
    GAMMA = np.empty((draws, C, p)); DELTA = np.empty((draws, C, J, L))
    ll_cur = _mnl_loglik(gamma, T, W)
    for it in range(draws + burn):
        # ---- 1. delta | T  (Dirichlet, conjugate) ----
        delta = np.empty((C, J, L))
        for c in range(C):
            Xc = X[T == c]
            for j in range(J):
                cnt = np.bincount(Xc[:, j], minlength=L) if len(Xc) else np.zeros(L)
                delta[c, j] = rng.dirichlet(1.0 + cnt)
        # ---- 2. T | delta, gamma  (Categorical, conjugate) ----
        logp = (W @ gamma.T)                         # membership log-weights (N,C)
        logp = logp - logsumexp(logp, axis=1, keepdims=True)
        logp = logp + _item_loglik(X, delta)
        P = _softmax_rows(logp)
        u = rng.random(N)
        T = (np.cumsum(P, axis=1) < u[:, None]).sum(1)
        T = np.clip(T, 0, C - 1)
        ll_cur = _mnl_loglik(gamma, T, W)
        # ---- 3. gamma | T  (per-class random-walk Metropolis) ----
        for c in range(C):
            prop = gamma.copy()
            prop[c] = gamma[c] + steps[c] * rng.standard_normal(p)
            ll_prop = _mnl_loglik(prop, T, W)
            logr = (ll_prop - ll_cur) - 0.5 * Bp * (prop[c] @ prop[c] - gamma[c] @ gamma[c])
            att[c] += 1
            if np.log(rng.random()) < logr:
                gamma = prop; ll_cur = ll_prop; acc[c] += 1
            if adapt and it < burn and att[c] >= 50:
                rate = acc[c] / att[c]
                steps[c] *= np.exp((rate - 0.3) * 0.5)
                steps[c] = min(max(steps[c], 1e-3), 2.0)
                acc[c] = 0; att[c] = 0
        if it >= burn:
            GAMMA[it - burn] = gamma; DELTA[it - burn] = delta
    if relabel:
        GAMMA, DELTA = _relabel(GAMMA, DELTA)
    return dict(gamma=GAMMA, delta=DELTA, W=W, L=L, accept=float((acc.sum() + 1) / (att.sum() + 1)))


def _relabel(GAMMA, DELTA):
    """Fix label switching: permute each draw's classes to best match a reference draw's
    item profiles (minimise over the C! permutations -- fine for small C)."""
    C = DELTA.shape[1]
    ref = DELTA[0]
    perms = [list(p) for p in permutations(range(C))]
    for d in range(len(DELTA)):
        best = None
        for pm in perms:
            dist = ((DELTA[d][pm] - ref) ** 2).sum()
            if best is None or dist < best[0]:
                best = (dist, pm)
        pm = best[1]
        DELTA[d] = DELTA[d][pm]; GAMMA[d] = GAMMA[d][pm]
    return GAMMA, DELTA


# --------------------------------------------------------------------------- #
#  simulation with a known membership regression                                #
# --------------------------------------------------------------------------- #

def simulate_lcareg(N, gamma, delta, rng, covariate=None):
    """Simulate concomitant-variable LCA data.
    gamma:(C,p) membership coeffs; delta:(C,J,L) item profiles.
    covariate: length-N array of the single covariate (a design column of ones is
    prepended to form W). Returns X:(N,J), W:(N,p), T:(N,)."""
    C, p = gamma.shape; J, L = delta.shape[1], delta.shape[2]
    if covariate is None:
        covariate = rng.integers(1, 8, N).astype(float)      # e.g. 7-point party id
    W = np.column_stack([np.ones(N), covariate]) if p == 2 else \
        np.column_stack([np.ones(N)] + [covariate] * (p - 1))
    P = _softmax_rows(W @ gamma.T)
    T = np.array([rng.choice(C, p=P[i]) for i in range(N)])
    X = np.empty((N, J), int)
    for i in range(N):
        for j in range(J):
            X[i, j] = rng.choice(L, p=delta[T[i], j])
    return X, W, T


def adjusted_rand(a, b):
    """Adjusted Rand index between two label vectors (class-recovery check)."""
    a = np.asarray(a); b = np.asarray(b); n = len(a)
    ca = {v: i for i, v in enumerate(np.unique(a))}; cb = {v: i for i, v in enumerate(np.unique(b))}
    M = np.zeros((len(ca), len(cb)))
    for i in range(n):
        M[ca[a[i]], cb[b[i]]] += 1
    su = M.sum(1); sv = M.sum(0)
    idx = (M * (M - 1) / 2).sum()
    ei = (su * (su - 1) / 2).sum(); ej = (sv * (sv - 1) / 2).sum()
    exp = ei * ej / (n * (n - 1) / 2); mx = 0.5 * (ei + ej)
    return (idx - exp) / (mx - exp) if mx != exp else 1.0

References

Dayton, C. M. & Macready, G. B. (1988). Concomitant-variable latent-class models. Journal of the American Statistical Association 83(401), 173–178. — the original formulation implemented here
Bandeen-Roche, K., Miglioretti, D. L., Zeger, S. L. & Rathouz, P. J. (1997). Latent variable regression for multiple discrete outcomes. Journal of the American Statistical Association 92(440), 1375–1386. — identification and estimation of covariate effects on class membership
Linzer, D. A. & Lewis, J. B. (2011). poLCA: An R package for polytomous variable latent class analysis. Journal of Statistical Software 42(10), 1–29. — the formula interface used as the reference implementation, and the source of the election data
Vermunt, J. K. (2010). Latent class modeling with covariates: two improved three-step approaches. Political Analysis 18(4), 450–469. — why the one-step model fitted here avoids the bias of post-hoc classify-then-regress
The American National Election Studies (2000). ANES 2000 Time Series Study. — the candidate-trait ratings and party-identification scale analysed here