diff --git a/R/pp_data.R b/R/pp_data.R index 4e42d009e..843b48977 100644 --- a/R/pp_data.R +++ b/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) } diff --git a/R/stan_clogit.R b/R/stan_clogit.R index ac2dec7b4..2917925be 100644 --- a/R/stan_clogit.R +++ b/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, diff --git a/tests/testthat/test_stan_clogit.R b/tests/testthat/test_stan_clogit.R index 532d4314f..049f651f9 100644 --- a/tests/testthat/test_stan_clogit.R +++ b/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)