Relaxing Local Independence (R)

Conditional dependence in LCA — from scratch, with a randomLCA fit check

The R counterpart to condlca_python.ipynb. We port both samplers from scratch in base R — the naive (local-independence) Gibbs and the random-effects (probit factor) Gibbs via Albert–Chib augmentation — so the bias and its fix reproduce in a second language. The randomLCA package (the standard R tool for random-effects LCA) then gives an independent BIC verdict on whether a random effect is warranted.

Recall the two models. Local independence assumes tests are independent given the class. The random-effects model adds a subject latent trait $b_i\sim N(0,1)$ with $\Pr(x_{ij}=1\mid c,b_i)=\Phi(a_{cj}+\beta_j b_i)$; the loading $\beta_j=0$ recovers local independence, and the marginal accuracies are $\text{Se}_j=\Phi(a_{1j}/\sqrt{1+\beta_j^2})$, $1-\text{Sp}_j=\Phi(a_{0j}/\sqrt{1+\beta_j^2})$. Self-contained (base R) for the samplers; randomLCA used only for the BIC check — and it must run with cores=1 (its default parallel cluster fails under the notebook kernel).

options(repr.plot.width=12, repr.plot.height=4.3)
.libPaths(Sys.getenv("R_LIBS_USER")); suppressMessages(library(randomLCA))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; GREY<-"#718096"; PURP<-"#6b46c1"

# naive local-independence 2-class Gibbs
naive_gibbs <- function(X, draws=3000, burn=1000){
  N<-nrow(X); J<-ncol(X); Se<-runif(J,.5,.9); Sp<-runif(J,.5,.9); pi<-0.4; SE<-matrix(0,draws,J); SP<-matrix(0,draws,J)
  for(t in 1:(draws+burn)){
    l1<-log(pi)+X%*%log(Se)+(1-X)%*%log(1-Se); l0<-log(1-pi)+X%*%log(1-Sp)+(1-X)%*%log(Sp)
    p1<-1/(1+exp(pmin(pmax(l0-l1,-700),700))); Tt<-as.integer(runif(N)<p1); dis<-Tt==1; hea<-!dis
    pi<-rbeta(1,1+sum(dis),1+sum(hea)); sd<-colSums(X[dis,,drop=FALSE]); nd<-sum(dis); Se<-rbeta(J,1+sd,1+nd-sd)
    sh<-colSums(X[hea,,drop=FALSE]); nh<-sum(hea); Sp<-rbeta(J,1+(nh-sh),1+sh)
    if(mean(Se)+mean(Sp)<1){ s2<-1-Sp; Sp<-1-Se; Se<-s2; pi<-1-pi }
    if(t>burn){ SE[t-burn,]<-Se; SP[t-burn,]<-Sp } }
  list(Se=colMeans(SE), Sp=colMeans(SP)) }

# random-effects (probit factor) LCA via Albert-Chib augmentation -- mirrors the Python engine
rtrunc <- function(mean, pos){ lo<-ifelse(pos,pnorm(-mean),0); hi<-ifelse(pos,1,pnorm(-mean))
  u<-lo+matrix(runif(length(mean)),nrow(mean),ncol(mean))*(hi-lo); mean+qnorm(pmin(pmax(u,1e-12),1-1e-12)) }
