Complex comparisons #1329

mattansb · 2025-01-14T08:39:08Z

mattansb
Jan 14, 2025

Currently, comparisons() supports comparing two distinct values of focal factor predictors. However, it is often of interest to compare
groups of levels or to compute otherwise arbitrary contrasts that are based on multiple levels of a factor.

The solution is to do this via the predictions() function with clever uses of hypothesis= and newdata= arguments.

Here are a few examples.

library(marginaleffects)

library(data.table)
library(tidyverse)
library(emmeans)

data("mtcars")

mtcars$carb <- factor(mtcars$carb)

mod <- glm(am ~ carb + mpg, 
           family = binomial("logit"), 
           data = mtcars)

plot_predictions(mod, condition = "carb") + 
  geom_line(aes(carb, estimate), group = 1)

Unit-level (conditional) estimates

We need two things - the first is a data grid that has a all the (relevant) levels of the focal predictor for each row. This is the
default behavior when supplying the variables= argument, but let’s do this manually with datagrid().

grid <- datagrid(model = mod, carb = levels(mtcars$carb), grid_type = "counterfactual")

subset(grid, rowidcf %in% 1:2) # each rowid has all the levels of "carb"
##     rowidcf am mpg carb
## 1         1  1  21    1
## 2         2  1  21    1
## 33        1  1  21    2
## 34        2  1  21    2
## 65        1  1  21    3
## 66        2  1  21    3
## 97        1  1  21    4
## 98        2  1  21    4
## 129       1  1  21    6
## 130       2  1  21    6
## 161       1  1  21    8
## 162       2  1  21    8

Next we need a function that can, for each rowidcf level compute the contrasts of interest. In this case, we are interested in the linear and quadratic trends:

w <- contr.poly(6)[,1:2] # weights
colnames(w) <- c("linear", "quadratic")
rownames(w) <- levels(mtcars$carb)
w
##       linear  quadratic
## 1 -0.5976143  0.5455447
## 2 -0.3585686 -0.1091089
## 3 -0.1195229 -0.4364358
## 4  0.1195229 -0.4364358
## 6  0.3585686 -0.1091089
## 8  0.5976143  0.5455447

`{data.table}`

custom_contrast <- function(x) {
  # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
  setorder(x, carb)[, .(term = colnames(w), 
                        estimate = as.vector(estimate %*% w)), 
                    keyby = "rowidcf"]
}

`{tidyverse}`

custom_contrast <- function(x) {
  x |> 
    # we need to sort because {marginaleffects} cannot be trusted to keep the
    # factor orders
    arrange(carb) |> 
    group_by(rowidcf) |> 
    reframe(
      term = colnames(w),
      estimate = as.vector(estimate %*% w)
    )
}

Bring it all together:

predictions(mod, variables = "carb",
            type = "response", 
            hypothesis = custom_contrast)
## 
##       Term Estimate Std. Error     z Pr(>|z|)    S   2.5 % 97.5 %
##  linear       0.865     0.1486  5.82   <0.001 27.3  0.5736  1.156
##  quadratic    0.185     0.1432  1.29    0.197  2.3 -0.0959  0.466
##  linear       0.865     0.1486  5.82   <0.001 27.3  0.5736  1.156
##  quadratic    0.185     0.1432  1.29    0.197  2.3 -0.0959  0.466
##  linear       0.662     0.1908  3.47   <0.001 10.9  0.2876  1.035
## --- 54 rows omitted. See ?avg_predictions and ?print.marginaleffects --- 
##  quadratic    0.198     0.1308  1.51    0.131  2.9 -0.0586  0.454
##  linear       0.963     0.0120 80.03   <0.001  Inf  0.9390  0.986
##  quadratic    0.395     0.0456  8.65   <0.001 57.5  0.3056  0.485
##  linear       0.829     0.1605  5.16   <0.001 22.0  0.5144  1.144
##  quadratic    0.191     0.1477  1.30    0.195  2.4 -0.0981  0.481
## Type:  response 
## Columns: rowidcf, term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

Average (marginal) estimates.

We can’t use avg_predictions() or predictions(by=) here, because the hypothesis= argument only get applied after the averaging.

So we need to wrap our custom_contrast() function in another function:

`{data.table}`

avg_custom_contrast <- function(x) {
  # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
  setorder(x, carb)[, .(term = colnames(w), 
                        estimate = as.vector(estimate %*% w)), 
                    keyby = "rowidcf"
  ][, .(estimate = mean(estimate)), 
    keyby = "term"]
}

`{tidyverse}`

avg_custom_contrast <- function(x) {
   x |> 
    # we need to sort because {marginaleffects} cannot be trusted to keep the
    # factor orders
    arrange(carb) |> 
    group_by(rowidcf) |> 
    reframe(
      term = colnames(w),
      estimate = as.vector(estimate %*% w)
    ) |> 
    group_by(term) |> 
    summarise(
      estimate = mean(estimate)
    )
}

predictions(mod, variables = "carb",
            type = "response", 
            hypothesis = avg_custom_contrast)
## 
##       Term Estimate Std. Error     z Pr(>|z|)     S 2.5 % 97.5 %
##  linear       0.758     0.0594 12.76   <0.001 121.5 0.642  0.874
##  quadratic    0.326     0.0728  4.48   <0.001  17.0 0.183  0.468
## 
## Type:  response 
## Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

Marginal effects at the mean

See also this amazingness.

This should closely match {emmeans}.

`{data.table}`

