Results (snapshots) have changed in _marginaleffects_ 0.15.0

[strengejacke]

Hi Vincent,
I just saw that some tests are failing for ggeffects, and it turned out that something has changed since the latest marginaleffects update to 0.15.0 - not sure if this was a bug fix or is a new bug.

See here results from version 0.15.0:

library(ggeffects)

set.seed(123)
dat <- data.frame(
  outcome = rbinom(n = 100, size = 1, prob = 0.35),
  var_binom = as.factor(rbinom(n = 100, size = 1, prob = 0.2)),
  var_cont = rnorm(n = 100, mean = 10, sd = 7),
  groups = sample(letters[1:4], size = 100, replace = TRUE)
)

m <- glm(outcome ~ var_binom * var_cont + groups,
  data = dat, family = binomial()
)

# new results since marginal effects 0.15.0
hypothesis_test(m, "var_binom", scale = "exp")
#> # Pairwise comparisons
#> 
#> var_binom | Contrast |       95% CI |     p
#> -------------------------------------------
#> 0-1       |     1.62 | [1.26, 2.10] | 0.901
#> 
#> Contrasts are presented as odds ratios.

# former results in marginal effects 0.14.x
#>    var_binom | Contrast |       95% CI |     p
#>    -------------------------------------------
#>    0-1       |     0.99 | [0.80, 1.22] | 0.902

d <- data_grid(m, "var_binom")
marginaleffects::predictions(
  m, newdata = d, by = "var_binom", hypothesis = "pairwise", transform = "exp"
)
#> 
#>   Term Estimate Pr(>|z|)   S 2.5 % 97.5 %
#>  0 - 1     1.62    0.901 0.2  1.26    2.1
#> 
#> Columns: term, estimate, p.value, s.value, conf.low, conf.high 
#> Type:  invlink(link)

packageVersion("marginaleffects")
#> [1] '0.15.0'

^{Created on 2023-09-15 with reprex v2.0.2}

And here for version 0.14.0:

# old results with marginal effects 0.14.0
hypothesis_test(m, "var_binom", scale = "exp")
#> # Pairwise comparisons
#> 
#> var_binom | Contrast |       95% CI |     p
#> -------------------------------------------
#> 0-1       |     0.99 | [0.80, 1.22] | 0.902
#> 
#> Contrasts are presented as odds ratios.

d <- data_grid(m, "var_binom")
marginaleffects::predictions(
  m, newdata = d, by = "var_binom", hypothesis = "pairwise", transform = "exp"
)
#> 
#>   Term Estimate Pr(>|z|)   S 2.5 % 97.5 %
#>  0 - 1    0.987    0.902 0.1 0.799   1.22
#> 
#> Columns: term, estimate, p.value, s.value, conf.low, conf.high

packageVersion("marginaleffects")
#> [1] '0.14.0'

Could you clarify whether the changes are intentional or a bug?

[vincentarelbundock]

This is not exactly a bug. It’s related to an improvement in functionality, but I did not expect it to break existing code. Sorry about that!

The error was not triggered when I checked rev-deps before releasing to CRAN. Is this snapshot in your test suite?

I have now added a snapshot test to my suite to make sure this doesn’t change: de7f8b6

Pre-0.15.0, the predictions() function would make predictions on the link scale and backtransform to the response scale if these three conditions held:

1. Suppported model (glm, gam, and a couple more)
2. type = NULL
3. transform = NULL

As soon as you called transform="exp", we would predict immediately on the response scale by flipping the argument to type="response" and not do any backtransformation.

This implicit behavior was hard to document and explain. In 0.15.0, I improved it by introducing a new explicit type = "invlink(link)". A side benefit of this approach is that we can now transform the backtransformed estimates, which wasn’t possible before.

This means that you now get different results if you use type="invlink(link)" (the default) or type="response":

