估计三重相互作用的边际效应

问题描述 投票：0回答：2

我有一个具有三重交互的模型，类似于：

m1 <- lm(mpg ~ am*cyl*hp, mtcars)

我试图展示

am
的效果如何根据
cyl
和
hp
的条件而变化。使用 
Effect()
库中的 
effects
函数，我可以显示 
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)

enter image description here

我尝试使用 

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))

enter image description here

有没有办法使用 

Effect()
显示预测值之间的估计差距？

r prediction effects interaction marginal-effects
2个回答
1
投票

制作所需绘图的所有相关信息都在 

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

然后，我们就可以制作情节了。 现在，每个点代表 

am==0
和 
am==1
的每个不同值的 
hp
和 
cyl
之间的差异：

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 包于 2022 年 6 月 1 日创建（v2.0.1）

1
投票

这里是使用 marginaleffects

 包
的替代方案，它被设计为 margins
 的“继承者”，具有更大的灵活性，更多支持的模型类型，并且通常更快
更快。 （免责声明：我是作者。） library(marginaleffects) m <- lm(mpg ~ am*cyl*hp, mtcars) plot_slopes(m, variables = "am", condition = c("hp", "cyl")) 


使用 
ggplot2
自定义绘图：
library(ggplot2)

plot_slopes(m, variables = "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


    
	

    

    


  

  

    
最新问题


  





  

    © www.soinside.com 2019 - 2024. All rights reserved.