re_gibbs <- function(X, draws=2500, burn=1000, a_sd=3, beta_sd=0.4){  # beta_sd regularises loadings (0.4 avoids the degenerate all-random-effect collapse)
  N<-nrow(X); J<-ncol(X); Tt<-as.integer(rowMeans(X)>mean(X)); alpha<-matrix(0,2,J); beta<-runif(J,.2,.6); b<-rnorm(N); pi<-0.5
  Aprec<-1/a_sd^2; Bprec<-1/beta_sd^2; SE<-matrix(0,draws,J); SP<-matrix(0,draws,J); BE<-matrix(0,draws,J); PR<-numeric(draws)
  for(it in 1:(draws+burn)){
    mu<-alpha[Tt+1,]+matrix(beta,N,J,byrow=TRUE)*b; z<-rtrunc(mu, X>0.5)
    r<-z-alpha[Tt+1,]; prec_b<-1+sum(beta^2); b<-as.vector(r%*%beta)/prec_b + rnorm(N)/sqrt(prec_b)
    for(c in 0:1){ m<-Tt==c; nc<-sum(m); resid<-z[m,,drop=FALSE]-outer(b[m],beta); prec<-Aprec+nc; alpha[c+1,]<-colSums(resid)/prec+rnorm(J)/sqrt(prec) }
    rj<-z-alpha[Tt+1,]; prec<-Bprec+sum(b^2); beta<-colSums(b*rj)/prec+rnorm(J)/sqrt(prec)
    d1<-z-(matrix(alpha[2,],N,J,byrow=TRUE)+outer(b,beta)); d0<-z-(matrix(alpha[1,],N,J,byrow=TRUE)+outer(b,beta))
    ll1<-log(pi)-0.5*rowSums(d1^2); ll0<-log(1-pi)-0.5*rowSums(d0^2)
    p1<-1/(1+exp(pmin(pmax(ll0-ll1,-700),700))); Tt<-as.integer(runif(N)<p1); pi<-rbeta(1,1+sum(Tt),1+N-sum(Tt))
    if(sum(beta)<0){ beta<--beta; b<--b }
    s<-sqrt(1+beta^2); Se<-pnorm(alpha[2,]/s); Sp<-1-pnorm(alpha[1,]/s)
    if(mean(Se)+mean(Sp)<1){ alpha<-alpha[2:1,]; Tt<-1-Tt; pi<-1-pi; Se<-pnorm(alpha[2,]/s); Sp<-1-pnorm(alpha[1,]/s) }
    if(it>burn){ i<-it-burn; SE[i,]<-Se; SP[i,]<-Sp; BE[i,]<-beta; PR[i]<-pi } }
  list(Se=colMeans(SE), Sp=colMeans(SP), beta=colMeans(BE), prev=mean(PR)) }

# marginal (observed-data) log-likelihood by 15-pt Gauss-Hermite quadrature, integrating out b
GHX <- c(-6.363948,-5.190094,-4.196208,-3.289082,-2.432437,-1.606710,-0.799129,0.000000,0.799129,1.606710,2.432437,3.289082,4.196208,5.190094,6.363948)
GHW <- c(0.00000000,0.00000060,0.00005642,0.00156736,0.01736577,0.08941780,0.23246229,0.31825952,0.23246229,0.08941780,0.01736577,0.00156736,0.00005642,0.00000060,0.00000000)
marg_ll <- function(X, prev, Se, Sp, beta){ J<-ncol(X); s<-sqrt(1+beta^2)
  a1<-qnorm(pmin(pmax(Se,1e-6),1-1e-6))*s; a0<-qnorm(pmin(pmax(1-Sp,1e-6),1-1e-6))*s; Li<-numeric(nrow(X))
  for(q in seq_along(GHX)){ p1<-pnorm(a1+beta*GHX[q]); p0<-pnorm(a0+beta*GHX[q])
    P1<-matrix(p1,nrow(X),J,byrow=TRUE); P0<-matrix(p0,nrow(X),J,byrow=TRUE)
    c1<-exp(rowSums(log(ifelse(X==1,P1,1-P1)))); c0<-exp(rowSums(log(ifelse(X==1,P0,1-P0))))
    Li<-Li+GHW[q]*(prev*c1+(1-prev)*c0) }
  sum(log(Li+1e-300)) }

# MULTI-RESTART: the factor-mixture posterior is multimodal; keep the run with the best marginal fit
re_ms <- function(X, restarts=5, draws=2000, burn=1000, beta_sd=0.4){ best<-NULL; bll<--Inf
  for(r in 1:restarts){ res<-re_gibbs(X, draws=draws, burn=burn, beta_sd=beta_sd)
    ll<-marg_ll(X, res$prev, res$Se, res$Sp, res$beta); if(ll>bll){ bll<-ll; best<-res } }
  best$marg_ll<-bll; best }