library(marginaleffects)
set.seed(123)
dat <- data.frame(
  outcome = rbinom(n = 100, size = 1, prob = 0.35),
  var_binom = as.factor(rbinom(n = 100, size = 1, prob = 0.2)),
  var_cont = rnorm(n = 100, mean = 10, sd = 7),
  groups = sample(letters[1:4], size = 100, replace = TRUE)
)
m <- glm(outcome ~ var_binom * var_cont + groups,
  data = dat, family = binomial()
)

d <- structure(list(var_binom = structure(1:2, levels = c("0", "1"
), class = "factor"), var_cont = c(9.24717241397544, 9.24717241397544
), groups = structure(c(1L, 1L), levels = "b", class = "factor")), class = "data.frame", row.names = c(NA, 
-2L))

# type = NULL or type = "invlink(link)"
predictions(
  m, newdata = d, by = "var_binom", hypothesis = "pairwise", transform = "exp"
)$estimate
# [1] 1.61968

# type = "response"
predictions(
  m, newdata = d, by = "var_binom", hypothesis = "pairwise", transform = "exp", type = "response"
)$estimate
# [1] 0.9867827

[strengejacke]

The error was not triggered when I checked rev-deps before releasing to CRAN. Is this snapshot in your test suite?

Yes, but snapshot tests are not run on CRAN, and furthermore, I think that suggested packages (marginaleffects is suggested by ggeffects) are ignored in rev-dep checks (i.e. reverse-suggests are not run for your package - I think).

Thank you for clarifying! One question remains: which of the options now returns contrasts on the odds ratio scale, and which on the probability (response) scale?

[vincentarelbundock]

Neither. Both options exponentiate a response scale. The transform argument is always the very last thing that the package's functions do. So what you are doing here:

1. Compute two predictions on the response (probability) scale. In the two commands, these predictions are arrived at slightly differently: directly on the response scale or via backtransformation. However, both commands produce the same scale of predictions.
2. Take the pairwise difference between your predictions.
3. Apply the transform argument (exponentiation) directly to the pairwise difference at the very last minute.

I'm not sure if that's what you want...

[strengejacke]

ok, got it, how transform works. But what is the default if response = NULL?

set.seed(1234)
dat <- data.frame(
  outcome = rbinom(n = 100, size = 1, prob = 0.35),
  x1 = as.factor(sample(1:3, size = 100, TRUE, prob = c(0.5, 0.2, 0.3))),
  x2 = rnorm(n = 100, mean = 10, sd = 7)
)

m <- glm(outcome ~ x1 + x2, data = dat, family = binomial())
d <- ggeffects::data_grid(m, "x1")

# what is this scale?
marginaleffects::predictions(
  m, newdata = d, by = "x1", hypothesis = "pairwise"
)$estimate
#> [1] 0.6268201 0.3798960 0.2672561

# contrasts as log-odds?
marginaleffects::predictions(
  m, newdata = d, by = "x1", hypothesis = "pairwise", type = "link"
)$estimate
#> [1]  0.5185991 -0.4899896 -1.0085887

# contrasts as probabilities?
marginaleffects::predictions(
  m, newdata = d, by = "x1", hypothesis = "pairwise", type = "response"
)$estimate
#> [1]  0.08277724 -0.10101103 -0.18378827

^{Created on 2023-09-15 with reprex v2.0.2}

[strengejacke]

In other words: If I now want the odds ratio scale, I would use

marginaleffects::predictions(
  m, newdata = d, by = "x1", hypothesis = "pairwise", type = "link", transform = "exp"
)

right?

[vincentarelbundock]

This command

marginaleffects::predictions(
  m, newdata = d, by = "x1", hypothesis = "pairwise", type = "link", transform = "exp"
)

gives you this:

$e^{(\beta_0 + \beta_1 \cdot 1) - (\beta_0 + \beta_1 \cdot 0)}$

This command

marginaleffects::predictions(
  m, newdata = d, by = "x1", hypothesis = "pairwise", type = "response", transform = "exp"
)

