I am trying to calculate the average marginal effect of a treatment (i.e., the average change in the response variable when participants are exposed to the treatment) per day of treatment exposure (which varies between participants).
I'm fitting the following Poisson regression with an offset in brms:
library(brms)
library(marginaleffects)
set.seed(42) # For reproducibility
# Number of observations
n <- 100
# Generate the data
DF_anly_1 <- data.frame(
catch_count = rpois(n, lambda = 2), # Simulated count data using a Poisson distribution
Treatment = sample(c('Control', 'Treatment'), n, replace = TRUE), # Treatment groups
village = sample(c('Village1', 'Village2', 'Village3'), n, replace = TRUE), # Village categories
phase_exp = sample(c('Phase1', 'Phase2'), n, replace = TRUE), # Experimental phases
boat = sample(c('Boat1', 'Boat2', 'Boat3'), n, replace = TRUE), # Random effects for boats
day_count = sample(1:30, n, replace = TRUE) # Days of count/exposure
)
# Model
m1.WF_brms <- brm(catch_count ~ Treatment + village + phase_exp + (1|boat) + offset(log(day_count)),
data = DF_anly_1, family = poisson(),
chains = 4, core = 4, iter = 4000)
I am calculating the average marginal effect of the treatment using the marginaleffects package on the response scale:
# Marginal treatment effect
AME <- comparisons(
m1.WF_brms,
variables = "Treatment",
re_formula = NULL,
type = "response") %>%
tidy()
My understanding is that this function calculates the treatment effect, holding constant the average effect of all other variables (i.e., the population-level effect of the treatment). Does this also apply to offsets? In other words, can I get the average marginal effect per day of treatment exposure by dividing the AME by the mean number of exposure days:
AME$estimate / mean(DF_anly_1$day_count)
I have implemented the advice here (see below): https://discourse.mc-stan.org/t/estimating-marginal-effects-of-offset-in-brms-poisson-regression-model/4463/4, but this appears to give the conditional rather than average marginal effect:
marginal_effects(m1.WF_brms,
conditions = data.frame(day_count = 1),
effects = "Treatment",
type = "response")
Sorry, but I don’t understand exactly what quantity you are looking for. But I’ll give you three options with marginaleffects and show manual computations so you know exactly what happens under the hood.
First, notice that if we call brms::posterior_epred() we get 8000 draws from the posterior distribution of predicted values (there are 100 points in your dataset):
library(brms)
library(marginaleffects)
p <- posterior_epred(m1.WF_brms, re_formula = NULL)
dim(p)
# [1] 8000 100
Now, for each row in the original data, we make predictions in the counterfactual worlds where Treatment is either “Control” or “Treatment”, and we take the difference between those predictions, and then take the median of the posterior draws. This gives us 100 estimates, one for each row:
d0 <- transform(DF_anly_1, Treatment = "Control")
d1 <- transform(DF_anly_1, Treatment = "Treatment")
p0 <- posterior_epred(m1.WF_brms, d0, re_formula = NULL)
p1 <- posterior_epred(m1.WF_brms, d1, re_formula = NULL)
effects <- apply(p1 - p0, 2, median)
head(effects)
# [1] -0.24480083 -0.14519547 -0.01207146 -0.09687572 -0.27200092 -0.29490911
This is equivalent to calling comparisons():
comparisons(
m1.WF_brms,
variables = "Treatment",
re_formula = NULL,
type = "response") |> head()
#
# Term Contrast Estimate 2.5 % 97.5 %
# Treatment Treatment - Control -0.2448 -1.081 0.4142
# Treatment Treatment - Control -0.1452 -0.575 0.2697
# Treatment Treatment - Control -0.0121 -0.054 0.0215
# Treatment Treatment - Control -0.0969 -0.389 0.1745
# Treatment Treatment - Control -0.2720 -1.201 0.4602
# Treatment Treatment - Control -0.2949 -1.293 0.5282
#
# Columns: rowid, term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx, catch_count, Treatment, village, phase_exp, day_count, boat
# Type: response
Now we can take the average of those 100 estimates. This is equivalent to calling avg_comparisons(). Warning: I use the avg_ prefix because the tidy() feature that your code relies on is now deprecated in marginaleffects 0.18.0:
median(apply(p1 - p0, 1, mean))
# [1] -0.2159907
avg_comparisons(
m1.WF_brms,
variables = "Treatment",
re_formula = NULL,
type = "response")
#
# Term Contrast Estimate 2.5 % 97.5 %
# Treatment Treatment - Control -0.216 -0.884 0.373
#
# Columns: term, contrast, estimate, conf.low, conf.high
# Type: response
Finally, you can also report aggregated estimates for each day of exposure using the by argument:
avg_comparisons(
m1.WF_brms,
by = "day_count",
variables = "Treatment",
re_formula = NULL,
type = "response")
#
# Term Contrast day_count Estimate 2.5 % 97.5 %
# Treatment mean(Treatment) - mean(Control) 1 -0.0143 -0.0582 0.0253
# Treatment mean(Treatment) - mean(Control) 2 -0.0253 -0.1075 0.0437
# Treatment mean(Treatment) - mean(Control) 3 -0.0388 -0.1676 0.0663
# Treatment mean(Treatment) - mean(Control) 4 -0.0536 -0.2318 0.0923
# Treatment mean(Treatment) - mean(Control) 5 -0.0702 -0.3078 0.1258
# Treatment mean(Treatment) - mean(Control) 6 -0.0822 -0.3423 0.1414
# Treatment mean(Treatment) - mean(Control) 7 -0.1010 -0.4138 0.1789
# Treatment mean(Treatment) - mean(Control) 8 -0.0966 -0.3994 0.1689
# Treatment mean(Treatment) - mean(Control) 9 -0.1457 -0.6396 0.2570
# Treatment mean(Treatment) - mean(Control) 10 -0.1409 -0.5689 0.2489
# Treatment mean(Treatment) - mean(Control) 11 -0.1428 -0.6072 0.2453
# Treatment mean(Treatment) - mean(Control) 12 -0.1599 -0.6642 0.2801
# Treatment mean(Treatment) - mean(Control) 13 -0.1775 -0.7256 0.3081
# Treatment mean(Treatment) - mean(Control) 14 -0.1822 -0.7467 0.3196
# Treatment mean(Treatment) - mean(Control) 15 -0.1954 -0.7978 0.3444
# Treatment mean(Treatment) - mean(Control) 16 -0.1968 -0.8214 0.3402
# Treatment mean(Treatment) - mean(Control) 17 -0.2415 -0.9836 0.4246
# Treatment mean(Treatment) - mean(Control) 18 -0.2448 -1.0810 0.4142
# Treatment mean(Treatment) - mean(Control) 19 -0.2328 -0.9672 0.3940
# Treatment mean(Treatment) - mean(Control) 20 -0.2533 -1.0799 0.4302
# Treatment mean(Treatment) - mean(Control) 21 -0.2893 -1.1689 0.5012
# Treatment mean(Treatment) - mean(Control) 22 -0.3194 -1.2644 0.5934
# Treatment mean(Treatment) - mean(Control) 23 -0.3380 -1.4055 0.5845
# Treatment mean(Treatment) - mean(Control) 24 -0.3073 -1.2686 0.5329
# Treatment mean(Treatment) - mean(Control) 25 -0.3229 -1.2969 0.5640
# Treatment mean(Treatment) - mean(Control) 26 -0.3215 -1.3020 0.5658
# Treatment mean(Treatment) - mean(Control) 27 -0.3411 -1.3967 0.5950
# Treatment mean(Treatment) - mean(Control) 28 -0.3539 -1.4158 0.6239
# Treatment mean(Treatment) - mean(Control) 29 -0.3764 -1.5477 0.6613
# Treatment mean(Treatment) - mean(Control) 30 -0.4241 -1.7160 0.7473
#
# Columns: term, contrast, day_count, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx
# Type: response