5

我无法在 R 中复制 Statamargins命令的特定用例： margins var1, over(var2) 我一直在尝试使用marginsR 中的包来这样做。

为了提供一个可重现的示例，我使用了 mtcars 数据集并将其从 R 导出到 Stata，因此我们在两个程序中使用相同的数据集：

代码：

library(foreign)
library(margins)
write.dta(mtcars, “mtcars.dta")

状态码：

use "mtcars.dta", clear

在两个程序中创建示例线性回归模型

状态码：

quietly regress mpg cyl i.am c.wt##c.hp

代码：

x <- lm(mpg ~ cyl + factor(am) + hp * wt, data = mtcars)

两个程序之间的模型输出（未显示）相同

比较模型中每个变量的平均边际效应表

状态代码和输出：

margins, dydx(*)

Average marginal effects                          Number of obs   =         32
Model VCE: OLS

Expression   : Linear prediction, predict() dy/dx w.r.t. : cyl 1.am wt hp

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |  -.3708001   .5293674    -0.70   0.490     -1.45893    .7173301
        1.am |  -.0709546   1.374981    -0.05   0.959    -2.897268    2.755359
          wt |  -3.868994   .9170145    -4.22   0.000    -5.753944   -1.984043
          hp |  -.0249882   .0120345    -2.08   0.048    -.0497254    -.000251
------------------------------------------------------------------------------ 
Note: dy/dx for factor levels is the discrete change from the base level.

R代码和输出：

xmarg <- margins(x)
summary(xmarg)

factor     AME     SE       z      p   lower   upper
    am1 -0.0710 1.3750 -0.0516 0.9588 -2.7659  2.6240
    cyl -0.3708 0.5294 -0.7005 0.4836 -1.4083  0.6667
     hp -0.0250 0.0120 -2.0764 0.0379 -0.0486 -0.0014
     wt -3.8690 0.9170 -4.2191 0.0000 -5.6663 -2.0717

如您所见，这两个输出非常相似，正如使用 Rmargins包所期望的那样。

问题 1：对变量值的边际预测

状态代码和输出：

margins, over(cyl)

Predictive margins                                Number of obs   =         32
Model VCE: OLS

Expression   : Linear prediction, predict()
over         : cyl

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         cyl |
          4  |   26.56699   .6390379    41.57   0.000     25.25342    27.88055
          6  |   20.04662   .5797511    34.58   0.000     18.85492    21.23831
          8  |   15.02406   .5718886    26.27   0.000     13.84853    16.19959
------------------------------------------------------------------------------

R代码和输出：

aggregate(fitted~cyl, data = xmarg, FUN = mean)
  cyl   fitted
1   4 26.56699
2   6 20.04662
3   8 15.02406

在上面的两个示例中，R 和 Stata 之间的边际预测是相同的。但是，有没有一种方法（没有手动完成）为每个边际预测生成 delta 方法标准误差，就像上面的 Stata 表中所做的那样？

问题 2：特定变量的边际预测：

状态代码和输出：

margins am

Predictive margins                                Number of obs   =         32
Model VCE    : OLS

Expression   : Linear prediction, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          am |
          0  |   20.11945   .6819407    29.50   0.000      18.7177     21.5212
          1  |    20.0485   .9052764    22.15   0.000     18.18767    21.90932
------------------------------------------------------------------------------

R代码和输出：

aggregate(fitted~am, data = xmarg, FUN = mean)
  am   fitted
1  0 17.14737
2  1 24.39231

在这个例子中，我们试图margins通过在预测后对数据集进行子集来复制命令中的 Stata 的“marginlist”参数。这似乎不是正确的方法。我们如何从 Stata 复制这些结果？

问题 3：一个变量对另一个变量值的边际预测

复制这个结果是我的主要目标！

状态代码和输出

margins am, over(cyl)

Predictive margins                                Number of obs   =         32
Model VCE    : OLS

Expression   : Linear prediction, predict()
over         : cyl

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      cyl#am |
        4 0  |   26.61859   1.246074    21.36   0.000     24.05725    29.17993
        4 1  |   26.54763   .7034599    37.74   0.000     25.10165    27.99362
        6 0  |   20.07703   .6449805    31.13   0.000     18.75125     21.4028
        6 1  |   20.00607   1.144518    17.48   0.000     17.65348    22.35866
        8 0  |    15.0342   .6228319    24.14   0.000     13.75395    16.31445
        8 1  |   14.96324   1.257922    11.90   0.000     12.37754    17.54894
------------------------------------------------------------------------------

R代码和输出：

aggregate(fitted ~ am + cyl, data = xmarg, FUN = mean)
  am cyl   fitted
1  0   4 22.83306
2  1   4 27.96721
3  0   6 19.06359
4  1   6 21.35732
5  0   8 15.08720
6  1   8 14.64519

如您所见，点估计值现在大不相同，并且再次没有 SE 表。解决上述问题 1 和问题 2 可能会解决问题 3。

4

2 回答 2

3

对于这些问题，您需要预测包，它是margins的一部分。目前无法获得平均预测的标准误差，但您至少可以使用以下方法获得与 Stata 相同的平均预测。

关于 Statamargins命令的关键直觉如下：

margins x1

相当于

margins, at(x1 = (...))

其中...是 的所有可能值x1。这些表达式中的任何一个都会生成反事实数据集，其中x1对于数据中的所有情况都固定为给定值，然后在数据集的这个临时的反事实版本上执行模型预测。

该over()选项是一个子集过程：

margins, over(x1)

根据 的值拆分数据x1，然后对每个子集执行模型预测。您可以将其与此结合使用，at但考虑起来会有些奇怪。例如：

margins, over(x1) at(x2 = (1 2))

将所有观测值固定x2为 1，然后将数据拆分为x1，然后为每个子集生成预测，并对它们进行平均。然后它对x2所有观察值设置为 2 的反事实版本重复此操作。

在 R 中，将为您提供using参数prediction::prediction()的等价物。它还将通过将数据子集传递给参数来为您提供等价物。at()atover()data

因此，对于您的问题 2：

> prediction::prediction(x, at = list(am = c(0,1)))
Average predictions for 32 observations:
 at(am) value
      0 20.12
      1 20.05

对于您的问题 3：

> prediction::prediction(x, at = list(am = c(0,1)), data = subset(mtcars, cyl == 4))
Average predictions for 11 observations:
 at(am) value
      0 26.62
      1 26.55
> prediction::prediction(x, at = list(am = c(0,1)), data = subset(mtcars, cyl == 6))
Average predictions for 7 observations:
 at(am) value
      0 20.08
      1 20.01
> prediction::prediction(x, at = list(am = c(0,1)), data = subset(mtcars, cyl == 8))
Average predictions for 14 observations:
 at(am) value
      0 15.03
      1 14.96

在这两种情况下，您都不能通过执行 a 并聚合预测来复制 Stata 的输出，predict(x)因为预测发生在反事实数据集上。

而且，目前尚未实施差异（截至 2018 年 8 月）。

于 2018-08-01T20:58:41.287 回答
1

我遇到了同样的问题，并找到了以下解决方法。线程当然是一个旧线程。但我认为将我的解决方案添加到此线程时会更容易找到。

我模拟了一个因变量的数据，该变量dv由变量level和treat以及它们的交互作用来解释。

  1. 数据模拟

    N <- 1000
uid <- rep(1:N)
treat <- rep(1:10, each = N/10)
level <- rep(1:100, each = N/100)
err <- rnorm(N, 0, 1)
hdv <- 40 + 2 * treat + .25 * level - .05 * treat * level + err
dv <- ifelse(hdv > 47, 1, 0)
dat <- data.frame(dv = dv, treat = treat, level = level, hdv = hdv)

  2. 估计

    由于因变量是二元的，我估计了一个 Logit 模型。众所周知，Logit 中的交互项（与任何非线性模型一样）不能直接解释。

    这就是为什么我想要“水平”而不是“治疗”的边际效应：

    logit <- glm(dv ~ treat*level, family = binomial(link = "logit"), data = dat)

  3. 边际效应

    在对数据进行子集化时，R 实际上可以通过置信区间恢复边际效应，如下所示：

    hmpr7 <- summary(margins(logit, variables = "level", data = dat[dat$treat == 7,]))

    以下是对所有治疗执行此操作的（有些涉及）方法：

    hmpr <- list()
for (i in 1:10) {
  hmpr[[i]] <- summary(margins(logit, variables = "level", data = dat[dat$treat == i,]))
}
# the result is a list. For further use it is transformed into a data.frame
mpr <- data.frame(matrix(unlist(hmpr), nrow=length(hmpr), byrow=T))
# in this process, all variables are classified as factors. This is changed here
mpr <- data.frame(lapply(mpr, function(x) as.numeric(as.character(x))))
# only the variables of interest for the graph are kept
mpr <- mpr[,c(2, 6, 7)]
# meaningful names are assigned to the variables
mpr <- setNames(mpr, c("pred", "lower", "upper")) 
# treatment classifier is added to rows
mpr$treat <- rep(1:10)

  4. 绘制结果（如在 Stata 中marginsplot）

    plot(mpr$pred ~ mpr$treat,
ylim = range(c(mpr$lower, mpr$upper)),
pch = 19, xlab = "treatment", ylab = "marginal effect + 95% CI",
main = "marginal effect of level per treatment")

arrows(mpr$treat, mpr$lower,
  mpr$treat, mpr$upper,
  length = .05, angle = 90, code = 3)

abline(h = 0, col = "red")
于 2019-07-20T17:46:00.443 回答