gives you this:

$e^{P(Y=1|X=1)-P(Y=1|X=0)}$

I'll need to keep this issue open to check what happens when type="invlink(link)" (the default), because I'm not sure if the invlink is applied to the individual estimates or to the difference (this would be a problem).

Unfortunately I don't have time to check now.

[strengejacke]

Ok, and what is the default for logistic regression when neither is provided?

# what is this scale?
marginaleffects::predictions(
  m, newdata = d, by = "x1", hypothesis = "pairwise"
)$estimate

I think then I have everything I need to know :-)

[vincentarelbundock]

I think you made me find a really major problem. I'll look at it tonight and will probably have to make an emergency release.

[vincentarelbundock]

OK, this should be sorted on Github. I will release a new version to CRAN as soon as possible. See the NEWS for a description of the bug.

Thanks again for the report!

Note that a warning is now issued in some cases. That warning can be suppressed by choosing the type value explicitly instead of letting it default to invlink(link).

library(marginaleffects)
set.seed(123)
dat <- data.frame(
  outcome = rbinom(n = 100, size = 1, prob = 0.35),
  var_binom = as.factor(rbinom(n = 100, size = 1, prob = 0.2)),
  var_cont = rnorm(n = 100, mean = 10, sd = 7),
  groups = sample(letters[1:4], size = 100, replace = TRUE)
)
m <- glm(outcome ~ var_binom * var_cont + groups,
  data = dat, family = binomial()
)

d <- structure(list(var_binom = structure(1:2, levels = c("0", "1"
), class = "factor"), var_cont = c(9.24717241397544, 9.24717241397544
), groups = structure(c(1L, 1L), levels = "b", class = "factor")), class = "data.frame", row.names = c(NA, 
-2L))

predictions(
  m, newdata = d, by = "var_binom", hypothesis = "pairwise", transform = "exp"
)$estimate
# Warning: The `type="invlink"` argument is not available unless `hypothesis` is
#   `NULL` or a single number. The value of the `type` argument was changed
#   to "response" automatically. To suppress this warning, use
#   `type="response"` explicitly in your function call.
# [1] 0.9867827

predictions(
  m, newdata = d, by = "var_binom", hypothesis = "pairwise", transform = "exp", type = "invlink(link)"
)$estimate
# Warning: The `type="invlink"` argument is not available unless `hypothesis` is
#   `NULL` or a single number. The value of the `type` argument was changed
#   to "response" automatically. To suppress this warning, use
#   `type="response"` explicitly in your function call.
# [1] 0.9867827

predictions(
  m, newdata = d, by = "var_binom", hypothesis = "pairwise", transform = "exp", type = "response",
)$estimate
# [1] 0.9867827

[strengejacke]

Great! To summarize:

• type = NULL returns on the response scale (e.g. probabilities for logistic regression)
• Else, type returns on the scale as defined
• transform applies transformation for predictions, based on what is returned by type

right?

As an example: for logistic regression, type = "link" and transform = "exp" would return on the odds ratio scale?

[vincentarelbundock]

transform applies to the final estimate after executing the hypothesis test. So if hypothesis="pairwise", we compute predicted probabilities, then take the pairwise differences, and only then apply transform="exp".

But you're right: The difference in two link scale predictions made by incrementing a single variable is the corresponding regression coefficient, and exponentiating a logit coefficient gives us an odds ratio.

[strengejacke]

btw, the changes in marginaleffects did not really "break" any code on my side, they just returned different results. So there's no need for me to update ggeffects as well, because results will be correct when the current dev of marginaleffects is installed. I just need to update the snapshot tests sometime...

[vincentarelbundock]

Yeah, but the 1.62 result is incorrect...

[strengejacke]

But that result should automatically change (to the correct one) once marginaleffects is updated, right?

[vincentarelbundock]

Absolutely!

[vincentarelbundock]

on cran now