eci*_*nex 5 r prediction stata marginal-effects
我无法在R中复制Stata margins命令的特定用例：
margins var1, over(var2)
我一直在尝试使用marginsR中的包进行复制。
为了提供可重现的示例，我使用了mtcars数据集并将其从R导出到Stata中，因此我们在两个程序中都使用了相同的数据集：
R代码：
library(foreign)
library(margins)
write.dta(mtcars, “mtcars.dta")
Stata代码：
use "mtcars.dta", clear
在两个程序中创建示例线性回归模型
Stata代码：
quietly regress mpg cyl i.am c.wt##c.hp
R代码：
x <- lm(mpg ~ cyl + factor(am) + hp * wt, data = mtcars)
两个程序之间的模型输出（未显示）相同
比较模型中每个变量的平均边际效应表
Stata代码和输出：
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
如您所见，这两个输出非常相似，就像使用R margins包所期望的那样。
问题1：边际预测超过变量的值
Stata代码和输出：
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之间的边际预测相同。但是，是否有一种方法（不需要手动完成），可以像上面的Stata表中那样为每个边际预测生成增量法标准误差？
问题2：特定变量的边际预测：
Stata代码和输出：
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：一个变量对另一个变量的边际预测
复制此结果是我的主要目标！
Stata代码和输出
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。
对于这些问题，您需要预测包，它是margins的一部分。目前无法获得平均预测的标准误差，但您至少可以使用以下方法获得与 Stata 相同的平均预测。
关于 Stata 命令的关键直觉margins如下：
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 中，prediction::prediction()将为您提供at()使用at参数的等效项。它还会为您提供相当于over()将数据子集传递给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
在这两种情况下，您都无法通过执行predict(x)并聚合预测来复制 Stata 的输出，因为预测是在反事实数据集上发生的。
而且，目前尚未实施差异（截至 2018 年 8 月）。