Stan_glmer Marginal Effects

Modeling
1

I am fitting a multi-level model in stan_glmer (rstanarm package). It is a logistic regression that predicts a binary outcome given two individual-level predictors and one variable with varying intercepts (50 levels).

Var1 is binary and var2 is continuous.

mod <- stan_glmer(outcome~var1+var2 + (1|state), data=d, binomial(link = ‘logit’),
prior_intercept = normal(0,1))

Is there a function that allows me to estimate marginal effects for var1 or var2 without making predictions for each state and then weighting the results by state observations? This would be especially helpful when I add more complexity to the model.

Thanks

2

I am not sure which meaning you are ascribing to “marginal effects”. If you mean, posterior predictions irrespective of which state the person lives, then you can use the syntax

posterior_predict(mod, re.form = ~ 0)

and similarly for posterior_linpred.

3

Thanks for the quick reply. What I meant by marginal effects is the effect of either var1 or var2 when everything else is held constant at some level. For a more straightforward example, without the state-varying intercept, I would hold var2 at its mean/median and then predict var1 based on the modeled association between var1 and outcome.

For the model with the state variable, I am not interested in predicting any particular state; rather, I want to essentially hold that constant, which usually requires (at least for the frequentist approach) making predictions on all levels of a factor and weighting by the number of state observations. A good example is the R package effects (e.g. effects(“var1”,mod)).

The tedious component comes into play when I have variables with many levels. Does posterior_predict(mod, re.form = ~0) accomplish that?

4

OK, the Stata definition of marginal effects is different from posterior_predict(mod, re.form = ~ 0). For that, you need to do

X <- model.matrix(mod)[ , c("var1", "var2")]
beta <- as.matrix(mod)[ , c("var1", "var2")]
Pr <- posterior_linpred(mod, transform = TRUE)
ME_var1 <- (beta[ , 1] / X[ , 1]) * Pr * (1 - Pr)
ME_var2 <- (beta[ , 2] / X[ , 2]) * Pr * (1 - Pr)
5

Great–thanks again for the help. I appreciate it.

6

I meant to say

ME_var1 <- beta[ , 1] * sweep(Pr * (1 - Pr), MARGIN = 2, STATS = X[ , 1], FUN = `/`)
ME_var2 <- beta[ , 2] * sweep(Pr * (1 - Pr), MARGIN = 2, STATS = X[ , 2], FUN = `/`)
7

Also relevant is this paper:
http://www.stat.columbia.edu/~gelman/research/published/ape17.pdf
The short answer is that there’s no unique way to frame the problem. In that paper we talk about what seems to me to be a natural framing, and it gets complicated!

2 Likes
8

Thank you for sending the article Prof Gelman–it was really informative.

bgoodri: for the ME_var1 (binary variable), did you mean the following:

ME_var1 <- beta[ , 1] * sweep(Pr * (1 - Pr), MARGIN = 2, STATS = X[ , 1], FUN = -)

Otherwise, I am getting a lot of -Inf.

9

That was assuming the predictors are continuous. If the predictor is binary (it is is a stretch to refer to things as “marginal”), I would say to just calculate the posterior distribution of the difference in probability

df <- model.frame(mod)
df$var1 <- 0
Pr_0 <- posterior_linpred(mod, data = df, transform = TRUE)
df$var1 <- 1
Pr_1 <- posterior_linpred(mod, data = df, transform = TRUE)
ME_var1 <- Pr_1 - Pr_0
1 Like