Monotonic effects bug

[DominiqueMakowski]

While we fixed get_datagrid to return a compatible dataframe (https://github.com/easystats/insight/issues/520), it turns out that the next bug happens in marginaleffects which treats the variable as continuous :)

library(brms)

model <- brm(mpg ~ mo(carb), data = mtcars, refresh = 0)
#> Compiling Stan program...
#> Start sampling

marginaleffects::marginaleffects(model, variable = "carb", insight::get_datagrid(model))
#> Error: Invalid values in variable 'carb': '1.00001', '2.00001', '3.00001', '4.00001', '6.00001', '8.00001'

^{Created on 2022-03-08 by the reprex package (v2.0.1)}

[vincentarelbundock]

Can you transform carb to a factor before fitting the model?

[vincentarelbundock]

@DominiqueMakowski, the devel version of marginaleffects now detects the mo() function in the formula and automatically treats the variable as categorical instead of continuous, that is, it computes the pairwise contrast instead of dydx:

library(brms)
library(marginaleffects)
mod <- brm(mpg ~ mo(carb), data = mtcars, refresh = 0, backend = "cmdstanr")
mfx <- marginaleffects(mod, insight::get_datagrid(mod))
summary(mfx)

    ## Average marginal effects 
    ##   Term Contrast  Effect   2.5 %   97.5 %
    ## 1 carb    2 - 1  -1.728  -4.924 -0.00671
    ## 2 carb    3 - 1  -5.090  -9.207 -0.62795
    ## 3 carb    4 - 1  -6.931 -11.477 -2.73500
    ## 4 carb    6 - 1  -9.133 -14.483 -3.60710
    ## 5 carb    8 - 1 -12.027 -19.240 -4.86427
    ## 
    ## Model type:  brmsfit 
    ## Prediction type:  response

[DominiqueMakowski]

Thanks!

But I wonder though, isn't the reason why people would use monotonic effects to still kinda interpret it as a (or model is as a pseudo) "continuous" (in the sense of ordered)? As far as I understand having dydx instead of contrasts against the lowest level would be what people would expect no? Or at least contrasts like 2-1, 3-2, 4-3 etc.

Tagging people that might have some more experience with monotonic @mattansb @bwiernik

[vincentarelbundock]

Interesting. Will be curious to hear what people think.

A real slope (dydx) won't be possible, because predict.brmsfit() doesn't work at x+eps. But perhaps a different contrast would make sense, and I'm very interested in implementing new contrasts.

[bwiernik]

sequential contrasts would make the most sense to me here

[vincentarelbundock]

I plan to create a new function specifically to estimate contrasts. It will be a bit more flexible and allow sequential contrasts.

What do you think of differences() as a function name?

I would like to avoid collisions with stats::contrasts and emmeans::contrast, but I also don't want to confuse users...

[bwiernik]

The MASS package has a subsequent diffs contrast that could be emulated https://www.rdocumentation.org/packages/MASS/versions/7.3-55/topics/contr.sdif

[bwiernik]

Differences or comparisons would make sense to me.

[vincentarelbundock]

I like comparisons a lot! Thanks!

[mattansb]

I agree that any mo() should be treated as discrete / a factor.

As for the best contrasts - either compare to baseline (for a cumulative effect) L1 - baseline, L2 - baseline, or sequential as Dom suggested.

comparisons() is great!

[vincentarelbundock]

Hi all,

I pushed a new comparisons() function to Github which supports a bunch of different contrasts for factors and numeric variables. This vignette illustrates how to use it:

https://vincentarelbundock.github.io/marginaleffects/articles/contrasts.html

I also paste @DominiqueMakowski's example below. Let me know if you get a chance to try it. It's still experimental, so I'm very eager to get feedback. But of course I know that everyone's busy, so please feel free to ignore this.

library(brms)
library(marginaleffects)
mod <- brm(mpg ~ hp + mo(carb), data = mtcars, backend = "cmdstanr")

comparisons(mod, newdata = datagrid())

##    rowid term contrast  comparison  conf.low conf.high       hp carb
## 1:     1   hp       +1 -0.07843755 -0.109624 -0.046889 146.6875    4
## 2:     1 carb    2 - 1  0.24414601 -1.137128  2.139320 146.6875    4
## 3:     1 carb    3 - 1  0.63031393 -2.221756  3.179065 146.6875    4
## 4:     1 carb    4 - 1  1.06382719 -2.964875  4.093695 146.6875    4
## 5:     1 carb    6 - 1  2.93805804 -4.856541  9.101595 146.6875    4
## 6:     1 carb    8 - 1  5.38878895 -6.303341 16.868125 146.6875    4

comparisons(mod,
            contrast_numeric = c(90, 110),
            contrast_factor = "sequential",
            newdata = datagrid())

##    rowid term contrast comparison   conf.low conf.high       hp carb
## 1:     1   hp 110 - 90 -1.5687510 -2.1924800 -0.937780 146.6875    4
## 2:     1 carb    2 - 1  0.2441460 -1.1371275  2.139320 146.6875    4
## 3:     1 carb    3 - 2  0.1879024 -1.1767482  1.818360 146.6875    4
## 4:     1 carb    4 - 3  0.2205211 -0.9333101  1.977351 146.6875    4
## 5:     1 carb    6 - 4  1.4909809 -1.9544407  6.450277 146.6875    4
## 6:     1 carb    8 - 6  1.7184008 -2.1166329  9.764168 146.6875    4

comparisons(mod,
            contrast_numeric = "iqr",
            contrast_factor = "pairwise") |>
    tidy()

##    term contrast   estimate   conf.low conf.high
## 1  carb    2 - 1  0.2441460 -1.1371275  2.139320
## 2  carb    3 - 1  0.6303139 -2.2217556  3.179065
## 3  carb    3 - 2  0.1879024 -1.1767482  1.818360
## 4  carb    4 - 1  1.0638272 -2.9648747  4.093695
## 5  carb    4 - 2  0.5990961 -2.0341544  2.905629
## 6  carb    4 - 3  0.2205211 -0.9333101  1.977351
## 7  carb    6 - 1  2.9380580 -4.8565411  9.101595
## 8  carb    6 - 2  2.4501333 -3.8982442  8.137011
## 9  carb    6 - 3  1.9847614 -2.8463971  7.274665
## 10 carb    6 - 4  1.4909809 -1.9544407  6.450277
## 11 carb    8 - 1  5.3887889 -6.3033411 16.868125
## 12 carb    8 - 2  4.8837559 -5.5709082 15.754891
## 13 carb    8 - 3  4.4405691 -4.1588337 15.162958
## 14 carb    8 - 4  3.9069607 -3.6354205 13.978538
## 15 carb    8 - 6  1.7184008 -2.1166329  9.764168
## 16   hp      IQR -6.5495354 -9.1536040 -3.915232