clmm2 WIP

CLAUDE.md

@@ -2,6 +2,8 @@
 
 This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
 
+Always use `pkgload::load_all()` never `library()`. This lets you use the new version.
+
 ## Overview
 
 `marginaleffects` is an R package that provides a unified interface for computing and plotting predictions, slopes, marginal means, and comparisons (contrasts, risk ratios, odds, etc.) for over 100 classes of statistical and machine learning models. The package supports both unit-level (conditional) and average (marginal) estimates with comprehensive uncertainty quantification.

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: marginaleffects
 Title: Predictions, Comparisons, Slopes, Marginal Means, and Hypothesis Tests
-Version: 0.30.0.4
+Version: 0.30.0.5
 Authors@R:
     c(person(given = "Vincent",
              family = "Arel-Bundock",

NAMESPACE

@@ -22,6 +22,7 @@ S3method(get_coef,betareg)
 S3method(get_coef,bracl)
 S3method(get_coef,brmsfit)
 S3method(get_coef,brmultinom)
+S3method(get_coef,clmm2)
 S3method(get_coef,comparisons)
 S3method(get_coef,data.frame)
 S3method(get_coef,default)
@@ -172,6 +173,7 @@ S3method(sanitize_model_specific,svyolr)
 S3method(set_coef,afex_aov)
 S3method(set_coef,aft)
 S3method(set_coef,betareg)
+S3method(set_coef,clmm2)
 S3method(set_coef,crch)
 S3method(set_coef,data.frame)
 S3method(set_coef,default)

NEWS.md

@@ -2,6 +2,12 @@
 
 ## Development
 
+New:
+
+* Add support for `ordinal::clmm2()` models (cumulative link mixed models).
+
+Misc:
+
 * Add support `MASS::lda()` models for estimates but not standard errors. Thanks to @friendly for feature request #1598.
 
 Bugs:

R/methods_ordinal.R

@@ -59,3 +59,29 @@ sanitize_model_specific.clm <- function(model, ...) {
     }
     return(model)
 }
+
+
+#' @include set_coef.R
+#' @rdname set_coef
+#' @export
+set_coef.clmm2 <- function(model, coefs, ...) {
+    # clmm2 models store coefficients in multiple places:
+    # - model$beta: fixed effects only (named vector)
+    # - model$Alpha (also Theta, xi): threshold parameters (named vector)
+    # - model$coefficients: all parameters including thresholds, fixed effects, and random SD
+    idx_alpha <- seq_len(length(model$Alpha))
+    idx_beta <- seq_len(length(model$beta)) + length(model$Alpha)
+    model$coefficients[seq_len(length(coefs))] <- coefs
+    model$Alpha <- coefs[idx_alpha]
+    model$beta <- coefs[idx_beta]
+    return(model)
+}
+
+
+#' @include get_coef.R
+#' @rdname get_coef
+#' @export
+get_coef.clmm2 <- function(model, ...) {
+    out <- c(model$Alpha, model$beta)
+    return(out)
+}

data-raw/supported_models.csv

@@ -64,6 +64,7 @@ nlme,lme
 nnet,multinom
 ordbetareg,ordbetareg
 ordinal,clm
+ordinal,clmm2
 partykit,glmtree
 partykit,lmtree
 phylolm,phyloglm

inst/tinytest/test-pkg-ordinal.R

@@ -149,3 +149,77 @@ mfx4 <- slopes(mod, variables = "hp", slope = "dyex")
 expect_equivalent(mfx2$estimate, mfx1$estimate * (mfx1$hp / p$estimate))
 expect_equivalent(mfx3$estimate, mfx1$estimate / p$estimate)
 expect_equivalent(mfx4$estimate, mfx1$estimate * mfx1$hp)
+
+
+# clmm2: Basic tests with binomial mixed model
+requiet("lme4")
+cbpp2 <- rbind(lme4::cbpp[,-(2:3)], lme4::cbpp[,-(2:3)])
+cbpp2 <- within(cbpp2, {
+    incidence <- as.factor(rep(0:1, each = nrow(lme4::cbpp)))
+    freq <- with(lme4::cbpp, c(incidence, size - incidence))
+})
+mod_clmm2 <- clmm2(incidence ~ period, random = herd, weights = freq, data = cbpp2, Hess = 1)
+
+# Test basic functionality
+expect_predictions(mod_clmm2)
+expect_comparisons(mod_clmm2)
+expect_slopes(mod_clmm2)
+
+# Test that predictions return reasonable values
+pred <- predictions(mod_clmm2)
+expect_true(all(pred$estimate >= 0 & pred$estimate <= 1))
+expect_true(all(!is.na(pred$estimate)))
+expect_true(all(!is.na(pred$std.error)))
+
+# Test avg_predictions
+avg_pred <- avg_predictions(mod_clmm2)
+expect_inherits(avg_pred, "predictions")
+
+# Test avg_comparisons
+avg_comp <- avg_comparisons(mod_clmm2)
+expect_inherits(avg_comp, "comparisons")
+
+
+# clmm2: Soup data with symmetric threshold (tests set_coef fix)
+dat_soup <- get_dataset("soup", "ordinal")
+dat_soup <- subset(dat_soup, as.numeric(dat_soup$RESP) <= 24)
+# Ensure variables are factors
+dat_soup$SURENESS <- factor(dat_soup$SURENESS, ordered = TRUE)
+dat_soup$RESP <- factor(dat_soup$RESP)
+dat_soup$RESP <- dat_soup$RESP[drop = TRUE]  # Drop unused levels
+dat_soup$PROD <- factor(dat_soup$PROD)  # Ensure PROD is a factor
+mod_soup <- clmm2(SURENESS ~ PROD, random = RESP, data = dat_soup,
+                  link = "probit", Hess = TRUE, method = "ucminf",
+                  threshold = "symmetric")
+
+# Test that set_coef properly updates Theta from Alpha
+expect_predictions(mod_soup)
+expect_comparisons(mod_soup)
+expect_slopes(mod_soup)
+
+# Test predictions
+pred_soup <- predictions(mod_soup)
+expect_true(all(pred_soup$estimate >= 0 & pred_soup$estimate <= 1))
+expect_true(all(!is.na(pred_soup$estimate)))
+expect_true(all(!is.na(pred_soup$std.error)))
+
+# Test avg_predictions by group
+avg_pred_soup <- avg_predictions(mod_soup, by = "PROD")
+expect_inherits(avg_pred_soup, "predictions")
+expect_true(nrow(avg_pred_soup) == 2)  # Two levels of PROD
+
+# Test avg_slopes
+avg_slopes_soup <- avg_slopes(mod_soup)
+expect_inherits(avg_slopes_soup, "slopes")
+expect_true(nrow(avg_slopes_soup) == 1)  # One term (PROD)
+
+# Test hypotheses
+hyp_soup <- hypotheses(mod_soup)
+expect_inherits(hyp_soup, "hypotheses")
+expect_true(nrow(hyp_soup) == 4)  # 3 threshold params + 1 beta
+
+# Test different quadrature methods (with Hess=TRUE for vcov)
+mod_soup_agq <- update(mod_soup, Hess = TRUE, nAGQ = 3)
+pred_agq <- predictions(mod_soup_agq)
+expect_inherits(pred_agq, "predictions")
+expect_true(all(!is.na(pred_agq$estimate)))

trash.qmd

@@ -0,0 +1,45 @@
+---
+title: "Comparing Inference Methods with clmm2 Models"
+format: gfm
+---
+
+Thanks for the nudge @carly-gray-19, I have added support for this model type. Note that the `inferences()` function can apply the same simulation-based inference as the `clarify` package, so we can compare different ways of computing uncertainty. Example:
+
+```{r}
+pkgload::load_all()
+library(ordinal)
+library(lme4)
+
+# Prepare data
+cbpp2 <- rbind(cbpp[,-(2:3)], cbpp[,-(2:3)])
+cbpp2 <- within(cbpp2, {
+    incidence <- as.factor(rep(0:1, each = nrow(cbpp)))
+    freq <- with(cbpp, c(incidence, size - incidence))
+})
+
+# Fit model
+m1 <- clmm2(
+    incidence ~ period,
+    random = herd,
+    weights = freq,
+    data = cbpp2,
+    Hess = 1
+)
+
+# Delta method (default)
+avg_predictions(m1, by = "period")
+
+# Bootstrap
+inferences(
+    avg_predictions(m1, by = "period"),
+    method = "boot",
+    R = 200
+)
+
+# Simulation (like clarify)
+inferences(
+    avg_predictions(m1, by = "period"),
+    method = "simulation",
+    R = 200
+)
+```