""" 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