Репликация аргумента Margins Stata с использованием пакета R Margins?

Я не могу воспроизвести в R конкретный вариант использования команды Stata margins: margins var1, over(var2) Я пытался сделать это с помощью пакета margins в R.

Чтобы предоставить воспроизводимый пример, я использовал набор данных mtcars и экспортировал его из R в Stata, поэтому мы используем один и тот же набор данных в обеих программах:

Код R:

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

Код статистики:

use "mtcars.dta", clear

Создайте пример модели линейной регрессии в обеих программах

Код статистики:

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

В этом примере мы пытаемся воспроизвести аргумент «маржинального списка» Stata в команде margins, разделив набор данных на подмножество после прогнозирования. Кажется, это неправильный путь. Как мы можем воспроизвести эти результаты 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.


r prediction stata marginal-effects
person ecidonex    schedule 31.07.2017    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Для решения этих проблем вам понадобится пакет прогнозов, который является частью полей. В настоящее время невозможно получить стандартные ошибки для средних прогнозов, но вы можете, по крайней мере, получить средние прогнозы, идентичные 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

Ни в одном из этих случаев вы не можете воспроизвести вывод Stata, просто выполнив predict(x) и агрегируя прогнозы, потому что прогнозы происходят на контрфактических наборах данных.

И, опять же, отклонения в настоящее время не реализованы (по состоянию на август 2018 г.).

person Thomas    schedule 01.08.2018
comment
Спасибо за подробный ответ! Цените помощь - person ecidonex; 18.08.2018

arrow_upward
1
arrow_downward

У меня была такая же проблема, и я нашел следующее решение. Нить конечно старая. Но я подумал, что мое решение будет легче найти, если добавить в эту ветку.

Я смоделировал данные зависимой переменной 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. Построение результата (как в 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")
person Christoph Engel    schedule 20.07.2019