clogit log_lik and prediction bug with group-specific terms

[avehtari]

The following bug report has been made with assistance from Claude AI. I have edited it and I understand everything in this issue.

stan_clogit() leaves a stale intercept column in glmod$X, misaligning coefficients for models with group-specific terms

Summary

stan_clogit() drops the global intercept before fitting, but for models with a
group-specific term it leaves an inconsistent pair of design matrices in the
returned object:

colnames(post$x)         #> "spontaneous" "induced"                 <- what was fitted
colnames(post$glmod$X)   #> "(Intercept)" "spontaneous" "induced"   <- stale

get_x() and model.matrix() read glmod$X for mer models, so they return one
more column than there are coefficients. ll_args.stanreg() and pp_eta() then
select coefficients positionally via beta_sel <- seq_len(ncol(x)), pairing
every column with the wrong coefficient.

Output is silently wrong from log_lik(), loo(), waic(), kfold(),
reloo(), posterior_linpred() and posterior_predict(). Fixed-effects-only
fits are unaffected — there get_x.default() correctly returns object$x.

Reproducible example

The example from the ?stan_clogit help page, unchanged apart from the number of
draws.

library(rstanarm)

dat  <- infert[order(infert$stratum), ]          # order by strata
post <- stan_clogit(case ~ spontaneous + induced + (1 | education),
                    strata = stratum, data = dat, subset = parity <= 2,
                    QR = TRUE, chains = 4, iter = 2000, seed = 1)

loo(post)

Computed from 4000 by 55 log-likelihood matrix.

         Estimate   SE
elpd_loo   -193.6 28.4
p_loo       149.5 27.1
looic       387.1 56.9

Pareto k diagnostic values:
                         Count Pct.
(-Inf, 0.7]   (good)     14    25.5%
   (0.7, 1]   (bad)       0     0.0%
   (1, Inf)   (very bad) 41    74.5%

Two independent signs this cannot be right:

• p_loo = 149.5 for a conditional likelihood with two identified parameters.
• Each stratum has 3 rows and 1 case, so every term is a 3-way softmax and a
  no-signal model scores 55 * log(1/3) = -60.4. The reported elpd_loo = -193.6
  is three times worse than random guessing.

Root cause

In R/stan_clogit.R the intercept is dropped from the local X passed to
stan_glm.fit(), but glmod$X — assigned earlier at X <- glmod$X — is never
updated, and the object stores both (x = X, glmod = glmod).
get_x.lmerMod() and model.matrix.stanreg() return object$glmod$X. Since
as.matrix(post) has no (Intercept) column, seq_len(ncol(x)) = 1:3 produces:

design matrix column	coefficient used
(Intercept)	spontaneous
spontaneous	induced
induced	b[(Intercept) education:0-5yrs]
0-5yrs	b[(Intercept) education:0-5yrs]
6-11yrs	b[(Intercept) education:6-11yrs]
12+ yrs	b[(Intercept) education:12+_yrs]

Within a stratum the (Intercept) and education columns cancel, so the effective
linear predictor collapses to spontaneous * beta_induced + induced * b[edu:0-5yrs].
That last term injects large draw-to-draw variance into the log-likelihood (range
widens from [-6.25, 0] to [-24.20, 0]), which is what produces the heavy-tailed
importance ratios and inflated Pareto-k values.

Note the polr branch of ll_args.stanreg() handles the same situation (no
intercept parameter) correctly, via .validate_polr_x() plus name-based selection
stanmat[, colnames(x), drop = FALSE]. The equivalent guard is missing for clogit.

Verification: the shifted mapping reproduces log_lik(post) exactly

X   <- get_x(post); Z <- as.matrix(get_z(post)); dr <- as.matrix(post)
yv  <- as.vector(get_y(post))
stv <- droplevels(factor(model.frame(post)[, "(weights)"]))   # strata live here

lse     <- function(z) { m <- max(z); m + log(sum(exp(z - m))) }
cond_ll <- function(eta) sapply(levels(stv), function(s) {
  i <- which(stv == s)
  eta[, i[yv[i] == 1]] - apply(eta[, i, drop = FALSE], 1, lse)
})

beta_wrong <- cbind(dr[, seq_len(ncol(X))], dr[, grep("^b\\[", colnames(dr))])
max(abs(cond_ll(tcrossprod(beta_wrong, cbind(X, Z))) - log_lik(post)))
#> 0

Second affected site: newdata

pp_eta() has the same positional selection, so posterior_linpred() is wrong by
up to 24.81 on the linear-predictor scale. Independently, .pp_data_mer_x()
(R/pp_data.R) rebuilds the design matrix from a bars-stripped formula, which
carries an implicit intercept regardless of glmod$X:

colnames(rstanarm:::.pp_data_mer_x(post, newdata = nd))
#> "(Intercept)" "spontaneous" "induced"

so the newdata path needs its own guard, mirroring the existing polr line in
.pp_data().

Fix

Two small changes: keep glmod$X in sync with the fitted X in stan_clogit()
(restoring the contrasts attribute, which column subsetting drops and
.pp_data_mer_x() later reads), and drop the implicit intercept in
.pp_data_mer_x() for clogit. Fixing glmod$X at the source repairs get_x(),
model.matrix(), ll_args.stanreg() and pp_eta() together, so log_lik.R and
posterior_predict.R need no changes.

Worth also switching the two beta_sel <- seq_len(ncol(x)) sites to name-based
selection, so a future mismatch fails loudly rather than silently.

Results after the fix

	before	after
elpd_loo	-193.6	-40.4
p_loo	149.5	2.1
Pareto k > 0.7	41/55 (74.5 %)	none
max abs(log_lik) change	—	24.11
max abs(posterior_linpred) change	—	24.81

p_loo = 2.1 matches the two identified coefficients, and loo() reports "All
Pareto k estimates are good". The matrix is 55 rather than 60 columns in both
cases because loo() already drops the 5 strata where spontaneous and induced
are constant across all three rows — there the conditional likelihood is exactly
1/3 for every draw. Adding their exact log(1/3) contributions back gives
elpd_loo = -45.9 over all 60 strata, comfortably better than the -65.9 null.

Aside: the help-page example

infert is matched on age, parity and education, so education is constant
within all 60 strata and the (1 | education) intercepts cancel identically from
the conditional likelihood — dropping b from the correct linear predictor
changes the log-likelihood by 2.7e-15, and those terms just reproduce their
prior. The example puts a hierarchical term on a matching variable, which carries
no information, and that unidentified term is what triggers the bug. Worth
changing independently of the fix.

Session info

R version 4.6.0 (2026-04-24), x86_64-pc-linux-gnu
rstanarm  GitHub master @ 97fdba4 (DESCRIPTION 2.32.2)
rstan     2.36.0.9000
loo       2.10.0