Comparison to `marginaleffects` (and a few questions)

vincentarelbundock commented 2 years ago

Hi @JWiley!

I just spent a half day playing with brmsmargins. Awesome, powerful stuff! There are some really neat design ideas in there.

Cards on the table: I was particularly interested in learning about brmsmargins because I am developing a similar R package: marginaleffects. As you know, this is kind of a crowded field, and people are often asking me how the packages relate to each others. So I started writing a vignette with short discussions and side-by-side syntax comparisons for different packages. Obviously, I know that I am fundamentally biased, but I am still trying hard to be as generous as possible in highlighting the benefits of other packages, and in giving them some exposure (we’re all on the same FOSS team!).

I would like to add brmsmargins to that comparisons page. Here’s a preliminary draft of the section:

Rendered HTML: http://www.arelbundock.com/brmsmargins.html
Rmarkdown: https://github.com/vincentarelbundock/marginaleffects/blob/main/sandbox/brmsmargins.Rmd

Note that the notebook requires the development version of marginaleffects:

remotes::install_github("vincentarelbundock/marginaleffects")

I was wondering if you could answer a couple questions I had while trying to replicate the analyses in your vignettes. No pressure if you don’t have time; feel to ignore and close this Issue:

Are there any cool brmsfit features that I ignored in the notebook and that you would like me to highlight?
Do we know why there are small numerical differences between our results in the “integrate out” example?

Cheers!

Vincent

vincentarelbundock commented 2 years ago

apologies for posting, editing, deleting. I noticed a ridiculous mistake in my code and have now updated the notebook.

JWiley commented 2 years ago

@vincentarelbundock I've seen marginaleffects and it looks great and very easy to use.

The main pain point I was trying to solve with brmsmargins was/is accurate point and uncertainty estimates for non linear transforms/link functions, in models with random effects.

In terms of your questions:

I'm not sure how cool it is to most people, but I've been trying to include marginalized coefficients: https://joshuawiley.com/brmsmargins/articles/mixed-effects-marginaleffects.html#marginal-coefficients based on: Hedeker, D., du Toit, S. H., Demirtas, H., & Gibbons, R. D. (2018). A note on marginalization of regression parameters from mixed models of binary outcomes. Biometrics, 74(1), 354-361.
I suspect the integrate out examples are the trickiest to match. So when integrating out REs to generate a marginal estimate, we basically need the posterior prediction based on the fixed effects and then to integrate over the (assumed) random effect distribution. I use numerical monte carlo integration because it makes it relatively easy to have multiple random effects and (down the road) for those not to follow a (multivariate) normal distribution. I'm particularly interested in supporting the t distribution, which brms already does to allow for heavier tailed REs. Anyway, my steps are, for each posterior draw are:

Generate the prediction on the link scale based on the fixed effects of the model
Take the estimated covariance matrix for all random effect blocks (i.e., in a random intercept only model, this is brms estimate of the variance of the intercept, not the estimated intercepts for specific people) --- this is important, it is not the same as generating predictions from brms using the REs, that will be in sample only and does not give the right estimate and uncertainty for the population marginal value.
Draw k random samples from the random effect blocks assumed distribution based on the distribution + covariance matrix
Add all k random samples to the fixed effects prediction
Back transform to the response scale
Average

In terms of the vignette you showed, one minor point, brmsmargins defaults to 99% CIs whereas I believe marginaleffects defaults to 95% CIs. There are a few points where you set CI = 0.95 but also a few where you do not (e.g., AMEs: Grand Mean).

I could not see how you integrate out REs in your package. Here is a small simulation study. I use the same marginaleffects code you did in your Rmd / html file under the integrate out RE section. You'll notice that while both packages have very similar and just 1-2 percent bias relative to the "true" synthetic population value, our CIs are markedly different. The average width of the CI interval is 0.124 for marginaleffects and 0.320 for brmsmargins. The synthetic population value is included in 99.6% of the 500 simulation runs for brmsmargins and 84.6% of the simulation runs for marginaleffects. Neither is quite perfect as in both cases, we are aspiring to a 95% interval, although it is only 500 simulation runs so some variation is expected.

library(brms)
library(brmsmargins)
library(data.table)
library(marginaleffects)

library(future)
library(doFuture)
library(foreach)
library(doRNG)

h <- 5e-5

set.seed(1234)
x.b0 <- rnorm(10000, mean = 0, sd = 3)
x.b1 <- rnorm(10000, mean = 2, sd = 3)

mu.0 <- x.b0 + x.b1 * 0
mu.1 <- x.b0 + x.b1 * 1

y.0a <- rbinom(n = 10000, size = 1, prob = plogis(mu.0))
y.0b <- rbinom(n = 10000, size = 1, prob = plogis(mu.0))
y.1a <- rbinom(n = 10000, size = 1, prob = plogis(mu.1))
y.1b <- rbinom(n = 10000, size = 1, prob = plogis(mu.1))

d <- data.table(
  ID = rep(1:10000, times = 4),
  y = c(y.0a, y.0b, y.1a, y.1b),
  x = rep(c(0, 0, 1, 1), each = 10000))
d[, x := factor(x)]

d[, .(Diff = mean(y[x == 1]) - mean(y[x == 0])), by = ID][, .(MDiff = mean(Diff))]
##     MDiff
##1: 0.16675

set.seed(1234)
IDs <- sample(1:10000, size = 100, replace = FALSE)

d <- as.data.frame(d)

