分位数回归中emmeans与手动计算多交互项差异不符问题

问题背景 #

在R v4.2.1中执行分位数回归（tau=0.5）时，使用emmeans v1.8.1-1提取成对差异可得到正确结果，但手动计算其他分位数的差异时，结果与emmeans输出不一致。具体场景：

• 变量：var1（二分类：A/B）、var2（二分类：High/Low）、var3_z（标准化，均值0、标准差1）
• 模型包含var1×var2、var1×var3_z交互项
• 对比发现：var2=Low时，手动计算var1A与B的差异为1.36，但emmeans显示为1.3（排除四舍五入问题）

模型与结果展示 #

分位数回归模型摘要 #

modelAll50 <- rq(output ~ var1 * var2 + var1 * var3_z, tau = 0.5, data = dfModelAllControl, method = "fn")
summary(modelAll50)

输出：

Call: rq(formula = output ~ var1 * var2 + var1 * var3_z, tau = 0.5, 
    data = dfModelAllControl, method = "fn")

tau: [1] 0.5

Coefficients:
               Value    Std. Error t value  Pr(>|t|)
(Intercept)     0.04322  0.01623    2.66359  0.00774
var1B           1.36359  0.19793    6.88936  0.00000
var2High        0.11678  0.04986    2.34223  0.01919
var3_z         -0.02829  0.01237   -2.28627  0.02226
var1B:var2High  6.60083  0.65356   10.09977  0.00000
var1B:var3_z   -0.18197  0.21099   -0.86245  0.38846

emmeans成对差异结果 #

em <- emmeans(modelAll50, pairwise ~ var1 | var2)
pairs(em) %>%  confint()

输出：

var2 = Low:
 contrast estimate    SE    df lower.CL upper.CL
 A - B        -1.3 0.207 10023    -1.70   -0.895

var2 = High:
 contrast estimate    SE    df lower.CL upper.CL
 A - B        -7.9 0.626 10023    -9.13   -6.673

Results are averaged over the levels of: var3_z 
Confidence level used: 0.95 

核心原因 #

你手动计算的是固定var3_z=0时的条件分位数差异，而emmeans输出的是对var3_z的全部分布进行边际平均后的分位数差异：

1. 当var1与var3_z存在交互项时，var1的效应会随var3_z的取值变化
2. emmeans默认对协变量（连续变量var3_z）的分布进行边际平均，而非固定在均值处
3. 即使var3_z已标准化，只要其分布不是单点0，交互项的存在就会让边际平均结果与固定条件结果产生偏差

手动修正计算方法 #

要得到与emmeans一致的结果，需按以下步骤计算：

1. 计算边际平均差异估计值 #

对数据中每个观测的var3_z值，计算var1=B与var1=A的分位数预测差异，再取均值：

# 提取模型系数
coefs <- coef(modelAll50)
# 计算每个观测的B-A差异（var2=Low时）
dfModelAllControl$diff_B_A <- coefs["var1B"] + coefs["var1B:var3_z"] * dfModelAllControl$var3_z
# 边际平均差异（对应emmeans中B-A的估计值，A-B为其相反数）
mean_diff <- mean(dfModelAllControl$diff_B_A)

2. 计算标准误 #

利用模型的协方差矩阵，结合var3_z的分布特征计算：

# 获取系数协方差矩阵
vcov_mat <- vcov(modelAll50)
# 由于var3_z标准化，mean(var3_z)=0、mean(var3_z²)=1，简化计算边际差异的方差
var_diff <- vcov_mat["var1B","var1B"] + vcov_mat["var1B:var3_z","var1B:var3_z"] + 2*vcov_mat["var1B","var1B:var3_z"]
# 标准误
se_diff <- sqrt(var_diff)

额外说明 #

当var1与var3_z无交互时，var1B:var3_z系数为0，此时边际平均差异等于固定var3_z=0的条件差异，因此手动计算与emmeans结果一致，这与你观察到的现象匹配。

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