与二元处理交互的多项式中某些（但不是全部）交互项的置信区间

Question

我正在寻求帮助来估计三阶多项式项的效果——在线性项之上在自变量

x

的不同值。虽然我可以估计点估计值，但我不确定如何计算置信区间。我认为边际效应和

marginaleffects

包可能是这里的答案——但这对我来说相对较新，我不知道如何将它们放在一起。

如果您有不同的包、公式或代码，看起来像是正确的答案，我也绝对对此感兴趣！这里的一般应用是使用社会科学调查数据的意外事件设计/中断时间序列。

这是一个可重现的示例。在此示例中，有一个运行变量

X

，范围介于

x = -2.5

和

x = 2.5

之间，以及结果变量

Y

。有两个感兴趣的量：

1. x = 0

   处的初始效果（即跳跃）
2. 衰减效应（即效应之前发生的线性斜率之上的额外效应

在此示例中，结果变量

Y

是

X

的逐步函数。我使用三阶多项式来近似这个函数，作为模拟的设计者我们知道，即使它是不正确的。

这是一个图（由下面的代码生成），它可视化了模拟数据（蓝点）、原始轨迹（橙色虚线）和三阶多项式（黑线）。

...这是代码

library(margins)
library(tidyverse)


set.seed(1900)

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

# Create outcome (stepwise outcome with slight upward slope, jump, gradual decrease back to original slope)
quad$y <- with(quad, 0.002 * x + 0.05 * post - 0.05 * x * post * (x <= 1) - 0.05 * (x > 1))

# Visualize outcome
plot(quad$x, quad$y)

# Specify model (third order polynomial)

mod_quad <- lm(formula = y ~ x + x_sq + x_th + post + post * x + post * x_sq + post * x_th, data = quad)

# Plot Predictions

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()

感谢您的帮助！如果有任何需要澄清的问题请告诉我。

Answer 1

您不应该对多项式项进行硬编码。为了计算标准误差，

marginaleffects

需要通过稍微调整

x

的值来获得数值导数。但是，当您对多项式项进行硬编码时，

marginaleffects

会调整

x

，但只有

x

列受到影响，而不是

x_sq

、

x_th

等。

相反，请在模型拟合调用中使用

poly()

函数，如评论中所建议：

library(dplyr)
library(ggplot2)
library(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 <- lm(formula = y ~ post * poly(x, 3), data = quad)

plot_predictions(mod, condition = "x")  +
  geom_point(data = quad, aes(x, y))