m1 <- brm(y ~ x + (1 + x | ID),
          data = d[d$ID %in% IDs, ], chains = 2, cores = 2,
          family = bernoulli(),
          iter = 4000, refresh = 0, silent = 2,
          backend = "cmdstanr")
summary(m1)

ame1 <- brmsmargins(m1,
                    at = data.frame(x = c("0", "1")),
                    contrasts = cbind("AME x" = c(-1, 1)),
                    CI = .95, CIType = "ETI",
                    effects = "integrateoutRE", k = 20L)
ame1$ContrastSummary

ame2 <- marginaleffects::marginaleffects(m1)
ame2sum <- as.data.table(summary(ame2, FUN = mean))

registerDoFuture()
plan(multisession, workers = 20)

set.seed(12345)
res <- foreach(i = 1:500, .combine = c) %dorng% {
  IDs <- sample(1:10000, size = 100, replace = FALSE)
  mfit <- update(m1, newdata = d[d$ID %in% IDs, ], refresh = 0, silent = 2)

  ame <- brmsmargins(mfit,
                     at = data.frame(x = c("0", "1")),
                     CI = .95, CIType = "ETI",
                     contrasts = cbind("AME x" = c(-1, 1)),
                     effects = "integrateoutRE", k = 20L)

  ame.use <- as.data.table(ame$ContrastSummary[, c("M", "LL", "UL")])

  amealt <- marginaleffects::marginaleffects(mfit)
  amealt.use <- as.data.table(summary(amealt, FUN = mean))[, .(estimate, conf.low, conf.high)]

  colnames(amealt.use) <- colnames(ame.use) <- c("Est", "LL", "UL")

  list(
    rbind(
      cbind(Type = "brmsmarginal", run = i, as.data.frame(ame.use)),
      cbind(Type = "marginaleffects", run = i, as.data.frame(amealt.use))))https://github.com/vincentarelbundock
}

test <- as.data.table(do.call(rbind, res))

true.diff <- as.data.table(d)[, .(Diff = mean(y[x == 1]) - mean(y[x == 0])), by = ID][, mean(Diff)]

test[, (mean(Est) - true.diff) / true.diff, by = Type]
##              Type          V1
##1:    brmsmarginal -0.02695898
##2: marginaleffects -0.01106552

test[, mean(LL <= true.diff & UL >= true.diff), by = Type]
##              Type    V1
##1:    brmsmarginal 0.996
##2: marginaleffects 0.846

test[, .(CIWidth = mean(UL - LL)), by = Type]
##              Type   CIWidth
##1:    brmsmarginal 0.3201659
##2: marginaleffects 0.1243733

Now let's re-run the simulation. This time I just use 20 simulations, but rather than sampling 100 people at a time, I sample 500 people at a time. Of course 20 simulations is very small. The trend it shows though is again similar, minimal bias in average estimate. The CI width for marginaleffects is now much smaller, the CI width for brmsmargins slightly smaller. Coverage will be noisy in just 20 simulations but is 100% for brmsmargins and 75% for marginaleffects.


set.seed(12345)
resbig <- foreach(i = 1:20, .combine = c) %dorng% {
  IDs <- sample(1:10000, size = 500, replace = FALSE)
  mfit <- update(m1, newdata = d[d$ID %in% IDs, ], refresh = 0, silent = 2)

  ame <- brmsmargins(mfit,
                     at = data.frame(x = c("0", "1")),
                     CI = .95, CIType = "ETI",
                     contrasts = cbind("AME x" = c(-1, 1)),
                     effects = "integrateoutRE", k = 20L)

  ame.use <- as.data.table(ame$ContrastSummary[, c("M", "LL", "UL")])

  amealt <- marginaleffects::marginaleffects(mfit)
  amealt.use <- as.data.table(summary(amealt, FUN = mean))[, .(estimate, conf.low, conf.high)]

  colnames(amealt.use) <- colnames(ame.use) <- c("Est", "LL", "UL")

  list(
    rbind(
      cbind(Type = "brmsmarginal", run = i, as.data.frame(ame.use)),
      cbind(Type = "marginaleffects", run = i, as.data.frame(amealt.use))))
}

testbig <- as.data.table(do.call(rbind, resbig))

testbig[, (mean(Est) - true.diff) / true.diff, by = Type]
##              Type         V1
##1:    brmsmarginal 0.01150610
##2: marginaleffects 0.01135871

testbig[, mean(LL <= true.diff & UL >= true.diff), by = Type]
##              Type   V1
##1:    brmsmarginal 1.00
##2: marginaleffects 0.75

testbig[, .(CIWidth = mean(UL - LL)), by = Type]
##              Type    CIWidth
##1:    brmsmarginal 0.29237039
##2: marginaleffects 0.05471667

vincentarelbundock commented 2 years ago

@JWiley thank you so much for taking the time to write this! It is tremendously helpful.

Yes, .99 vs. .95. explained most of the discrepancies.

I could not see how you integrate out REs in your package

A fair question! Especially since I do not integrate out REs in the package. Your transparent and useful description of the procedure also makes it clear that, for now at least, this is really outside the scope of my package.

I have added your package to my comparison page with a selection of models I was able to match exactly, and a note highlighting what seem to be two areas where brmsmargins offers features which are not available at all in marginaleffects: (1) marginal coefficients, and (2) integrating out the REs.

Thanks again for clarifying all this! I learned a lot from your post and from reading through parts of your code. (This is fun.)

Edit link: https://vincentarelbundock.github.io/marginaleffects/articles/alternative_software.html

JWiley / brmsmargins

Comparison to `marginaleffects` (and a few questions) #8