估计三重交互的边际效应
Estimate marginal effect in triple interaction
我有一个具有三重交互的模型，类似于：
m1 <- lm(mpg ~ am*cyl*hp, mtcars)
我试图展示 am 的效果如何根据 cyl 和 hp 的条件而变化。使用 effects 库中的 Effect() 函数，我可以显示 mpg 的预测值。这对我的数据集运行良好且速度相当快。但是，我想显示每种情况下点之间的差距大小。
library(effects)
p1 <- as.data.frame(Effect(m1, focal.predictors = c("am", "cyl", "hp"), xlevels = list(am=c(0, 1), cyl=c(4,8), hp = c(100, 200))))
library(ggplot2)
ggplot(p1, aes(cyl, fit, color = as.factor(am))) +
geom_point(position = position_dodge(0.5)) +
geom_errorbar(aes(ymin=lower, ymax=upper), width = 0, position = position_dodge(0.5)) +
facet_grid(~hp)
我尝试使用 margins 库中的 margins() 函数。如下所示。这显示了平均边际效应 (AME)，我想这就是我想要展示的。但是，我的数据集花费了大量时间，因为我控制了与年份和一个自变量相互作用的国家/地区固定效应。
p2 <- margins(m1, at=list(cyl = c(4, 8), hp = c(100, 200)), variables = "am")
p2 <- summary(p2)
ggplot(p2, aes(cyl, AME, color = as.factor(hp))) +
geom_point(position = position_dodge(0.5)) +
geom_errorbar(aes(ymin=lower, ymax=upper), width = 0, position = position_dodge(0.5))
有什么方法可以使用 Effect() 来显示预测值之间的估计差距？
制作所需情节的所有相关信息都在 Effect() 函数的输出中。首先，我们 运行 模型并生成效果对象。
library(effects)
#> Loading required package: carData
#> lattice theme set by effectsTheme()
#> See ?effectsTheme for details.
library(ggplot2)
data(mtcars)
m1 <- lm(mpg ~ am*cyl*hp + wt + drat, mtcars)
E <- Effect(m1,
focal.predictors = c("am", "cyl", "hp"),
xlevels = list(am=c(0, 1),
cyl=c(4,8),
hp = seq(52, 335, length=10)))
p1 <- as.data.frame(E)
现在，要获得 am == 0 和 am == 1 之间的差异，我们需要在数据框中识别这些值 p1。
w0 <- which(p1$am == 0)
w1 <- which(p1$am == 1)
然后我们可以制作一个矩阵，我们将使用它来产生效果的差异。我们将它初始化为全零值，它需要有 nrow(p1) 行和 length(w0) 列：
D <- matrix(0, nrow = nrow(p1), ncol = length(w0))
现在，每一列对应于 am == 0 和 am == 1 预测值之间的差异，具体值为 hp 和 cal。对于 am == 0 的值，我们需要矩阵具有 -1 值，对于 am == 1，我们需要它具有值 1。我们可以按如下方式完成此操作：
D[cbind(w0, 1:ncol(D))] <- -1
D[cbind(w1, 1:ncol(D))] <- 1
然后我们可以这样得到差异：
diffs <- c(t(p1$fit) %*% D)
为了确保我们得到正确的数字，让我们看一下 diffs 的前两个值：
diffs[1:2]
#> [1] 3.2716241 -0.8526864
p1[1:4, ]
#> am cyl hp fit se lower upper
#> 1 0 4 52 24.74134 2.784239 18.967181 30.51550
#> 2 1 4 52 28.01296 2.203560 23.443061 32.58287
#> 3 0 8 52 19.12017 3.466758 11.930556 26.30979
#> 4 1 8 52 18.26749 4.455793 9.026738 27.50823
p1$fit[2]-p1$fit[1]
#> [1] 3.271624
p1$fit[4]-p1$fit[3]
#> [1] -0.8526864
您可以看到 diffs 的前两个值与我们从 p1 计算的值相同。现在，我们需要计算差异的方差，我们可以这样做：
v_diffs <- t(D) %*% vcov(E) %*% D
接下来，我们制作一个数据集，让我们能够绘制出这些差异。我们将数据保存在 am == 0 的位置，这样我们就不会为每次比较复制行。然后，我们添加差异、标准误差和置信区间。
p11 <- p1[p1$am == 0, ]
p11$diff <- diffs
p11$se_diff <- sqrt(diag(v_diffs))
p11$lwr <- p11$diff - 1.96*p11$se_diff
p11$upr <- p11$diff + 1.96*p11$se_diff
然后，我们就可以制作剧情了。现在，对于 hp 和 cyl 的每个不同值，每个点代表 am==0 和 am==1 之间的差异：
ggplot(p11, aes(x=hp, y=diff, ymin=lwr, ymax=upr)) +
geom_pointrange() +
geom_hline(yintercept=0, linetype=2, col="red") +
facet_wrap(~as.factor(cyl)) +
theme_bw() +
theme(panel.grid=element_blank())
由 reprex package (v2.0.1)
创建于 2022-06-01
这是一个使用 the marginaleffects package 的替代方案，它被设计为 margins 的“继承者”，具有更大的灵活性，支持更多的模型类型，并且通常 多 更快。 （免责声明：我是作者。）
library(marginaleffects)
m <- lm(mpg ~ am*cyl*hp, mtcars)
plot_cme(m, effect = "am", condition = c("hp", "cyl"))
自定义绘图 ggplot2:
library(ggplot2)
plot_cme(m, effect = "am", condition = c("hp", "cyl"), draw = FALSE) |>
ggplot(aes(condition1, dydx, ymin = conf.low, ymax = conf.high)) +
geom_ribbon(alpha = .2) +
geom_line() +
facet_wrap(~condition2) +
theme_classic()
只是数字结果：
mfx <- marginaleffects(m)
summary(mfx, by = "cyl")
#> Average marginal effects
#> Term cyl Effect Std. Error z value Pr(>|z|) 2.5 % 97.5 %
#> 1 am 6 2.00971 1.95284 1.0291 0.3034237 -1.81779 5.837211
#> 2 am 4 5.08709 1.76482 2.8825 0.0039453 1.62811 8.546073
#> 3 am 8 2.15232 2.23291 0.9639 0.3350919 -2.22410 6.528733
#> 4 cyl 6 -0.94326 0.71546 -1.3184 0.1873742 -2.34554 0.459026
#> 5 cyl 4 -2.06503 0.85842 -2.4056 0.0161455 -3.74751 -0.382553
#> 6 cyl 8 0.47177 1.71136 0.2757 0.7828002 -2.88243 3.825974
#> 7 hp 6 -0.05691 0.02592 -2.1956 0.0281202 -0.10772 -0.006108
#> 8 hp 4 -0.09174 0.03417 -2.6852 0.0072497 -0.15870 -0.024776
#> 9 hp 8 -0.03583 0.01893 -1.8930 0.0583573 -0.07293 0.001267
#>
#> Model type: lm
#> Prediction type: response
我有一个具有三重交互的模型，类似于：
m1 <- lm(mpg ~ am*cyl*hp, mtcars)
我试图展示 am 的效果如何根据 cyl 和 hp 的条件而变化。使用 effects 库中的 Effect() 函数，我可以显示 mpg 的预测值。这对我的数据集运行良好且速度相当快。但是，我想显示每种情况下点之间的差距大小。
library(effects)
p1 <- as.data.frame(Effect(m1, focal.predictors = c("am", "cyl", "hp"), xlevels = list(am=c(0, 1), cyl=c(4,8), hp = c(100, 200))))
library(ggplot2)
ggplot(p1, aes(cyl, fit, color = as.factor(am))) +
geom_point(position = position_dodge(0.5)) +
geom_errorbar(aes(ymin=lower, ymax=upper), width = 0, position = position_dodge(0.5)) +
facet_grid(~hp)
我尝试使用 margins 库中的 margins() 函数。如下所示。这显示了平均边际效应 (AME)，我想这就是我想要展示的。但是，我的数据集花费了大量时间，因为我控制了与年份和一个自变量相互作用的国家/地区固定效应。
p2 <- margins(m1, at=list(cyl = c(4, 8), hp = c(100, 200)), variables = "am")
p2 <- summary(p2)
ggplot(p2, aes(cyl, AME, color = as.factor(hp))) +
geom_point(position = position_dodge(0.5)) +
geom_errorbar(aes(ymin=lower, ymax=upper), width = 0, position = position_dodge(0.5))
有什么方法可以使用 Effect() 来显示预测值之间的估计差距？
制作所需情节的所有相关信息都在 Effect() 函数的输出中。首先，我们 运行 模型并生成效果对象。
library(effects)
#> Loading required package: carData
#> lattice theme set by effectsTheme()
#> See ?effectsTheme for details.
library(ggplot2)
data(mtcars)
m1 <- lm(mpg ~ am*cyl*hp + wt + drat, mtcars)
E <- Effect(m1,
focal.predictors = c("am", "cyl", "hp"),
xlevels = list(am=c(0, 1),
cyl=c(4,8),
hp = seq(52, 335, length=10)))
p1 <- as.data.frame(E)
现在，要获得 am == 0 和 am == 1 之间的差异，我们需要在数据框中识别这些值 p1。
w0 <- which(p1$am == 0)
w1 <- which(p1$am == 1)
然后我们可以制作一个矩阵，我们将使用它来产生效果的差异。我们将它初始化为全零值，它需要有 nrow(p1) 行和 length(w0) 列：
D <- matrix(0, nrow = nrow(p1), ncol = length(w0))
现在，每一列对应于 am == 0 和 am == 1 预测值之间的差异，具体值为 hp 和 cal。对于 am == 0 的值，我们需要矩阵具有 -1 值，对于 am == 1，我们需要它具有值 1。我们可以按如下方式完成此操作：
D[cbind(w0, 1:ncol(D))] <- -1
D[cbind(w1, 1:ncol(D))] <- 1
然后我们可以这样得到差异：
diffs <- c(t(p1$fit) %*% D)
为了确保我们得到正确的数字，让我们看一下 diffs 的前两个值：
diffs[1:2]
#> [1] 3.2716241 -0.8526864
p1[1:4, ]
#> am cyl hp fit se lower upper
#> 1 0 4 52 24.74134 2.784239 18.967181 30.51550
#> 2 1 4 52 28.01296 2.203560 23.443061 32.58287
#> 3 0 8 52 19.12017 3.466758 11.930556 26.30979
#> 4 1 8 52 18.26749 4.455793 9.026738 27.50823
p1$fit[2]-p1$fit[1]
#> [1] 3.271624
p1$fit[4]-p1$fit[3]
#> [1] -0.8526864
您可以看到 diffs 的前两个值与我们从 p1 计算的值相同。现在，我们需要计算差异的方差，我们可以这样做：
v_diffs <- t(D) %*% vcov(E) %*% D
接下来，我们制作一个数据集，让我们能够绘制出这些差异。我们将数据保存在 am == 0 的位置，这样我们就不会为每次比较复制行。然后，我们添加差异、标准误差和置信区间。
p11 <- p1[p1$am == 0, ]
p11$diff <- diffs
p11$se_diff <- sqrt(diag(v_diffs))
p11$lwr <- p11$diff - 1.96*p11$se_diff
p11$upr <- p11$diff + 1.96*p11$se_diff
然后，我们就可以制作剧情了。现在，对于 hp 和 cyl 的每个不同值，每个点代表 am==0 和 am==1 之间的差异：
ggplot(p11, aes(x=hp, y=diff, ymin=lwr, ymax=upr)) +
geom_pointrange() +
geom_hline(yintercept=0, linetype=2, col="red") +
facet_wrap(~as.factor(cyl)) +
theme_bw() +
theme(panel.grid=element_blank())
由 reprex package (v2.0.1)
创建于 2022-06-01这是一个使用 the marginaleffects package 的替代方案，它被设计为 margins 的“继承者”，具有更大的灵活性，支持更多的模型类型，并且通常 多 更快。 （免责声明：我是作者。）
library(marginaleffects)
m <- lm(mpg ~ am*cyl*hp, mtcars)
plot_cme(m, effect = "am", condition = c("hp", "cyl"))
自定义绘图 ggplot2:
library(ggplot2)
plot_cme(m, effect = "am", condition = c("hp", "cyl"), draw = FALSE) |>
ggplot(aes(condition1, dydx, ymin = conf.low, ymax = conf.high)) +
geom_ribbon(alpha = .2) +
geom_line() +
facet_wrap(~condition2) +
theme_classic()
只是数字结果：
mfx <- marginaleffects(m)
summary(mfx, by = "cyl")
#> Average marginal effects
#> Term cyl Effect Std. Error z value Pr(>|z|) 2.5 % 97.5 %
#> 1 am 6 2.00971 1.95284 1.0291 0.3034237 -1.81779 5.837211
#> 2 am 4 5.08709 1.76482 2.8825 0.0039453 1.62811 8.546073
#> 3 am 8 2.15232 2.23291 0.9639 0.3350919 -2.22410 6.528733
#> 4 cyl 6 -0.94326 0.71546 -1.3184 0.1873742 -2.34554 0.459026
#> 5 cyl 4 -2.06503 0.85842 -2.4056 0.0161455 -3.74751 -0.382553
#> 6 cyl 8 0.47177 1.71136 0.2757 0.7828002 -2.88243 3.825974
#> 7 hp 6 -0.05691 0.02592 -2.1956 0.0281202 -0.10772 -0.006108
#> 8 hp 4 -0.09174 0.03417 -2.6852 0.0072497 -0.15870 -0.024776
#> 9 hp 8 -0.03583 0.01893 -1.8930 0.0583573 -0.07293 0.001267
#>
#> Model type: lm
#> Prediction type: response