custom_contrast <- function(x) {
  # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
  setorder(x, carb)[, .(term = colnames(w), 
                        estimate = as.vector(estimate %*% w))]
}

`{tidyverse}`

custom_contrast <- function(x) {
  x |> 
    # we need to sort because {marginaleffects} cannot be trusted to keep the
    # factor orders
    arrange(carb) |> 
    reframe(
      term = colnames(w),
      estimate = as.vector(estimate %*% w)
    )
}

predictions(mod, variables = "carb",
            type = "response",
            newdata = "mean", hypothesis = custom_contrast)
## 
##       Term Estimate Std. Error    z Pr(>|z|)    S   2.5 % 97.5 %
##  linear       0.924      0.119 7.74   <0.001 46.6  0.6901  1.158
##  quadratic    0.188      0.135 1.40    0.161  2.6 -0.0752  0.452
## 
## Type:  response 
## Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

emmeans(mod, ~ carb, regrid = TRUE) |> 
  contrast(method = "poly", max.deg = 2)
##  contrast  estimate    SE  df z.ratio p.value
##  linear        7.73 0.998 Inf   7.745  <.0001
##  quadratic     1.73 1.230 Inf   1.401  0.1612

Note that {emmeans} scales the estimate differently, but the z- and p-values are the same.

More…

We can further expand these ideas to get conditional contrasts…

`{data.table}`

custom_conditional_contrast <- function(x) {
  # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
  setorder(x, carb, mpg)[, .(term = colnames(w), 
                             estimate = as.vector(estimate %*% w)), 
                         keyby = c("rowidcf", "mpg")
  ][, .(estimate = mean(estimate)), 
    keyby = c("term", "mpg")]
}

`{tidyverse}`

custom_conditional_contrast <- function(x) {
  x |> 
  # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
    arrange(carb, mpg) |> 
    group_by(rowidcf, mpg) |> 
    reframe(
      term = colnames(w), 
      estimate = as.vector(estimate %*% w)
    ) |> 
    group_by(term, mpg) |> 
    summarise(
      estimate = mean(estimate)
    )
}

predictions(mod, variables = list("carb" = unique, "mpg" = "sd"),
            type = "response", 
            hypothesis = custom_conditional_contrast)
## 
##       Term  mpg Estimate Std. Error     z Pr(>|z|)     S   2.5 % 97.5 %
##  linear    17.1    0.970     0.0367 26.43  < 0.001 509.1  0.8981  1.042
##  linear    23.1    0.619     0.1940  3.19  0.00141   9.5  0.2391  0.999
##  quadratic 17.1    0.322     0.0876  3.67  < 0.001  12.0  0.1499  0.493
##  quadratic 23.1    0.263     0.1601  1.64  0.10086   3.3 -0.0511  0.576
## 
## Type:  response 
## Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, mpg

… and interaction contrasts:

`{data.table}`

custom_interaction_contrast <- function(x) {
    # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
  setorder(x, carb, mpg)[, .(term = colnames(w), 
                             estimate = as.vector(estimate %*% w)), 
                         keyby = c("rowidcf", "mpg")
  ][, .(estimate = mean(estimate)), 
    keyby = c("term", "mpg")
  ][, .(estimate = diff(estimate)), 
    keyby = "term"
  ][, term := paste0("(carb:", term, ") (mpg: hi - lo)")]
}

`{tidyverse}`

custom_interaction_contrast <- function(x) {
  x |> 
  # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
    arrange(carb, mpg) |> 
    group_by(rowidcf, mpg) |> 
    reframe(
      term = colnames(w), 
      estimate = as.vector(estimate %*% w)
    ) |> 
    group_by(term, mpg) |> 
    summarise(
      estimate = mean(estimate)
    ) |> 
    group_by(term) |> 
    summarise(
      estimate = diff(estimate)
    ) |> 
    ungroup() |> 
    mutate(
      term = paste0("(carb:", term, ") (mpg: hi - lo)")
    )
}

predictions(mod, variables = list(carb = unique, mpg = "sd"),
            type = "response", 
            hypothesis = custom_interaction_contrast)
## 
##                             Term Estimate Std. Error      z Pr(>|z|)   S  2.5 %
##  (carb:linear) (mpg: hi - lo)     -0.3508      0.181 -1.942   0.0522 4.3 -0.705
##  (carb:quadratic) (mpg: hi - lo)  -0.0589      0.157 -0.376   0.7066 0.5 -0.366
##   97.5 %
##  0.00329
##  0.24793
## 
## Type:  response 
## Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

etc…

This code gives only differences between aggregated means, but the hypothesis=<function> functionality can truly be expanded to arbitrary estimates.

For example, compute the harmonic means, then compare their ratios:

`{data.table}`

avg_custom_contrast <- function(x) {
  setorder(x, carb)[, .(term = "contrast",
                        estimate = mean(estimate[4:6]^-1)^-1 / mean(estimate[1:3]^-1)^-1),
                    keyby = "rowidcf"
  ][, .(estimate = mean(estimate)), keyby = "term"]
}

`{tidyverse}`

avg_custom_contrast <- function(x) {
  x |> 
    arrange(carb) |> 
    group_by(rowidcf) |> 
    summarise(
      term = "contrast",
      estimate = mean(estimate[4:6]^-1)^-1 / mean(estimate[1:3]^-1)^-1
    ) |> 
    group_by(term) |> 
    summarise(
      estimate = mean(estimate)
    )
}

predictions(mod, 
            newdata = datagrid(carb = unique, grid_type = "counterfactual"),
            type = "response", 
            hypothesis = avg_custom_contrast)