sim_dep <- function(N, prev, Se, Sp, beta){ J<-length(Se); s<-sqrt(1+beta^2); a1<-qnorm(Se)*s; a0<-qnorm(1-Sp)*s
  Tt<-as.integer(runif(N)<prev); b<-rnorm(N); A<-ifelse(matrix(Tt,N,J)==1, matrix(a1,N,J,byrow=TRUE), matrix(a0,N,J,byrow=TRUE))
  matrix(as.integer(runif(N*J)<pnorm(A+matrix(beta,N,J,byrow=TRUE)*b)),N,J) }
cat("from-scratch R naive + random-effects (Albert-Chib) samplers ready\n")

Warning message:
"package 'randomLCA' was built under R version 4.6.1"

from-scratch R naive + random-effects (Albert-Chib) samplers ready

1. The bias and the fix, on simulated data

Five tests, two classes; tests 1–3 share a random effect (correlated within class), 4–5 independent. The naive model inflates the correlated tests; the random-effects model recovers the true accuracies, and its loadings flag which tests are dependent.

set.seed(1)
Se_t<-c(.90,.85,.80,.75,.70); Sp_t<-c(.95,.90,.88,.85,.90); beta_t<-c(1.2,1.0,0.9,0,0)
X<-sim_dep(4000, 0.35, Se_t, Sp_t, beta_t)
nv<-naive_gibbs(X, draws=2500, burn=800); re<-re_ms(X, restarts=5, draws=2000, burn=800)
tab<-data.frame(true_Se=Se_t, naive_Se=round(nv$Se,2), RE_Se=round(re$Se,2), true_Sp=Sp_t, naive_Sp=round(nv$Sp,2), RE_Sp=round(re$Sp,2))
rownames(tab)<-paste0("test",1:5); print(tab)
par(mfrow=c(1,2), mar=c(4,4,3,1)); xp<-1:5
plot(xp, Se_t, type="b", pch=19, lwd=1, ylim=c(.6,1), xlab="test", ylab="Se", main="Sensitivity: naive inflates correlated tests")
lines(xp, nv$Se, type="b", pch=15, lty=2, col=RED); lines(xp, re$Se, type="b", pch=17, lty=3, col=BLUE)
legend("bottomleft", c("true","naive (local-indep)","random-effects"), col=c("black",RED,BLUE), lwd=2, pch=c(19,15,17), bty="n", cex=.8)
barplot(re$beta, names.arg=paste0("t",1:5), col=ifelse(re$beta>0.5,PURP,GREY), main="Estimated loadings (flag dependent tests)", ylab="loading beta")
par(mfrow=c(1,1))
cat(sprintf("Naive over-states the correlated tests (Se1 %.2f vs true %.2f); random-effects recovers them (%.2f), loadings large for tests 1-3.\n", nv$Se[1], Se_t[1], re$Se[1]))

      true_Se naive_Se RE_Se true_Sp naive_Sp RE_Sp
test1    0.90     0.93  0.90    0.95     0.97  0.96
test2    0.85     0.90  0.87    0.90     0.92  0.90
test3    0.80     0.84  0.80    0.88     0.89  0.87
test4    0.75     0.73  0.75    0.85     0.86  0.87
test5    0.70     0.68  0.70    0.90     0.89  0.90

Naive over-states the correlated tests (Se1 0.93 vs true 0.90); random-effects recovers them (0.90), loadings large for tests 1-3.

2. Carcinoma — do the pathologists violate local independence?

The naive and random-effects samplers on the carcinoma ratings, plus an independent BIC test from randomLCA: if adding a random effect lowers the BIC, the readers are conditionally dependent.

