estimate_contrasts: custom treatment perhaps needed for response models (ordinal (cumulative); categorical; etc) models

[DominiqueMakowski]

Some models return predicted values for each response. It's the case for categorical, as well as cumulative models.

Cumulative

estimate_means works as expected:

library(brms)
library(ggplot2)

N <- 400
x <- sample(c(0, 1), size = N, replace = TRUE) 

true_beta <- 1.5           # The true effect of x
tau1 <- -1.0               # Threshold 1 (separates category 1 and 2)
tau2 <- 0.5                # Threshold 2 (separates category 2 and 3)
tau3 <- 2.0                # Threshold 3 (separates category 3 and 4)

latent_y <- true_beta * x + rlogis(N, location = 0, scale = 1)

y_num <- ifelse(latent_y <= tau1, 1,
                ifelse(latent_y <= tau2, 2,
                       ifelse(latent_y <= tau3, 3, 4)))

y <- factor(y_num, ordered = TRUE)
df <- data.frame(x = as.factor(x), y = y)

m <- brm(
  formula = y ~ x,
  data = df,
  family = cumulative(),
  iter = 1000, 
  backend = "cmdstanr",
  algorithm = "pathfinder"
)

modelbased::estimate_means(m, by = "x", predict = "response", backend = "marginaleffects")

x | Response | Median |       95% CI |   pd |          ROPE | % in ROPE
-----------------------------------------------------------------------
0 |        1 | 0.25   | [0.20, 0.32] | 100% | [-0.10, 0.10] |        0%
0 |        2 | 0.38   | [0.32, 0.44] | 100% | [-0.10, 0.10] |        0%
0 |        3 | 0.25   | [0.20, 0.30] | 100% | [-0.10, 0.10] |        0%
0 |        4 | 0.11   | [0.08, 0.15] | 100% | [-0.10, 0.10] |    25.13%
1 |        1 | 0.06   | [0.04, 0.08] | 100% | [-0.10, 0.10] |      100%
1 |        2 | 0.18   | [0.14, 0.22] | 100% | [-0.10, 0.10] |        0%
1 |        3 | 0.35   | [0.30, 0.40] | 100% | [-0.10, 0.10] |        0%
1 |        4 | 0.41   | [0.35, 0.48] | 100% | [-0.10, 0.10] |        0%

But estimate_contrasts combines the response within the level.

modelbased::estimate_contrasts(m, contrast = "x", predict = "response", backend = "marginaleffects")

Level1 | Level2 |    Median |         95% CI |     pd |          ROPE | % in ROPE
---------------------------------------------------------------------------------
1 1    | 1 0    |     -0.20 | [-0.25, -0.15] |   100% | [-0.10, 0.10] |        0%
2 0    | 1 0    |      0.13 | [ 0.02,  0.22] | 98.90% | [-0.10, 0.10] |    27.37%
2 1    | 1 0    |     -0.07 | [-0.15,  0.00] | 97.10% | [-0.10, 0.10] |    76.63%
3 0    | 1 0    | -8.91e-03 | [-0.10,  0.08] | 58.00% | [-0.10, 0.10] |    99.68%
3 1    | 1 0    |      0.09 | [ 0.01,  0.16] | 98.30% | [-0.10, 0.10] |    60.63%
4 0    | 1 0    |     -0.15 | [-0.23, -0.06] |   100% | [-0.10, 0.10] |    11.99%

Categorical

# 1. Define sample size and predictor
N <- 400
x <- sample(c(0, 1), size = N, replace = TRUE) 

# 2. Define true parameters for Category 2 and 3 (relative to Category 1)
# Category 1 is the reference, so its parameters are fixed at 0.
alpha2 <- -0.5   # Intercept for Category 2
beta2  <- 1.0    # Effect of x on Category 2

alpha3 <- 0.5    # Intercept for Category 3
beta3  <- -1.5   # Effect of x on Category 3

# 3. Calculate the linear predictor (log-odds) for each category
eta1 <- rep(0, N)               # Reference category is always 0
eta2 <- alpha2 + beta2 * x
eta3 <- alpha3 + beta3 * x

# 4. Convert log-odds to probabilities using the softmax function
denominator <- exp(eta1) + exp(eta2) + exp(eta3)
p1 <- exp(eta1) / denominator
p2 <- exp(eta2) / denominator
p3 <- exp(eta3) / denominator

# 5. Sample the actual categories based on these probabilities
y_num <- sapply(1:N, function(i) {
  sample(1:3, size = 1, prob = c(p1[i], p2[i], p3[i]))
})

# Convert to an UNORDERED factor (required for categorical models)
y <- factor(y_num, ordered = FALSE)

# Combine into a dataframe
df <- data.frame(x = as.factor(x), y = y)


m <- brm(
  formula = y ~ x,
  data = df,
  family = categorical(),
  iter = 1000, 
  backend = "cmdstanr",
  algorithm = "pathfinder"
)

modelbased::estimate_means(m, by = "x", predict = "response", backend = "marginaleffects")
modelbased::estimate_contrasts(m, contrast = "x", predict = "response", backend = "marginaleffects")

Question

Should we either present that better, or only run the contrasts within every response category? (by default, and allow "Response" to be added to the contrast if someone one also the pairwise contrasts between response levels?)

[DominiqueMakowski]

Change 1: add a Response1 and Response2 columns
Change 2: whether we additionally want to filter by default to keep only Response1 == Response2