Please kindly help me accomplish this:
I want to get the true average marginal effects using the function marginal_effects in brms that in fact marginalize (not condition) on the covariates that are not the treatment.
Specifically, I want to get the following using brms and marginal_effects:
"(…), average fitted values, calculates the value ˆ y for every case in the
data and averages the resulting fitted values. The interpretation of this quantity is as
the average (or typical) outcome we would expect to observe were our model an accurate
representation of the data-generating process for the outcome.
(…) AMEs calculate marginal effects at every observed value of X and average across the resulting effect estimates. " as described here by Thomas J. Leeper.
This is what the documentation of brms says about marginal_effects:
When creating marginal_effects for a particular predictor (or interaction of two predictors), one has to choose the values of all other predictors to condition on. By default, the mean is used for continuous variables and the reference category is used for factors, but you may change these values via argument conditions. (…) Since we condition on rather than actually marginalizing variables, the name marginal_effects is possibly not ideally chosen in retrospect.
The option “to choose the values of all other predictors to condition on” is ok, but I would like to additionally be able to marginalize (not condition) on all the other covariates that I did not condition on (those that I did not choose their values), instead of doing this “[b]y default, the mean is used for continuous variables and the reference category is used for factors” for those variables whose values I did not specify.
Here are the sample data and code:
dt = read.table(header = TRUE, text = "
n r r/n group treat c2 c1 w
62 3 0.048387097 1 0 0.1438 1.941115288 1.941115288
96 1 0.010416667 1 0 0.237 1.186583128 1.186583128
17 0 0 0 0 0.2774 1.159882668 3.159882668
41 2 0.048780488 1 0 0.2774 1.159882668 3.159882668
212 170 0.801886792 0 0 0.2093 1.133397521 1.133397521
143 21 0.146853147 1 1 0.1206 1.128993008 1.128993008
143 0 0 1 1 0.1707 1.128993008 2.128993008
143 33 0.230769231 0 1 0.0699 1.128993008 1.128993008
73 62 1.260273973 0 1 0.1351 1.121927228 1.121927228
73 17 0.232876712 0 1 0.1206 1.121927228 1.121927228")
mod<-brm(r | trials(n) + weights(w) ~ treat*c2+(1|group), data=dt, family=binomial(link=logit))
Thank you in advance for any help.