三阶多项式与二元处理交互项的部分置信区间求解

问题背景 #

在社会科学调查的意外事件设计/中断时间序列分析场景中，需要在自变量x的不同取值下，估算三阶多项式项（超出线性项的部分）的效应。目前已能得到点估计，但不确定如何计算置信区间，考虑使用marginaleffects包但不熟悉用法，也希望了解其他合适的R包、公式或代码。

核心关注两个关键效应：

• x = 0处的初始效应（即跳变）：post=1与post=0在x=0时的预测值差异
• 衰减效应（即超出效应发生前线性斜率的额外效应）：post=1组的斜率与post=0组斜率的差值

可复现模拟数据与模型 #

library(tidyverse)
library(marginaleffects) # plot_predictions属于marginaleffects包

set.seed(1900)

# 创建数据
N <- 200
quad <- data.frame(x = rnorm(N)) %>% 
  mutate(post = as.integer(x > 0),
         x_sq = x^2,
         x_th = x^3)

# 生成因变量：分段函数（小幅上升斜率+跳变+逐步回落至原斜率）
quad$y <- with(quad, 0.002 * x + 0.05 * post - 0.05 * x * post * (x <= 1) - 0.05 * (x > 1))

# 拟合三阶多项式模型
mod_quad <- lm(formula = y ~ x + x_sq + x_th + post + post * x + post * x_sq + post * x_th, data = quad)

# 可视化预测值、原始数据与基线轨迹
plot_predictions(mod_quad, by = "x", vcov = TRUE) + 
  geom_point(data = quad, aes(x = x, y = y), color = "blue", alpha = 0.5) +
  geom_line(data = tibble(x = seq(min(quad$x), max(quad$x), 0.01), y = x * 0.002),
            aes(x = x, y = y), linetype = "dashed", color = "darkorange") +
  theme_bw()

解决方案 #

方法1：使用marginaleffects包（推荐） #

marginaleffects是计算边际效应、对比值和置信区间的高效工具，支持线性模型的复杂效应计算。

1.1 计算x=0处的初始跳变效应

直接对比post=1和post=0在x=0时的预测值，自动生成置信区间：

# 计算跳变效应及置信区间
jump_effect <- comparisons(
  mod_quad,
  variables = list(post = c(0, 1)),
  newdata = data.frame(x = 0, x_sq = 0, x_th = 0)
)

# 查看结果（包含点估计、标准误、置信区间）
print(jump_effect)

1.2 计算不同x值下的衰减效应（额外斜率）

衰减效应是两组斜率的差值，可通过emtrends（marginaleffects兼容）直接计算：

# 生成x序列
x_seq <- seq(min(quad$x), max(quad$x), length.out = 100)

# 用emtrends计算两组斜率的对比
attenuation_effect <- emtrends(
  mod_quad,
  pairwise ~ post,
  var = "x",
  at = list(x = x_seq, x_sq = x_seq^2, x_th = x_seq^3)
)$contrasts

# 转换为tibble方便可视化
attenuation_df <- as_tibble(attenuation_effect) %>%
  rename(x = x, effect = estimate, se = SE, lower_ci = lower.CL, upper_ci = upper.CL)

# 可视化衰减效应及置信区间
ggplot(attenuation_df, aes(x = x, y = effect)) +
  geom_line(color = "black") +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2) +
  labs(title = "衰减效应（额外斜率）", x = "运行变量x", y = "post=1与post=0的斜率差") +
  theme_bw()

方法2：手动计算置信区间（基于方差-协方差矩阵） #

对于线性模型，置信区间可通过回归系数的线性组合结合t分布计算，适合理解底层逻辑。

2.1 初始跳变效应的置信区间

当x=0时，跳变效应等于post的系数（所有交互项为0）：

# 提取post系数及标准误
beta_post <- coef(mod_quad)["post"]
se_post <- sqrt(vcov(mod_quad)["post", "post"])
df_resid <- mod_quad$df.residual

# 计算95%置信区间
ci_jump <- c(
  beta_post - qt(0.975, df = df_resid) * se_post,
  beta_post + qt(0.975, df = df_resid) * se_post
)
names(ci_jump) <- c("lower_ci", "upper_ci")
print(ci_jump)

2.2 衰减效应的置信区间

衰减效应的表达式为：Δ斜率 = β5 + 2β6x + 3β7x²（β5是post:x系数，β6是post:x_sq系数，β7是post:x_th系数）：

# 示例：计算x=0.5处的衰减效应置信区间
x_val <- 0.5
# 线性组合权重
weights <- c(0, 0, 0, 0, 1, 2*x_val, 3*x_val^2)

# 点估计
atten_est <- sum(coef(mod_quad) * weights)
# 标准误
atten_se <- sqrt(t(weights) %*% vcov(mod_quad) %*% weights)
# 95%置信区间
ci_atten <- c(
  atten_est - qt(0.975, df = df_resid) * atten_se,
  atten_est + qt(0.975, df = df_resid) * atten_se
)
names(ci_atten) <- c("lower_ci", "upper_ci")
print(ci_atten)

方法3：使用emmeans包 #

emmeans专注于边际均值和趋势的估计，适合复杂模型的效应对比：

install.packages("emmeans")
library(emmeans)

# 计算x=0处的跳变效应
emmeans(mod_quad, pairwise ~ post, at = list(x=0, x_sq=0, x_th=0))

# 计算不同x下的衰减效应
emm_trends <- emtrends(mod_quad, pairwise ~ post, var = "x", at = list(x = x_seq))
attenuation_emmeans <- as_tibble(emm_trends$contrasts)

内容的提问来源于stack exchange，提问作者Ryan Baxter-King