我正在用 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)
我正在使用响应量表上的边际效应包计算治疗的平均边际效应：
# Marginal treatment effect
AME <- comparisons(
m1.WF_brms,
variables = "Treatment",
re_formula = NULL,
type = "response") %>%
tidy()
我的理解是，这个函数计算治疗效果，保持所有其他变量的平均效果（即治疗的人群水平效果）不变。这也适用于抵消吗？换句话说，我可以通过将 AME 除以平均暴露天数来获得每天治疗暴露的平均边际效应吗：
AME$estimate / mean(DF_anly_1$day_count)
我已经在这里实施了建议（见下文）：https://discourse.mc-stan.org/t/estimating-marginal-effects-of-offset-in-brms-poisson-regression-model/4463/4 ，但这似乎给出了条件边际效应而不是平均边际效应：
marginal_effects(m1.WF_brms,
conditions = data.frame(day_count = 1),
effects = "Treatment",
type = "response")
抱歉，我不明白您要找的具体数量。但我会用
marginaleffects首先，请注意，如果我们调用
brms::posterior_epred()library(brms)
library(marginaleffects)
p <- posterior_epred(m1.WF_brms, re_formula = NULL)
dim(p)
# [1] 8000 100
现在，对于原始数据中的每一行，我们在反事实世界中进行预测，其中
Treatmentd0 <- 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
这相当于调用
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
现在我们可以取这 100 个估计值的平均值。这相当于调用
avg_comparisons()avg_tidy()marginaleffectsmedian(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
最后，您还可以使用
byavg_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