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

发布于 2025-02-03 12:46:34 字数 1541 浏览 3 评论 0原文

我有一个具有三重互动的模型，类似于以下方式：

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

我试图根据am的效果如何根据cyl 和hp <am更改。 /代码>。使用效果（）函数从效果库中，我可以显示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))

有没有一种方法可以使用feform（）显示预测值之间的估计差距？

原文

I have a model with a triple interaction, similar to this:

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

I am trying to show how the effect of am changes based on the conditions of cyl and hp. Using the Effect() function from the effects library, I can show the predicted values of mpg. This works well and fairly quickly with my dataset. However, I want to show the size of the gap between the points in each case.

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)

I have tried using the margins() function from the margins library. As shown below. This shows the average marginal effect (AME), which I suppose is what I am trying to show. However, it takes an exorbitant amount of time with my dataset because I control for country fixed effects interacted with year and one of the independent variables.

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

Is there a way I can use Effect() to show the estimated gap between predicted values?

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

晚风撩人 2025-02-10 12:46:34

使您想要的绘图的所有相关信息都来自feform（）函数的输出。首先，让我们运行模型并生成效应对象。

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 < /code>：

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

^{在2022-06-01创建的}

All the relevant information to make the plot you want is in the output from the Effect() function. First, let's run the model and generate the effect object.

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)

Now, to get the difference between am == 0 and am == 1, we'll need to identify those values in the data frame p1.

w0 <- which(p1$am == 0)
w1 <- which(p1$am == 1)

We can then make a matrix that we will use to generate the differences in effects. We initialize it to have all zero values and it needs to have nrow(p1) rows and length(w0) columns:

D <- matrix(0, nrow = nrow(p1), ncol = length(w0))

Now, each column corresponds to a difference between the am == 0 and am == 1 prediction for a particular value of hp and cal. For the values where am == 0, we need the matrix to have -1 values and for am == 1, we need it to have values of 1. We can accomplish this as follows:

D[cbind(w0, 1:ncol(D))] <- -1
D[cbind(w1, 1:ncol(D))] <- 1

We can then get the differences this way:

diffs <- c(t(p1$fit) %*% D)

Just to make sure we got the right numbers, let's look at the first two values of 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

You can see that the first two values of diffs are the same as the ones we would calculate from p1. Now, we need to calculate the variance of the differences, we can do this as follows:

v_diffs <- t(D) %*% vcov(E) %*% D

Next, we make a dataset that will allow us to plot these differences. We keep the data where am == 0 just so we don't have replicated rows for each comparison. Then, we add in the differences, standard errors and confidence bounds.

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

Then, we can make the plot. Now, each point represents the difference between am==0 and am==1 for each different value of hp and 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())

^{Created on 2022-06-01 by the reprex package (v2.0.1)}

回复收藏 0 原文

凡尘雨 2025-02-10 12:46:34

这是使用 marginaleffects package 的替代方案。 “继任者”到Margins，具有更高的灵活性，更受支持的模型类型，并且通常更快地。（免责声明：我是作者。）

library(marginaleffects)

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

plot_slopes(m, variables = "am", condition = c("hp", "cyl"))

>自定义图

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

用ggplot2：

仅数字结果：

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

Here is an alternative using the marginaleffects package, which was designed as a “successor” to margins, with more flexibility, more supported model types, and is often much faster. (Disclaimer: I am the author.)

library(marginaleffects)

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

plot_slopes(m, variables = "am", condition = c("hp", "cyl"))

Customize the plot with 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()

Just the numeric results:

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

回复收藏 0 原文

~没有更多了~