fix clogit loglik and predictions

R/pp_data.R

@@ -170,6 +170,9 @@ pp_data <-
     NULL else lapply(mf[isFac], levels)
   mfnew <- model.frame(delete.response(Terms), newdata, xlev = orig_levs)
   x <- model.matrix(RHS, data = mfnew, contrasts.arg = attr(x, "contrasts"))
+  # RHS has an implicit intercept but stan_clogit has no intercept parameter
+  if (is_clogit(object))
+    x <- x[, colnames(x) != "(Intercept)", drop = FALSE]
   return(x)
 }
 

R/stan_clogit.R

@@ -153,9 +153,13 @@ stan_clogit <- function(formula, data, subset, na.action = NULL, contrasts = NUL
   xint <- match("(Intercept)", colnames(X), nomatch = 0L)
   if (xint > 0L) {
     X <- X[, -xint, drop = FALSE]
+    attr(X, "contrasts") <- contrasts
     # I cannot remember why I was calling drop.terms() to get rid of the intercept
     # mt <- drop.terms(mt, dropx = xint)
     attr(mt, "intercept") <- 0L
+    # get_x() and model.matrix() read glmod$X for mer models, so it must not
+    # keep the intercept column that was just dropped from the fitted X
+    if (has_bars) glmod$X <- X
   }
   f <- binomial(link = "logit")
   stanfit <- stan_glm.fit(x = X, y = Y, weights = weights, 

tests/testthat/test_stan_clogit.R

@@ -68,6 +68,61 @@ test_that("loo/waic for stan_clogit works", {
   expect_identical(ll_fun(fit), rstanarm:::.ll_clogit_i)
 })
 
+SW(fit_mer <- stan_clogit(case ~ spontaneous + induced + (1 | education),
+                          strata = stratum,
+                          data = infert[order(infert$stratum), ],
+                          QR = TRUE, init_r = 0.5,
+                          chains = CHAINS, iter = ITER, seed = SEED, refresh = 0))
+
+# linear predictor built by matching coefficients to columns by name
+clogit_eta_ref <- function(object) {
+  mat <- as.matrix(object)
+  x <- get_x(object)
+  eta <- tcrossprod(mat[, colnames(x), drop = FALSE], x)
+  b <- grep("^b\\[", colnames(mat))
+  if (length(b))
+    eta <- eta + tcrossprod(mat[, b, drop = FALSE], as.matrix(get_z(object)))
+  eta
+}
+
+# one conditional log-likelihood term per stratum; every stratum in infert has
+# exactly one case, so the denominator reduces to a log-sum-exp
+clogit_ll_ref <- function(object) {
+  eta <- clogit_eta_ref(object)
+  y <- as.vector(get_y(object))
+  g <- droplevels(as.factor(model.frame(object)[, "(weights)"]))
+  vapply(levels(g), FUN.VALUE = numeric(nrow(eta)), FUN = function(s) {
+    j <- which(g == s)
+    e <- eta[, j, drop = FALSE]
+    mx <- apply(e, 1, max)
+    eta[, j[y[j] == 1]] - (mx + log(rowSums(exp(e - mx))))
+  })
+}
+
+test_that("stan_clogit with group-specific terms has no intercept column", {
+  expect_false("(Intercept)" %in% colnames(get_x(fit_mer)))
+  expect_identical(colnames(get_x(fit_mer)), colnames(fit_mer$x))
+  expect_identical(colnames(model.matrix(fit_mer)), colnames(fit_mer$x))
+})
+
+test_that("log_lik matches coefficients to columns for stan_clogit with group terms", {
+  ll <- log_lik(fit_mer)
+  g <- droplevels(as.factor(model.frame(fit_mer)[, "(weights)"]))
+  expect_equal(ncol(ll), nlevels(g))
+  expect_equivalent(ll, clogit_ll_ref(fit_mer))
+})
+
+test_that("posterior_linpred matches coefficients to columns for stan_clogit with group terms", {
+  expect_equivalent(posterior_linpred(fit_mer), clogit_eta_ref(fit_mer))
+  expect_equivalent(posterior_linpred(fit_mer, newdata = infert[order(infert$stratum), ]),
+                    clogit_eta_ref(fit_mer))
+})
+
+test_that("loo for stan_clogit with group terms gives a sensible p_loo", {
+  SW(loo_mer <- loo(fit_mer))
+  expect_lt(loo_mer$estimates["p_loo", "Estimate"], 10)
+})
+
 context("posterior_predict (stan_clogit)")
 test_that("compatible with stan_clogit", {
   PPD1 <- posterior_predict(fit)