使用
我正在尝试用emmeans::emtrends()函数估计连续变量的联合交互作用，但我在这样做时遇到了困难。如果有任何帮助，我们将不胜感激。请参阅下面的示例
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
#>
#> (Dispersion parameter for gaussian family taken to be 4.08059)
#>
#> Null deviance: 3217872 on 9999 degrees of freedom
#> Residual deviance: 40790 on 9996 degrees of freedom
#> AIC: 42447
#>
#> Number of Fisher Scoring iterations: 2
library("emmeans")
#mean difference in y for 1 unit increase in age_c across levels of bmi_c
(beta_age_c_bmi_c0 = 2)
#> [1] 2
(beta_age_c_bmi_c1 = 2 + 3.99)
#> [1] 5.99
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
#mean difference in y for 1 unit increase in bmi_c across levels of age_c
(beta_bmi_c_age_c0 = 3)
#> [1] 3
(beta_bmi_c_age_c1 = 3 + 3.99)
#> [1] 6.99
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
我想得到年龄_c增加1个单位和bmi_c增加1个单元的y的平均差。我想使用emtrends()获得它，但我不确定命令是什么。我正在尝试获得低于的输出
(beta_bmi_c1_age_c1 = 3 + 2 + 3.99)
#> [1] 8.99
创建于2022-10-05，reprex v2.0.2
我不知道你说的"联合交互"；，但从你问题的底线来看，你似乎只想要(1,1(和(0,0(的估计值之间的差异，其中坐标指的是(age_c, bmi_c)。获得这一点的一种方法是
emm <- emmeans(model_y, ~ age_c * bmi_c, at = list(age_c = 0:1, bmi_c = 0:1))
emm
contrast (emm, list(est = c(-1, 0, 0, 1)))
你的例子是不可复制的，因为你没有设定种子。所以我的结果有些不同：
> emm
age_c bmi_c emmean SE df lower.CL upper.CL
0 0 2.031 0.02004 9996 1.991 2.070
1 0 4.042 0.02241 9996 3.998 4.086
0 1 5.044 0.02236 9996 5.001 5.088
1 1 11.052 0.02502 9996 11.003 11.101
Confidence level used: 0.95
> contrast (emm, list(est = c(-1, 0, 0, 1)))
contrast estimate SE df t.ratio p.value
est 9.021 0.01501 9996 600.895 <.0001
使用emtrends()执行
使用emtrends()可以通过组合两个斜率来做同样的事情，并跟踪起点：
> (age_c.tr <- emtrends(model_y, "bmi_c", var = "age_c",
+ at = list(bmi_c = 0)))
bmi_c age_c.trend SE df lower.CL upper.CL
0 2.01 0.01 9996 1.99 2.03
Confidence level used: 0.95
> (bmi_c.tr <- emtrends(model_y, "age_c", var = "bmi_c",
+ at = list(age_c = 1)))
age_c bmi_c.trend SE df lower.CL upper.CL
1 7.01 0.01124 9996 6.988 7.032
Confidence level used: 0.95
> # Now sum these two results
> contrast(rbind(age_c.tr, bmi_c.tr), list(comb.trend = c(1, 1)))
contrast estimate SE df t.ratio p.value
comb.trend 9.021 0.01501 9996 600.895 <.0001
但我认为第一种方法更简单，也不太容易发生危险