我想测试 R 中不同回归的系数是否相同,但我不确定应该如何做。
让我举例说明。我有一个感兴趣的自变量,称为
treat
我想看看治疗对两个不同的结果 y1 和 y2 是否有相同的效果。
我该怎么办?
假设我有一个如下所示的数据集。
# Set seed for reproducibility
set.seed(123)
# Create example dataset
n <- 500 # Number of observations
groups <- 20 # Number of groups
years <- 10 # Number of years
# Simulate group and year identifiers
group <- rep(1:groups, each = n / groups)
yr <- rep(1:years, length.out = n)
# Simulate treatment assignment
treat <- sample(c(0, 1), size = n, replace = TRUE)
# Create covariates and treatment effects
X <- rnorm(n) # Random covariate
epsilon1 <- rnorm(n, sd = 2) # Error term for Y1
epsilon2 <- rnorm(n, sd = 2) # Error term for Y2
# Generate dependent variables
Y1 <- 5 + 2 * treat + 0.5 * X + epsilon1 # Y1: treatment effect = 2
Y2 <- 5 + 1 * treat + 0.5 * X + epsilon2 # Y2: treatment effect = 1
# Combine into a data frame
dt <- data.frame(group, yr, treat, X, Y1, Y2)
我的目标是评估
treat
相对于 Y1
对 Y2
的影响是否显着不同。
因此我估计了我的两个模型。在这种情况下,我使用两种方式的固定效果。
# Estimate the model for Y1
twfe_Y1 <- feols(Y1 ~ i(treat, ref = FALSE) | yr + group, data = dt, vcov = ~group)
# Estimate the model for Y2
twfe_Y2 <- feols(Y2 ~ i(treat, ref = FALSE) | yr + group, data = dt, vcov = ~group)
然后我可以使用 etable 比较模型
etable(twfe_Y1,twfe_Y2)
twfe_Y1 twfe_Y2
Dependent Var.: Y1 Y2
treat = 1 2.079*** (0.2368) 1.062*** (0.1361)
Fixed-Effects: ----------------- -----------------
yr Yes Yes
group Yes Yes
_______________ _________________ _________________
S.E.: Clustered by: group by: group
Observations 500 500
R2 0.22430 0.11058
Within R2 0.19325 0.05948
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
如何测试 Y1 的治疗系数是否与 Y2 的系数显着不同?
请注意,我无法堆叠 Y 并估计治疗组和 Y 组之间相互作用的模型,因为在实际分析中,我在估计过程中使用了一些新的差异包中的差异,并且尚不完全清楚它们如何能够包括与治疗的相互作用。
所以我需要一种事后获取系数和标准误差并测试系数是否相等的方法。
任何有关如何实施这一点的指导都会很棒。
正如 Vincent 所建议的,github.com/kylebutts/vcovSUR 应该能够满足您的需求。 这是如何用你的例子来做到这一点
# install.packages("devtools")
devtools::install_github("kylebutts/vcovSUR")
library(fixest)
library(vcovSUR)
sur_hypotheses(
list(twfe_Y1, twfe_Y2), cluster = "rowid",
hypothesis = "`treat::1_1` = `treat::1_2`" # mind the quotes!
)
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
1.02 0.261 3.9 <0.001 13.3 0.506 1.53
我也在此处阅读了一个非常相似的主题https://stats.stackexchange.com/questions/93540/testing-equality-of-coefficients-from-two- Different-regressions,建议使用 Z统计,我认为可以计算如下
beta1 <- as.numeric(twfe_Y1$coefficients)
se1<- as.numeric(twfe_Y1$se)
beta2 <- as.numeric(twfe_Y2$coefficients)
se2<- as.numeric(twfe_Y2$se)
Z <- (beta1 - beta2)/ sqrt(se1^2 + se2^2)
我没有深入研究它,但有趣的是,我注意到这两种方法返回非常相似(但不相同)的 p 值
pnorm(Z, lower.tail = F)
9.784923e-05
sur_hypotheses(
list(twfe_Y1, twfe_Y2), cluster = "rowid",
hypothesis = "`treat::1_1` = `treat::1_2`"
)$p.value
9.743827e-05