如何用emmeans::emtrends()计算连续变量交互的联合趋势效应？

我正在尝试用emmeans::emtrends()函数估计连续变量的交互效应，目前已经能计算单个变量在另一变量不同水平下的趋势，但不知道如何通过该函数获取age_c和bmi_c同时增加1单位时y的均值变化。以下是我的模拟数据与代码示例：

library("tibble")
n=1e4
simdata <- tibble(
  age = rnorm(n, mean=30, sd=2),
  bmi = rnorm(n, mean=22, sd=2),
  age_c = age - mean(age),
  bmi_c = bmi - mean(bmi),
  y = rnorm(n, mean=2 + 2*age_c + 3*bmi_c + 4*age_c*bmi_c, sd=2)
)
model_y <- glm(y~age_c*bmi_c, family = gaussian(link = "identity"), data=simdata)
summary(model_y)
#> 
#> Call:
#> glm(formula = y ~ age_c * bmi_c, family = gaussian(link = "identity"), 
#>     data = simdata)
#> 
#> Deviance Residuals: 
#>    Min      1Q  Median      3Q     Max  
#> -7.773  -1.382   0.014   1.365   7.950  
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 1.990902   0.020204   98.54   <2e-16 ***
#> age_c       1.999852   0.010100  198.00   <2e-16 ***
#> bmi_c       3.000901   0.009994  300.27   <2e-16 ***
#> age_c:bmi_c 4.003130   0.005000  800.55   <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

我已经成功用emtrends()计算了单个变量的趋势：

library("emmeans")
# 不同bmi_c水平下，age_c每增加1单位时y的均值变化
emmeans::emtrends(model_y, 
         specs = "bmi_c",
         var = "age_c", 
         at = list(bmi_c= c(0, 1)))
#>  bmi_c age_c.trend      SE   df lower.CL upper.CL
#>      0       2.000 0.01010 9996    1.980    2.020
#>      1       6.003 0.01112 9996    5.981    6.025
#> 
#> Confidence level used: 0.95

# 不同age_c水平下，bmi_c每增加1单位时y的均值变化
emmeans::emtrends(model_y, 
         specs = "age_c",
         var = "bmi_c", 
         at = list(age_c= c(0, 1)))
#>  age_c bmi_c.trend       SE   df lower.CL upper.CL
#>      0       3.001 0.009994 9996    2.981    3.020
#>      1       7.004 0.011227 9996    6.982    7.026
#> 
#> Confidence level used: 0.95

现在我需要获取age_c和bmi_c同时增加1单位时y的均值变化，预期结果为2 + 3 + 4.003 ≈ 9.003（对应模型系数的组合）。

解决方案 #

emtrends()主要用于计算单个变量的趋势斜率，要获取两个连续变量同时变化的均值变化，更直接的方式是计算两个特定点的边际均值差值：

方法1：用emmeans() + contrast()计算差值 #

# 生成(0,0)和(1,1)两个点的边际均值
emm_points <- emmeans(model_y, ~ age_c * bmi_c, at = list(age_c = c(0,1), bmi_c = c(0,1)))

# 计算(1,1)与(0,0)的均值差
contrast(emm_points, list(simultaneous_change = c(-1, 0, 0, 1)))

运行后会得到包含均值变化量、标准误和置信区间的结果，这个差值就是两个变量同时增加1单位时y的变化，对应模型中age_c + bmi_c + age_c:bmi_c的系数和。

方法2：直接通过模型系数计算 #

如果只需要数值结果，也可以直接提取模型系数求和：

# 计算age_c、bmi_c及其交互项的系数和
coef(model_y)["age_c"] + coef(model_y)["bmi_c"] + coef(model_y)["age_c:bmi_c"]

这个结果和方法1完全一致，因为边际均值的差值本质上就是这三个系数的和。

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