## 
##  Estimate Std. Error        z Pr(>|z|)   S     2.5 %   97.5 %
##  10746577   3.97e+10 0.000271        1 0.0 -7.78e+10 7.78e+10
## 
## Term: contrast
## Type:  response 
## Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

Joint tests

Conditional

For the effect of carb

w <- contr.helmert(6)

`{data.table}`

avg_custom_contrast <- function(x) {
  # we need to sort because {marginaleffects} cannot be trusted to keep the
  # factor orders
  setorder(x, carb)[, .(term = 1:5, 
                        estimate = as.vector(estimate %*% w))]
}

`{tidyverse}`

avg_custom_contrast <- function(x) {
   x |> 
    # we need to sort because {marginaleffects} cannot be trusted to keep the
    # factor orders
    arrange(carb) |> 
    reframe(
      term = 1:5,
      estimate = as.vector(estimate %*% w)
    )
}

predictions(mod, variables = list("carb" = unique, mpg = 14),
            type = "response", newdata = "mean",
            hypothesis = avg_custom_contrast) |> 
  hypotheses(joint = TRUE)
## 
## 
## Joint hypothesis test:
##  1 = 0
##  2 = 0
##  3 = 0
##  4 = 0
##  5 = 0
##  
##        F Pr(>|F|) Df 1 Df 2
##  1.9e+07   <0.001    5   25
## 
## Columns: statistic, p.value, df1, df2

predictions(mod, variables = list("carb" = unique, mpg = 26),
            type = "response", newdata = "mean",
            hypothesis = avg_custom_contrast) |> 
  hypotheses(joint = TRUE)
## 
## 
## Joint hypothesis test:
##  1 = 0
##  2 = 0
##  3 = 0
##  4 = 0
##  5 = 0
##  
##     F Pr(>|F|) Df 1 Df 2
##  1559   <0.001    5   25
## 
## Columns: statistic, p.value, df1, df2

joint_tests(mod, by = "mpg", regrid = TRUE, at = list(mpg = c(14, 26)))
## mpg = 14:
##  model term df1 df2      F.ratio    Chisq p.value
##  carb         5 Inf 18356583.000 91782913  <.0001
## 
## mpg = 26:
##  model term df1 df2      F.ratio    Chisq p.value
##  carb         5 Inf     1559.000     7795  <.0001

vincentarelbundock · 2025-01-14T11:37:47Z

vincentarelbundock
Jan 14, 2025
Maintainer

Ooooh, this looks super cool!

Life is crazy right now, but I look forward to checking this out carefully.

0 replies

vincentarelbundock · 2025-01-18T21:33:53Z

vincentarelbundock
Jan 18, 2025
Maintainer

That's interesting. I'll take a close look when I have some time.

But just to fix expectations, the design I currently have in mind is to expand the formula interface to cover most of the usual contrasts, and to steer other users toward custom functions. I think the interface for specify_hypothesis() is extremely complex, and I don't want to add that complexity to the interface.

5 replies

mattansb Jan 18, 2025
Author

I think the current interface for the hypothesis= argument is super flexible, even if we completely ignore the arbitrary function option.

