lm_robust模型A与V交互项边际效应与Stata结果不一致求助

核心差异原因 #

你在R和Stata中得到不同结果的核心原因是边际效应的计算逻辑与对象不一致：

• Stata的margins A#V输出的是每个A-V组合的调整后平均预测值（即固定A和V的水平，其他变量取均值/模式时的平均预测结果）；
• 你在R中用subset(df, V == 0)计算的是V=0子集中A的平均边际效应（即A各水平相对于基准的预测值差异在V=0样本中的平均值）。

此外，两者对其他控制变量的调整方式不同：Stata默认将控制变量调整到均值/模式，而你R中用子集计算时保留了原始样本的控制变量取值，这也会导致结果偏差。

步骤1：确认模型拟合一致性 #

首先需确保R和Stata的模型系数完全一致，这是后续计算的基础：

1. 在R中统一因子变量的基准水平（与Stata默认的最小水平对齐）：

library(estimatr)
library(margins)

# 将所有分类变量转为因子，指定水平顺序（与Stata一致）
df$A <- factor(df$A, levels = 1:5)
df$V <- factor(df$V, levels = 0:1)
df$Y <- factor(df$Y)
df$W <- factor(df$W)
df$T <- factor(df$T)

# 重新拟合模型
model1 <- lm_robust(P ~ (A + Y + W + T)*V,
                    se_type = "stata", cluster = ID, data = df)

# 查看系数
summary(model1)$coefficients

2. 在Stata中运行回归后输出系数，对比两者的主效应、交互项系数是否完全一致。若系数不一致，需检查数据编码、缺失值处理是否统一。

步骤2：在R中复现Stata margins A#V结果 #

若要得到与Stata完全一致的A-V组合平均预测值，需模拟Stata的调整逻辑（控制变量取均值/模式，计算所有A-V组合的平均预测值）：

# 定义获取众数的函数（分类变量用众数调整）
get_mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

# 生成控制变量的调整值（分类变量取众数）
adj_params <- list(
  Y = get_mode(df$Y),
  W = get_mode(df$W),
  T = get_mode(df$T)
)

# 生成所有A-V组合的网格
at_grid <- expand.grid(
  A = levels(df$A),
  V = levels(df$V),
  adj_params
)

# 计算调整后平均预测值
margins_result <- margins(model1, at = at_grid)
margins_summary <- summary(margins_result)

# 整理结果格式，匹配Stata输出
final_result <- margins_summary[, c("factor", "AME", "SE", "z", "p", "lower", "upper")]
colnames(final_result) <- c("A#V", "Margin", "Std. Err.", "z", "P>|z|", "[95% Conf. Interval Lower]", "[95% Conf. Interval Upper]")

print(final_result)

步骤3：计算A在不同V下的边际效应（若为真实需求） #

如果你实际需要的是A在V=0和V=1下的边际效应（即A各水平相对于基准的效应差异），对应Stata的margins, dydx(A) at(V=(0 1))，可使用以下代码：

# 计算V=0时A的边际效应
margins_v0 <- margins(model1, at = list(V = "0"), variables = "A")
# 计算V=1时A的边际效应
margins_v1 <- margins(model1, at = list(V = "1"), variables = "A")

# 合并结果
effect_result <- rbind(
  cbind(summary(margins_v0), V = "0"),
  cbind(summary(margins_v1), V = "1")
)

# 重命名列名
colnames(effect_result)[colnames(effect_result) == "AME"] <- "Marginal Effect"
print(effect_result)

内容的提问来源于stack exchange，提问作者Saidbf