In the past, I have used the multcomp package to compute inferences for my statistical models in R. However, the marginaleffects package seems to be the new kid in town, and I wanted to learn how it works. In this post, I tried to familiarize myself with the marginaleffects syntax to compute different statistical estimands for multilevel models fitted with the lme4 and brms packages.
1 Simulate multilevel data
I will use simulated data from a Generalized Linear Mixed Model (GLMM) in this post. The simulation code was inspired by the online documentation of the faux R package. For this example, imagine that n_subjects subjects respond to n_items stimuli in a diagnostic decision task. The binary response variable y_bin reflects whether the diagnostic decision is correct or not. For some trials, the participants are presented with advice (advice_present = 1) that should help them making the correct diagnostic decision.
simulate <- function(n_subjects = 100, n_items = 50,
b_0 = 0.8, b_a = 1,
sd_u0s = 0.5, sd_u0i = 0.5, sd_u1s = 0.5, ...){
require(dplyr)
require(faux)
# simulate design
dat <- add_random(subject = n_subjects, item = n_items) %>%
mutate(advice_present = rbinom(n(), 1, prob = 2/3)) %>%
# add random effects
add_ranef("subject", u0s = sd_u0s) %>%#
add_ranef("subject", u1s = sd_u1s) %>%
add_ranef("item", u0i = sd_u0i) %>%
# compute dependent variable
mutate(linpred = b_0 + u0i + u0s +
(b_a + u1s) * advice_present) %>%
mutate(y_prob = plogis(linpred)) %>%
mutate(y_bin = rbinom(n = n(), size = 1, prob = y_prob))
dat
}
set.seed(1)
dat <- simulate()2 Fit multilevel models
I will fit a multilevel model with both the lme4 and the brms package. The specified model is equal to the true model from which the data have been simulated.
f <- y_bin ~ 1 + advice_present +
(1 + advice_present || subject) + (1|item)
fit_lme4 <- glmer(f, data = dat, family = "binomial")
set.seed(1)
fit_brms <- brm(f, data = dat, family = "bernoulli",
backend = "cmdstanr",
chains = 4, cores = 4)Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2 finished in 35.7 seconds.
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 finished in 35.8 seconds.
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 finished in 37.4 seconds.
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1 finished in 38.2 seconds.
All 4 chains finished successfully.
Mean chain execution time: 36.8 seconds.
Total execution time: 38.3 seconds.summary(fit_lme4)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: y_bin ~ 1 + advice_present + (1 + advice_present || subject) +
(1 | item)
Data: dat
AIC BIC logLik deviance df.resid
4997.9 5030.5 -2494.0 4987.9 4995
Scaled residuals:
Min 1Q Median 3Q Max
-5.4937 0.2292 0.3969 0.5510 1.4685
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 0.2173 0.4661
subject.1 advice_present 0.2519 0.5018
item (Intercept) 0.1949 0.4415
Number of obs: 5000, groups: subject, 100; item, 50
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.75270 0.09540 7.89 3.03e-15 ***
advice_present 1.02600 0.09085 11.29 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
advic_prsnt -0.339summary(fit_brms) Family: bernoulli
Links: mu = logit
Formula: y_bin ~ 1 + advice_present + (1 + advice_present || subject) + (1 | item)
Data: dat (Number of observations: 5000)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~item (Number of levels: 50)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.46 0.06 0.35 0.60 1.00 1318 1884
~subject (Number of levels: 100)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.48 0.07 0.36 0.62 1.01 1248 2241
sd(advice_present) 0.53 0.10 0.33 0.72 1.00 899 1800
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.75 0.10 0.56 0.94 1.00 1682 2416
advice_present 1.03 0.10 0.85 1.23 1.00 2441 2573
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).When looking at the model outputs, we can see that both models have very similar parameter estimates, and parameter estimates closely match the true values we specified in the simulation. The diagnostics of the brms model (Rhats < 0, decent ESS, and no divergent transitions) suggest that the model is identified and has successfully converged, which is what we expect when fitting the true model to a large enough sample.
3 Estimate different contrasts with marginaleffects
In the following sections, I compute predictions from these models with the marginaleffects package and estimate contrasts for different estimands. There are usually several different ways how to compute the same estimates with marginaleffects, and I will convince myself that they produce similar results. Great resources on these topics are the documentation of the marginaleffects package, as well as the excellent blog-posts (1, 2, 3) by Andrew Heiss.
3.1 Some important function arguments and options in marginaleffects
When using the datagrid function inside one of the many marginaleffects functions, be aware that by default all variables not explicitly specified are set to the mean or mode (depending on the variable type). Better check the result of datagrid to make sure you know which predictor values your predictions are actually based on!
When using type = response in marginaleffects for multilevel models fitted with lme4 or brms, this will produce estimates for the conditional expected value of the response \(E(Y|x, u)\) (and does not simulate individual response values from the posterior predictive distribution).
By default, marginaleffects averages the posterior draws of brms models using the median. However, we might prefer using the mean to assure that the order of aggregation does not matter.
The other option specifies the type of posterior intervals (equal-tailed intervals vs. highest density intervals). We use the default "eti" but list this option here as a reminder that it exist.
3.2 Treatment effect for an average person and an average item
Estimand:
\[ \begin{aligned} & P(Y = 1 | advice\_present = 1, u_{0s} = 0, u_{1s} = 0, u_{0i} = 0) \\ & \quad - P(Y = 1 | advice\_present = 0, u_{0s} = 0, u_{1s} = 0, u_{0i} = 0) \end{aligned} \]
3.2.1 lme4
The option re.form = NA specifies that all random effects are set to 0 when computing predictions. The option re.form = NULL specifies that all random effects are always included.
avg_predictions(fit_lme4,
variables = list(advice_present = 0:1),
re.form = NA, type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))
Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
b2-b1=0 0.176 0.0166 10.6 <0.001 84.8 0.143 0.208
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response avg_comparisons(fit_lme4,
variables = "advice_present",
re.form = NA, type = "response")
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 %
advice_present mean(1) - mean(0) 0.176 0.0166 10.6 <0.001 84.8 0.143
97.5 %
0.208
Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
Type: response predictions(fit_lme4,
newdata = datagrid(advice_present = 0:1),
re.form = NA, type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))
Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
b2-b1=0 0.176 0.0166 10.6 <0.001 84.8 0.143 0.208
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response comparisons(fit_lme4,
newdata = datagrid(advice_present = 1),
re.form = NA, type = "response")
Term Contrast advice_present Estimate Std. Error z Pr(>|z|) S
advice_present 1 - 0 1 0.176 0.0166 10.6 <0.001 84.8
2.5 % 97.5 % subject item
0.143 0.208 subject001 item01
Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, advice_present, predicted_lo, predicted_hi, predicted, subject, item, y_bin
Type: response 3.2.2 brms
The option re_formula = NA specifies that all random effects are set to 0 when computing predictions. The option re_formula = NULL specifies that all random effects are always included.
avg_predictions(fit_brms,
variables = list(advice_present = 0:1),
re_formula = NA, type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))
Term Estimate 2.5 % 97.5 %
b2-b1=0 0.177 0.144 0.21
Columns: term, estimate, conf.low, conf.high
Type: response avg_comparisons(fit_brms,
variables = "advice_present",
re_formula = NA, type = "response")
Term Contrast Estimate 2.5 % 97.5 %
advice_present mean(1) - mean(0) 0.177 0.144 0.21
Columns: term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx
Type: response predictions(fit_brms,
newdata = datagrid(advice_present = 0:1),
re_formula = NA, type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))
Term Estimate 2.5 % 97.5 %
b2-b1=0 0.177 0.144 0.21
Columns: term, estimate, conf.low, conf.high
Type: response comparisons(fit_brms,
newdata = datagrid(advice_present = 1),
re_formula = NA, type = "response")
Term Contrast advice_present Estimate 2.5 % 97.5 % subject item
advice_present 1 - 0 1 0.177 0.144 0.21 subject001 item01
Columns: rowid, term, contrast, estimate, conf.low, conf.high, advice_present, predicted_lo, predicted_hi, predicted, tmp_idx, subject, item, y_bin
Type: response 3.3 Treatment effect averaged across the actually observed persons and items
Estimand:
\[ \begin{aligned} \frac{1}{S \cdot I} \sum_{s} \sum_{i} & \quad P(Y = 1 | advice\_present = 1, u_{0s}, u_{1s}, u_{0i}) \\ & \quad - P(Y = 1 | advice\_present = 0, u_{0s}, u_{1s}, u_{0i}) \end{aligned} \]
3.3.1 lme4
avg_predictions(fit_lme4,
newdata = datagrid(advice_present = 0:1, subject = unique, item = unique),
by = "advice_present",
re.form = NULL, type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
b2-b1=0 0.167 0.0162 10.3 <0.001 80.1 0.135 0.199
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response avg_comparisons(fit_lme4,
variables = "advice_present",
re.form = NULL, type = "response")Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 %
advice_present mean(1) - mean(0) 0.167 0.0162 10.3 <0.001 80.1 0.135
97.5 %
0.199
Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
Type: response 3.3.2 brms
avg_predictions(fit_brms,
newdata = datagrid(advice_present = 0:1, subject = unique, item = unique),
by = "advice_present",
re_formula = NULL, type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))
Term Estimate 2.5 % 97.5 %
b2-b1=0 0.164 0.139 0.19
Columns: term, estimate, conf.low, conf.high
Type: response avg_comparisons(fit_brms,
variables = "advice_present",
re_formula = NULL, type = "response")
Term Contrast Estimate 2.5 % 97.5 %
advice_present mean(1) - mean(0) 0.164 0.139 0.19
Columns: term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx
Type: response 3.4 Treatment effect averaged across new persons and new items
Estimand:
\[ \begin{aligned} \frac{1}{S_{new} \cdot I_{new}} \sum_{s_{new}} \sum_{i_{new}} & \quad P(Y = 1 | advice\_present = 1, u_{0s_{new}}, u_{1s_{new}}, u_{0i_{new}}) \\ & \quad - P(Y = 1 | advice\_present = 0, u_{0s_{new}}, u_{1s_{new}}, u_{0i_{new}}) \end{aligned} \]
3.4.1 lme4
This cannot actually be done with lme4, which cannot sample new subjects or items (at least not with its predict function; it would work with the simulatefunction in lme4, which cannot be used by the marginaleffects package)! The code below produces the same results as the lme4 estimates for the Treatment effect for an average person and an average item!
avg_predictions(fit_lme4,
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = -1:-20),
by = "advice_present",
re.form = NULL, allow.new.levels = TRUE,
type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
b2-b1=0 0.176 0.0166 10.6 <0.001 84.8 0.143 0.208
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response avg_comparisons(fit_lme4,
variables = "advice_present",
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = -1:-20),
re.form = NULL, allow.new.levels = TRUE,
type = "response")Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 %
advice_present mean(1) - mean(0) 0.176 0.0166 10.6 <0.001 84.8 0.143
97.5 %
0.208
Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
Type: response 3.4.2 brms
The option sample_new_levels = "gaussian" specifies that for new factor levels, random effects are drawn from the estimated (multivariate) normal distribution of random effects. Note that this setting is not the default! Also note that allow_new_levels = TRUE will make the brms predictions for new levels non-deterministic!
The following two options produce slightly different results, and I currently do not know why! It might only be randomness introduced by different seeds.
set.seed(1)
avg_predictions(fit_brms,
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = -1:-20),
by = "advice_present",
re_formula = NULL, allow_new_levels = TRUE, sample_new_levels = "gaussian",
type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))
Term Estimate 2.5 % 97.5 %
b2-b1=0 0.162 0.121 0.203
Columns: term, estimate, conf.low, conf.high
Type: response set.seed(1)
avg_comparisons(fit_brms,
variables = "advice_present",
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = -1:-20),
re_formula = NULL, allow_new_levels = TRUE, sample_new_levels = "gaussian",
type = "response")
Term Contrast Estimate 2.5 % 97.5 %
advice_present mean(1) - mean(0) 0.163 0.122 0.203
Columns: term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx
Type: response 3.5 Treatment effect averaged across new persons but the actually observed items
Estimand:
\[ \begin{aligned} \frac{1}{S_{new} \cdot I} \sum_{s_{new}} \sum_{i} & \quad P(Y = 1 | advice\_present = 1, u_{0s_{new}}, u_{1s_{new}}, u_{0i}) \\ & \quad - P(Y = 1 | advice\_present = 0, u_{0s_{new}}, u_{1s_{new}}, u_{0i}) \end{aligned} \]
3.5.1 lme4
Because lme4 cannot sample new levels, what it actually does for a new subject is to set all random effects for the subject to 0.
avg_predictions(fit_lme4,
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = unique),
by = "advice_present",
re.form = NULL, allow.new.levels = TRUE,
type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
b2-b1=0 0.177 0.0164 10.8 <0.001 87.7 0.145 0.209
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response avg_comparisons(fit_lme4,
variables = "advice_present",
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = unique),
re.form = NULL, allow.new.levels = TRUE,
type = "response")Warning: For this model type, `marginaleffects` only takes into account the
uncertainty in fixed-effect parameters. You can use the `re.form=NA`
argument to acknowledge this explicitly and silence this warning.
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 %
advice_present mean(1) - mean(0) 0.177 0.0164 10.8 <0.001 87.7 0.145
97.5 %
0.209
Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
Type: response 3.5.2 brms
The following two options produce slightly different results, and I currently do not know why! It might only be randomness introduced by different seeds.
set.seed(1)
avg_predictions(fit_brms,
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = unique),
by = "advice_present",
re_formula = NULL, allow_new_levels = TRUE, sample_new_levels = "gaussian",
type = "response") %>%
hypotheses(hypothesis = c("b2 - b1 = 0"))
Term Estimate 2.5 % 97.5 %
b2-b1=0 0.163 0.126 0.199
Columns: term, estimate, conf.low, conf.high
Type: response set.seed(1)
avg_comparisons(fit_brms,
variables = "advice_present",
newdata = datagrid(advice_present = 0:1, subject = -1:-50, item = unique),
re_formula = NULL, allow_new_levels = TRUE, sample_new_levels = "gaussian",
type = "response")
Term Contrast Estimate 2.5 % 97.5 %
advice_present mean(1) - mean(0) 0.163 0.125 0.198
Columns: term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx
Type: response Reuse
Citation
@online{pargent2024,
author = {Pargent, Florian},
title = {Computing Predictions for Multilevel Models with the
Marginaleffects Package},
date = {2024-05-28},
url = {https://FlorianPargent.github.io//posts/002_marginaleffects/002_marginaleffects.html},
langid = {en}
}