d<-read.csv("carcinoma.csv"); Xc<-as.matrix(d)-1; raters<-names(d)
set.seed(2); nvC<-naive_gibbs(Xc, draws=4000, burn=1500); reC<-re_ms(Xc, restarts=5, draws=3000, burn=1200)
res<-data.frame(Se_naive=round(nvC$Se,2), Se_RE=round(reC$Se,2), Sp_naive=round(nvC$Sp,2), Sp_RE=round(reC$Sp,2), loading=round(reC$beta,2), row.names=raters)
print(res); cat(sprintf("\nmean |loading| = %.2f (large => strong conditional dependence)\n", mean(abs(reC$beta))))
# independent BIC check with randomLCA (single shared loading for stability; cores=1)
bic<-tryCatch({ f0<-randomLCA(as.data.frame(Xc), nclass=2, random=FALSE, cores=1)
  f1<-randomLCA(as.data.frame(Xc), nclass=2, random=TRUE, quadpoints=31, cores=1)
  c(BIC(f0), BIC(f1)) }, error=function(e) c(NA,NA))
if(!any(is.na(bic))) cat(sprintf("randomLCA BIC: local-independence %.1f vs random-effects %.1f  (drop of %.1f => the random effect is warranted)\n", bic[1], bic[2], bic[1]-bic[2]))
par(mfrow=c(1,2), mar=c(4,4,3,1)); y<-7:1
plot(nvC$Se, y+.12, xlim=c(.2,1), pch=15, col=RED, yaxt="n", ylab="", xlab="probability", main="Naive vs dependence-corrected")
points(reC$Se, y+.12, pch=19, col=BLUE); points(nvC$Sp, y-.12, pch=15, col=ORANGE); points(reC$Sp, y-.12, pch=19, col=GREEN)
axis(2, at=y, labels=raters, las=1); legend("bottomleft", c("Se naive","Se RE","Sp naive","Sp RE"), col=c(RED,BLUE,ORANGE,GREEN), pch=c(15,19,15,19), bty="n", cex=.7)
barplot(reC$beta, names.arg=raters, col=PURP, horiz=TRUE, las=1, main="Reader loadings on shared trait", xlab="loading beta")
par(mfrow=c(1,1))
cat("The loadings are large and randomLCA's BIC prefers the random-effects model: the seven pathologists are strongly\n")
cat("conditionally dependent (= Project 2's third class). The corrected sensitivities are lower -- local independence had\n")
cat("over-credited shared agreement to individual accuracy.\n")

  Se_naive Se_RE Sp_naive Sp_RE loading
A     0.95  0.93     0.88  0.90    0.69
B     0.97  0.97     0.68  0.71    0.77
C     0.71  0.68     0.98  1.00    0.69
D     0.51  0.47     0.98  1.00    0.70
E     0.96  0.94     0.81  0.82    0.76
F     0.40  0.36     0.98  0.99    0.67
G     0.98  0.95     0.91  0.92    0.88

mean |loading| = 0.74 (large => strong conditional dependence)

randomLCA BIC: local-independence 706.1 vs random-effects 668.9  (drop of 37.2 => the random effect is warranted)

The loadings are large and randomLCA's BIC prefers the random-effects model: the seven pathologists are strongly

conditionally dependent (= Project 2's third class). The corrected sensitivities are lower -- local independence had

over-credited shared agreement to individual accuracy.

3. Summary

The from-scratch naive and random-effects (Albert–Chib) samplers reproduce the Python conclusion in base R: ignoring conditional dependence inflates the sensitivity and specificity of correlated tests, and a random-effects model — a latent trait inside the latent class — removes the bias while its loadings reveal which tests are dependent. On the carcinoma ratings the loadings are large and the randomLCA package's BIC prefers the random-effects model, so the seven pathologists genuinely violate local independence (the source of Project 2's third class), and their dependence-corrected accuracies are the more honest audit.

randomLCA (random=TRUE, cores=1) is the standard R tool; for a few known correlated pairs, an explicit fixed conditional-covariance term (Dendukuri–Joseph) is the alternative. This closes the conditional-dependence project and the five-part latent class arc, all built on one conjugate data-augmentation core.