But, for more complex contrasts (as the ones I've outlined here), there is currently only one option - writing a function. I think the hypotheses_spec() supports this - simplifying complex manipulations by breaking them down into smaller bits, and allowing them to work conditionally (as in, group-wise).

If you find that you instead want to expand the formula interface, and send anyone wanting a more complex than that subset of operations to write their own functions, I think the original post still has that covered 100% (and can be polished up to become a "case study"?)

vincentarelbundock Jan 18, 2025
Maintainer

Yes, I do want to make the formula interface more powerful. What other contrast styles do people need?

And you're right that the spec_hypotheses() function is cool and powerful. Maybe it would be worth it to write a vignette about something like that.

If I'm honest, I think that part of my resistance to exposing these features is kind of "ideological." I keep repeating to everyone who will listen that they shouldn't report odds ratios, because they are unintelligible. But then, the website shows how to compute group-wise ratios between cross-cubic contrast and a linear difference, or whatever...

Deep down, I feel that quantities like these are fundamentally uninterpretable. That comment probably just exposes my own limitations, but I simply can't wrap my brain around them. And I feel that things like these kind of contradict the spirit of "What is your estimand?"

That said, I understand that other people find them useful. I'm just giving you personal background on the "intuitive" and "ideological" and probably terrible reasons why I've been hesitant to build infrastructure in that direction.

mattansb Jan 19, 2025
Author

Maybe the polynomial example is drawing too much focus here -
I come from an experimental background, heavy on multi-leveled categorical predictors, where being able to estimate (and test) compound comparisons (comparing more than 2 weighted estimates) is very common and useful.

Some examples:

If we have 3 groups: control, treatment1, and treatment 2, the compound contrast [-1, 0.5, 0.5] gives an estimate of an effect of treatment (any treatment).
The polynomial contrasts are typically used with discrete ordinal predictors - ones that aren't quantitative, but their order has meaning. Outside Bayesian contexts, ordinal predictors are treated as categorical predictors, but we can get estimates for the monotonic effect with a linear contrast (the exact scaling of the weights can be used to give more or less specific meaning).

We can conditional effects like these, and compare them between sub-populations (interaction contrasts). All very useful (and have I mentioned common) operations.

I understand your resistance is due to the opaqueness of compound comparisons, but as a user of such comparisons I see no fundamental difference between them and simple (two level at a time) comparisons (comparisons()).

vincentarelbundock Jan 19, 2025
Maintainer

Right, makes sense. That's totally fair.

I was looking at the hypothesis code yesterday, and it has grown into a horrible monstrosity. I was hoping to publish version 1.0.0 soon, but hypothesis needs a complete re-write before that.

I'll give this a shot when I have time. Then, I'll implement more contrasts via the formula interface.

You [-1, 0.5, 0.5] example is a very good one, and I think there should be an easy way to specify that.

I agree that polynomials should also be easy.

When those points are all sorted, I'll circle back to the more complex cases.

Thanks for your patience!!

mattansb Jan 19, 2025
Author

Cool (: Let me know if you want some feedback on user-facing things. I still think the hypotheses_spec() is something worth adding, but just like as.emmGrid(), I might develop it just for my own use (;

mattansb · 2025-01-22T11:15:05Z

mattansb
Jan 22, 2025
Author

Updating this example - it works great! (still missing #1344)

@vincentarelbundock I can probably cook up a better example that actually makes sense that I would be happy to contribute as a "case study". Let me know what you think (:

This is essentially a show-and-tell of solutions to #319, #479 &
#1232.

Complex contrasts

Currently, comparisons() supports comparing two distinct values of
focal factor predictors. However, it is often of interest to compare
groups of levels or to compute otherwise arbitrary contrasts that
are based on multiple levels of a factor.

The solution is to do this via the predictions() function with clever
uses of hypothesis= and newdata= arguments.

Here are a few examples.

# pak::pak("vincentarelbundock/marginaleffects")
library(marginaleffects)
packageVersion("marginaleffects")
## [1] '0.24.0.15'

library(emmeans)

data("mtcars")

mtcars$carb <- factor(mtcars$carb)

mod <- glm(am ~ carb + mpg,
           family = binomial("logit"),
           data = mtcars)

plot_predictions(mod, condition = "carb") +
  ggplot2::geom_line(ggplot2::aes(carb, estimate), group = 1)

Unit-level (conditional) estimates

We need two things - the first is a data grid that has a all the
(relevant) levels of the focal predictor for each row. This is the
default behavior when supplying the variables= argument, but let’s do
this manually with datagrid().

grid <- datagrid(model = mod, carb = levels(mtcars$carb), grid_type = "counterfactual")

subset(grid, rowidcf %in% 1:2) # each rowid has all the levels of "carb"
##     rowidcf am mpg carb
## 1         1  1  21    1
## 2         2  1  21    1
## 33        1  1  21    2
## 34        2  1  21    2
## 65        1  1  21    3
## 66        2  1  21    3
## 97        1  1  21    4
## 98        2  1  21    4
## 129       1  1  21    6
## 130       2  1  21    6
## 161       1  1  21    8
## 162       2  1  21    8

Next we need a function that can, for each rowidcf level compute the
contrasts of interest. In this case, we are interested in the linear and
quadratic trends:

w <- contr.poly(6)[,1:2] # weights
colnames(w) <- c("linear", "quadratic")
rownames(w) <- levels(mtcars$carb)
w
##       linear  quadratic
## 1 -0.5976143  0.5455447
## 2 -0.3585686 -0.1091089
## 3 -0.1195229 -0.4364358
## 4  0.1195229 -0.4364358
## 6  0.3585686 -0.1091089
## 8  0.5976143  0.5455447

custom_contrast <- function(x) {
  setNames(
    as.vector(x %*% w),
    nm = colnames(w)
  )
}

Bring it all together:

h1 <- predictions(mod, variables = list("carb" = levels),
                  type = "response",
                  hypothesis = ~ I(custom_contrast(x)) | rowidcf)

nrow(mtcars)
## [1] 32
nrow(h1) # two contrasts per obs.
## [1] 64

h1
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)    S   2.5 % 97.5 %
##   linear       0.865     0.1486  5.82   <0.001 27.3  0.5736  1.156
##   quadratic    0.185     0.1432  1.29    0.197  2.3 -0.0959  0.466
##   linear       0.865     0.1486  5.82   <0.001 27.3  0.5736  1.156
##   quadratic    0.185     0.1432  1.29    0.197  2.3 -0.0959  0.466
##   linear       0.662     0.1908  3.47   <0.001 10.9  0.2876  1.035
## --- 54 rows omitted. See ?print.marginaleffects --- 
##   quadratic    0.198     0.1308  1.51    0.131  2.9 -0.0586  0.454
##   linear       0.963     0.0120 80.03   <0.001  Inf  0.9390  0.986
##   quadratic    0.395     0.0456  8.65   <0.001 57.5  0.3056  0.485
##   linear       0.829     0.1605  5.16   <0.001 22.0  0.5144  1.144
##   quadratic    0.191     0.1477  1.30    0.195  2.4 -0.0981  0.481
## Type:  response

Average (marginal) estimates.

We can’t use avg_predictions() or predictions(by=) here, because the
hypothesis= argument only get applied after the averaging.

But we can pass the results fo further processing to hypotheses()

h2 <- hypotheses(h1, hypothesis = ~ I(mean(x)) | hypothesis)

h2
## 
##       Term Estimate Std. Error     z Pr(>|z|)     S 2.5 % 97.5 %
##  linear       0.758     0.0594 12.76   <0.001 121.5 0.642  0.874
##  quadratic    0.326     0.0728  4.48   <0.001  17.0 0.183  0.468

Marginal effects at the mean

More…

conditional contrasts

We can further expand these ideas to get conditional contrasts…

h3 <- predictions(mod, variables = list("carb" = levels, "mpg" = "sd"),
                  type = "response",
                  hypothesis = ~ I(custom_contrast(x)) | rowidcf + mpg)

h3
## 
##   mpg Hypothesis Estimate Std. Error     z Pr(>|z|)     S   2.5 % 97.5 %
##  17.1  linear       0.970     0.0367 26.43  < 0.001 509.1  0.8981  1.042
##  17.1  quadratic    0.322     0.0876  3.67  < 0.001  12.0  0.1499  0.493
##  23.1  linear       0.619     0.1940  3.19  0.00141   9.5  0.2391  0.999
##  23.1  quadratic    0.263     0.1601  1.64  0.10086   3.3 -0.0511  0.576
##  17.1  linear       0.970     0.0367 26.43  < 0.001 509.1  0.8981  1.042
## --- 118 rows omitted. See ?print.marginaleffects --- 
##  23.1  quadratic    0.263     0.1601  1.64  0.10086   3.3 -0.0511  0.576
##  17.1  linear       0.970     0.0367 26.43  < 0.001 509.1  0.8981  1.042
##  17.1  quadratic    0.322     0.0876  3.67  < 0.001  12.0  0.1499  0.493
##  23.1  linear       0.619     0.1940  3.19  0.00141   9.5  0.2391  0.999
##  23.1  quadratic    0.263     0.1601  1.64  0.10086   3.3 -0.0511  0.576
## Type:  response

h4 <- hypotheses(h3, hypothesis = ~ I(mean(x)) | hypothesis + mpg)
h4
## 
##       Term Estimate Std. Error     z Pr(>|z|)     S   2.5 % 97.5 %
##  linear       0.970     0.0367 26.43  < 0.001 509.1  0.8981  1.042
##  linear       0.619     0.1940  3.19  0.00141   9.5  0.2391  0.999
##  quadratic    0.322     0.0876  3.67  < 0.001  12.0  0.1499  0.493
##  quadratic    0.263     0.1601  1.64  0.10086   3.3 -0.0511  0.576

interaction contrasts

… and interaction contrasts:

custom_interaction_contrast <- function(x) {
  setNames(
    diff(x),
    "mpg: hi - lo"
  )
}

hypotheses(h4, hypothesis = ~ I(custom_interaction_contrast(x)) | term)
## 
##  Estimate Std. Error      z Pr(>|z|)   S  2.5 %  97.5 %
##   -0.3508      0.181 -1.942   0.0522 4.3 -0.705 0.00329
##   -0.0589      0.157 -0.376   0.7066 0.5 -0.366 0.24793
## 
## Term: mpg: hi - lo

etc…

Joint tests

Conditional

For the effect of carb

w <- contr.helmert(6)

avg_custom_contrast <- function(x) {
  setNames(
    as.vector(x %*% w),
    nm = 1:5
  )
}

h5 <- predictions(mod, variables = list("carb" = levels, mpg = c(14, 26)),
                  type = "response", newdata = "mean",
                  hypothesis = ~ I(avg_custom_contrast(x)) | mpg)

h5
## 
##  mpg Hypothesis Estimate Std. Error      z Pr(>|z|)    S   2.5 %  97.5 %
##   14          1 -0.00151    0.00567 -0.266    0.790  0.3 -0.0126  0.0096
##   14          2 -0.00635    0.01449 -0.438    0.662  0.6 -0.0348  0.0221
##   14          3  0.16852    0.21821  0.772    0.440  1.2 -0.2592  0.5962
##   14          4  3.93536    0.08672 45.380   <0.001  Inf  3.7654  4.1053
##   14          5  3.93537    0.08642 45.539   <0.001  Inf  3.7660  4.1047
##   26          1 -0.06812    0.21718 -0.314    0.754  0.4 -0.4938  0.3575
##   26          2 -1.65814    0.30747 -5.393   <0.001 23.8 -2.2608 -1.0555
##   26          3  1.31185    0.27035  4.852   <0.001 19.6  0.7820  1.8417
##   26          4  1.35184    0.32103  4.211   <0.001 15.3  0.7226  1.9811
##   26          5  1.35184    0.32103  4.211   <0.001 15.3  0.7226  1.9811
## 
## Type:  response

emc <- emmeans(mod, ~ carb + mpg, regrid = TRUE, at = list(mpg = c(14, 26))) |>
  contrast(method = as.data.frame(w), by = "mpg")
emc
## mpg = 14:
##  contrast estimate      SE  df z.ratio p.value
##  V1       -0.00151 0.00567 Inf  -0.266  0.7901
##  V2       -0.00634 0.01450 Inf  -0.438  0.6616
##  V3        0.16852 0.21800 Inf   0.772  0.4399
##  V4        3.93536 0.08670 Inf  45.381  <.0001
##  V5        3.93537 0.08640 Inf  45.539  <.0001
## 
## mpg = 26:
##  contrast estimate      SE  df z.ratio p.value
##  V1       -0.06813 0.21700 Inf  -0.314  0.7538
##  V2       -1.65814 0.30700 Inf  -5.393  <.0001
##  V3        1.31185 0.27000 Inf   4.852  <.0001
##  V4        1.35184 0.32100 Inf   4.211  <.0001
##  V5        1.35184 0.32100 Inf   4.211  <.0001

hypotheses(h5, joint = 1:5)
## 
## 
## Joint hypothesis test:
##  b1 = 0
##  b2 = 0
##  b3 = 0
##  b4 = 0
##  b5 = 0
##  
##        F Pr(>|F|) Df 1 Df 2
##  1.9e+07   <0.001    5   25
hypotheses(h5, joint = 6:10)
## 
## 
## Joint hypothesis test:
##  b6 = 0
##  b7 = 0
##  b8 = 0
##  b9 = 0
##  b10 = 0
##  
##     F Pr(>|F|) Df 1 Df 2
##  1559   <0.001    5   25

joint_tests(mod, by = "mpg", regrid = TRUE, at = list(mpg = c(14, 26)))
## mpg = 14:
##  model term df1 df2      F.ratio    Chisq p.value
##  carb         5 Inf 18356583.000 91782913  <.0001
## 
## mpg = 26:
##  model term df1 df2      F.ratio    Chisq p.value
##  carb         5 Inf     1559.000     7795  <.0001

0 replies

vincentarelbundock · 2025-01-22T12:47:35Z

vincentarelbundock
Jan 22, 2025
Maintainer

Yep, I think it's a good idea to add a couple examples of the custom contrast functions via formula interface. Note that we already have a bunch of examples of grouping with formula here:

https://marginaleffects.com/bonus/hypothesis.html#formulas-group-wise

The website code is in a private repo because I'm working on the book there, but I paste the relevant section code (folded below) in case you want to play with some options.

Two minor style things:

We shouldn't assign or modify built in datasets like mtcars. Create a new dat instead.
It's good to add examples, but if the website is too complex, people just won't read it. So I would like to keep the examples very minimalist. Maybe 2 or 3 examples of custom functions, but let's strip down other concepts to the maximum. For example, the idea of a counterfactual grid can be combined for good effect, but it is another layer of complexity that we don't need just to illustrate the functionality. There can be blog posts elsewhere to go deep in certain use-cases (like Heiss's posts).

Thanks for engaging! I'm excited about the new features we ended up with.

````r # Formulas (Group-wise)

Since version 0.20.1.5 of marginaleffects, the hypothesis argument also accepts a formula to specify hypotheses to be tested on every row of the output, or within subsets. Consider this set of average predictions:

dat <- within(mtcars, {
    gear = factor(gear)
    cyl = factor(cyl)
})
dat <- sort_by(dat, ~ gear + cyl + am)
mod <- lm(mpg ~ am * wt * factor(cyl), data = dat)

avg_predictions(mod, by = "gear")

We can make "sequential", "reference", or "meandev" comparisons, either as differences (the default) or ratios:

avg_predictions(mod, 
    hypothesis = ~ sequential,
    by = "gear")

avg_predictions(mod, 
    hypothesis = ~ meandev,
    by = "gear")

avg_predictions(mod, 
    hypothesis = ratio ~ sequential,
    by = "gear")

We can also test hypotheses by subgroup. For example, in this case, we want to compare every estimate to the reference category (first estimate) in each subset of am:

avg_predictions(mod, 
    by = c("am", "gear")
)

avg_predictions(mod, 
    hypothesis = ~ sequential | am,
    by = c("am", "gear")
)

The grouping component of the formula is particularly useful when conducting tests on contrasts:

avg_comparisons(mod, 
    variables = "cyl",
    by = "am")

avg_comparisons(mod, 
    variables = "cyl",
    by = "am", 
    hypothesis = ~ sequential | contrast)

avg_comparisons(mod, 
    variables = "cyl",
    by = "am", 
    hypothesis = ~ sequential | am)

It is also a powerful tool in categorical outcome models, to compare estimates for different outcome levels:

library("MASS")

mod <- polr(gear ~ cyl + hp, dat)

avg_comparisons(mod, variables = "cyl")

Is the "8 vs. 4" contrast equal to the "6 vs. 4" contrast, within each outcome level?

avg_comparisons(mod, 
    variables = "cyl",
    hypothesis = ~ sequential | group)

Is the "8 vs. 4" contrast for outcome level 3 equal to the "8 vs. 4" contrast for outcome level 4?

avg_comparisons(mod, 
    variables = "cyl",
    hypothesis = ~ sequential | contrast)

Functions

Hypothesis tests can also be conducted using arbitrary R functions. This allows users to test hypotheses on complex aggregations or transformations. To achieve this, we define a custom function which accepts a fitted model or marginaleffects objects, and returns a data frame with at least two columns: term (or hypothesis) and estimate.

Predictions

When supplying a function to the hypothesis argument, that function must accept an argument x which is a data frame with columns rowid and estimate (and optional columns for other elements of newdata). That function must return a data frame with columns term (or hypothesis) and estimate.

In this example, we test if the mean predicted value is different from 2:

dat <- transform(mtcars, gear = factor(gear), cyl = factor(cyl))

mod <- lm(wt ~ mpg * hp * cyl, data = dat)

hyp <- function(x) {
    data.frame(
        hypothesis = "Avg(Ŷ) = 2",
        estimate = mean(x$estimate) - 2
    )
}
predictions(mod, hypothesis = hyp)

In this ordinal logit model, the predictions() function returns one row per observation and per level of the outcome variable:

#| messages: false
library(MASS)
library(dplyr)

mod <- polr(gear ~ cyl + hp, dat)

avg_predictions(mod)

We can use a function in the hypothesis argument to collapse the rows, displaying the average predicted values in groups 3-4 vs. 5:

fun <- function(x) {
    out <- x |> 
        mutate(term = ifelse(group %in% 3:4, "3 & 4", "5")) |>
        summarize(estimate = mean(estimate), .by = term)
    return(out)
}
avg_predictions(mod, hypothesis = fun)

And we can compare the two categories by doing:

fun <- function(x) {
    out <- x |> 
        mutate(term = ifelse(group %in% 3:4, "3 & 4", "5")) |>
        summarize(estimate = mean(estimate), .by = term) |>
        summarize(estimate = diff(estimate), term = "5 - (3 & 4)")
    return(out)
}
avg_predictions(mod, hypothesis = fun)

</details>

7 replies

mattansb Jan 22, 2025
Author

@vincentarelbundock How's this? To be added at the Formulas section and the end of the lm() example and at the end of the polr() example.

We can also pass simple functions to manipulate the vector of estimates
via the right hand side of the formula. For example we can compare the
estimated effect of wt for cyl = 8 to the average estimated effect
for cyl = c(4, 6). We want both the difference and the ratio. (Note
that we ordered dat so we know the which level of cyl each
estimate corresponds to):

less_vs_8 <- function(x) {
  c("{8} - {4, 6}" = x[3] - mean(x[1:2]),
    "{8} / {4, 6}" = x[3] / mean(x[1:2]))
}

avg_slopes(mod, 
    variables = "wt",
    by = "cyl",
    hypothesis = ~ I(less_vs_8(x)))
## 
##    Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##  {8} - {4, 6}    3.957      4.805 0.823    0.410 1.3 -5.461  13.37
##  {8} / {4, 6}    0.375      0.357 1.052    0.293 1.8 -0.324   1.07
## 
## Type:  response

We can also combine these functions with a grouping component:

avg_slopes(mod, 
    variables = "wt",
    by = c("cyl", "am"),
    hypothesis = ~ I(less_vs_8(x)) | am)
## 
##  am   Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##   0 {8} - {4, 6}    2.597       6.43 0.404    0.686 0.5 -10.01  15.20
##   0 {8} / {4, 6}    0.484       0.64 0.757    0.449 1.2  -0.77   1.74
##   1 {8} - {4, 6}    2.243      10.90 0.206    0.837 0.3 -19.12  23.61
##   1 {8} / {4, 6}    0.471       2.15 0.219    0.826 0.3  -3.74   4.69
## 
## Type:  response

It is also a powerful tool in categorical outcome models […]

Here too we can a custom function. For example, we can conditionally
estimate the average rank why treating out estimate as the relative
“weight” of each level:

mean_rank <- function(x) {
  weighted.mean(c(3, 4, 5), w = x)
}

ranks <- avg_predictions(mod,
           variables = "cyl", 
           hypothesis = ~ I(mean_rank(x)) | cyl)
## 
## Re-fitting to get Hessian
ranks
## 
##  cyl Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
##    4     4.64     0.0885 52.4   <0.001   Inf  4.46   4.81
##    6     4.09     0.1156 35.4   <0.001 908.0  3.86   4.32
##    8     3.12     0.0375 83.2   <0.001   Inf  3.04   3.19
## 
## Type:  probs

These can then be further compared using the hypotheses() function:

hypotheses(ranks, ratio ~ pairwise)
## 
##   Term Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
##  6 / 4    0.882     0.0296 29.8   <0.001 647.9 0.824  0.940
##  8 / 4    0.672     0.0150 44.7   <0.001   Inf 0.643  0.702
##  4 / 6    1.134     0.0380 29.8   <0.001 647.9 1.059  1.208
##  8 / 6    0.762     0.0234 32.6   <0.001 770.5 0.716  0.808
##  4 / 8    1.487     0.0333 44.7   <0.001   Inf 1.422  1.553
##  6 / 8    1.312     0.0403 32.6   <0.001 770.5 1.233  1.391

vincentarelbundock Jan 22, 2025
Maintainer

FYI, I plan to release the new features to CRAN in the next couple weeks.

vincentarelbundock Jan 22, 2025
Maintainer

Yes, I think that's great!

Maybe show a simpler example of unnamed function first, like hypothesis = I(mean(x)[1:2]) | groupid

mattansb Jan 22, 2025
Author

Good idea.

Here is an updated version, to be added at the Formulas section and the end of the lm() example and at the end of the polr() example.

I've called the stuff inside I() "expressions".

We can also pass custom expressions to manipulate the vector of
estimates by using I() in the right hand side of the formula. For
example, if we want the mean slope for cyl = c(6, 8) we can express
this as mean(x[2:3]) (with x representing the vector of estimates):

avg_slopes(mod,
    variables = "wt",
    by = "cyl",
    hypothesis = ~ I(mean(x[2:3])))
## 
##  Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
##     -5.07       4.53 -1.12    0.263 1.9 -13.9    3.8
## 
## Type:  response

(Note that we ordered dat so we know the which level of cyl each
estimate corresponds to):

These expressions can return multiple, named, estiamtes. Often is more
conviniant to wrap these in a function that will then be called in the
I() expression. For example we can compare the estimated effect of
wt for cyl = 8 to the average estimated effect for cyl = c(4, 6).
We want both the difference and the ratio.

less_vs_8 <- function(x) {
  c("{8} - {4, 6}" = x[3] - mean(x[1:2]),
    "{8} / {4, 6}" = x[3] / mean(x[1:2]))
}

avg_slopes(mod,
    variables = "wt",
    by = "cyl",
    hypothesis = ~ I(less_vs_8(x)))
## 
##    Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##  {8} - {4, 6}    3.957      4.805 0.823    0.410 1.3 -5.461  13.37
##  {8} / {4, 6}    0.375      0.357 1.052    0.293 1.8 -0.324   1.07
## 
## Type:  response

We can also combine these expressions with a grouping component:

avg_slopes(mod,
    variables = "wt",
    by = c("cyl", "am"),
    hypothesis = ~ I(less_vs_8(x)) | am)
## 
##  am   Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##   0 {8} - {4, 6}    2.597       6.43 0.404    0.686 0.5 -10.01  15.20
##   0 {8} / {4, 6}    0.484       0.64 0.757    0.449 1.2  -0.77   1.74
##   1 {8} - {4, 6}    2.243      10.90 0.206    0.837 0.3 -19.12  23.61
##   1 {8} / {4, 6}    0.471       2.15 0.219    0.826 0.3  -3.74   4.69
## 
## Type:  response

It is also a powerful tool in categorical outcome models […]

Here too we can a custom function within an I() expression. For
example, we can conditionally estimate the average rank by treating our
estimate as the relative “weight” of each level:

mean_rank <- function(x) {
  weighted.mean(c(3, 4, 5), w = x)
}

ranks <- avg_predictions(mod,
           variables = "cyl",
           hypothesis = ~ I(mean_rank(x)) | cyl)
## 
## Re-fitting to get Hessian
ranks
## 
##  cyl Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
##    4     4.64     0.0885 52.4   <0.001   Inf  4.46   4.81
##    6     4.09     0.1156 35.4   <0.001 908.0  3.86   4.32
##    8     3.12     0.0375 83.2   <0.001   Inf  3.04   3.19
## 
## Type:  probs

These can then be further compared using the hypotheses() function:

hypotheses(ranks, ratio ~ pairwise)
## 
##  Hypothesis Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
##       6 / 4    0.882     0.0296 29.8   <0.001 647.9 0.824  0.940
##       8 / 4    0.672     0.0150 44.7   <0.001   Inf 0.643  0.702
##       4 / 6    1.134     0.0380 29.8   <0.001 647.9 1.059  1.208
##       8 / 6    0.762     0.0234 32.6   <0.001 770.5 0.716  0.808
##       4 / 8    1.487     0.0333 44.7   <0.001   Inf 1.422  1.553
##       6 / 8    1.312     0.0403 32.6   <0.001 770.5 1.233  1.391

vincentarelbundock Jan 23, 2025
Maintainer

Thanks! This is awesome!

I added them to the website.

I'm really quite pleased with this new arbitrary functions in formula solution. It feels very elegant, powerful, and consistent with the rest of the interface.

vincentarelbundock · 2025-01-23T03:18:46Z

vincentarelbundock
Jan 23, 2025
Maintainer

Let's close this discussion now, but feel free to open others to discuss more specific points.

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Complex comparisons #1329

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 5 comments 12 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Functions

Predictions

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Complex comparisons #1329

mattansb Jan 14, 2025

Unit-level (conditional) estimates

{data.table}

{tidyverse}

Average (marginal) estimates.

{data.table}

{tidyverse}

Marginal effects at the mean

{data.table}

{tidyverse}

More…

{data.table}

{tidyverse}

{data.table}

{tidyverse}

{data.table}

{tidyverse}

Joint tests

Conditional

{data.table}

{tidyverse}

Replies: 5 comments · 12 replies

vincentarelbundock Jan 14, 2025 Maintainer

vincentarelbundock Jan 18, 2025 Maintainer

mattansb Jan 18, 2025 Author

vincentarelbundock Jan 18, 2025 Maintainer

mattansb Jan 19, 2025 Author

vincentarelbundock Jan 19, 2025 Maintainer

mattansb Jan 19, 2025 Author

mattansb Jan 22, 2025 Author

Complex contrasts

Unit-level (conditional) estimates

Average (marginal) estimates.

Marginal effects at the mean

More…

conditional contrasts

interaction contrasts

Joint tests

Conditional

vincentarelbundock Jan 22, 2025 Maintainer

Functions

Predictions

mattansb Jan 22, 2025 Author

vincentarelbundock Jan 22, 2025 Maintainer

vincentarelbundock Jan 22, 2025 Maintainer

mattansb Jan 22, 2025 Author

vincentarelbundock Jan 23, 2025 Maintainer

vincentarelbundock Jan 23, 2025 Maintainer

mattansb
Jan 14, 2025

`{data.table}`

`{tidyverse}`

`{data.table}`

`{tidyverse}`

`{data.table}`

`{tidyverse}`

`{data.table}`

`{tidyverse}`

`{data.table}`

`{tidyverse}`

`{data.table}`

`{tidyverse}`

`{data.table}`

`{tidyverse}`

Replies: 5 comments 12 replies

vincentarelbundock
Jan 14, 2025
Maintainer

vincentarelbundock
Jan 18, 2025
Maintainer

mattansb Jan 18, 2025
Author

vincentarelbundock Jan 18, 2025
Maintainer

mattansb Jan 19, 2025
Author

vincentarelbundock Jan 19, 2025
Maintainer

mattansb Jan 19, 2025
Author

mattansb
Jan 22, 2025
Author

vincentarelbundock
Jan 22, 2025
Maintainer

mattansb Jan 22, 2025
Author

vincentarelbundock Jan 22, 2025
Maintainer

vincentarelbundock Jan 22, 2025
Maintainer

mattansb Jan 22, 2025
Author

vincentarelbundock Jan 23, 2025
Maintainer

vincentarelbundock
Jan 23, 2025
Maintainer