Replicating the Stata marginlist argument using R margins package?

Question

I am unable to replicate in R a particular use case of the Stata margins command: margins var1, over(var2) I've been trying to do so using the margins package in R.

To provide a reproducible example, I used the mtcars dataset and exported it from R into Stata so we are using the same dataset in both programs:

R code:

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

Stata code:

use "mtcars.dta", clear

Create an example linear regression model in both programs

Stata code:

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

R code:

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

The model output (not shown) is identical between the two programs

Compare the average marginal effects table for each variable in the model

Stata code and output:

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 code and output:

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

As you can see, these two outputs are very similar to one another, as expected using the R margins package.

Problem 1: Marginal predictions OVER the value of a variable

Stata Code and Output:

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 Code and Output:

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

In the two examples above, the marginal prediction is identical between R and Stata. However, is there a way (short of doing it by hand) to generate the delta-method standard error for each marginal prediction as is done in the Stata table above?

Problem 2: Marginal predictions for a specific variable:

Stata Code and Output:

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 Code and Output:

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

In this example, we are trying to replicate Stata’s “marginlist” argument in the margins command by subsetting the dataset after prediction. This does not seem to be the right way. How can we replicate these results from Stata?

Problem 3: Marginal prediction of one variable over the value of another

Replicating this result is my main goal!

Stata Code and output

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 Code and Output:

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

As you can see, the point estimates are now substantially different and again there is no SE table. Solving Problem 1 and Problem 2 above will likely allow the solution to Problem 3.

Answer 1

For these problems you want the prediction package, which is part of margins. It's not currently possible to get standard errors for average predictions, but you can at least get average predictions identical to Stata using the following.

The key intuition about Stata's margins command is the following:

margins x1

is equivalent to

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

where ... is all possible values of x1. Either of these expressions produces counterfactual datasets where x1 is fixed at a given value for all cases in the data, and then the model prediction is performed on this temporary, counterfactual version of the dataset.

The over() option is a subsetting procedure:

margins, over(x1)

splits the data based the value of x1 and then performs model prediction on each subset. You can combine this with at but it gets a bit odd to think about. For example:

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

fixes x2 to 1 for all observations, then splits the data by x1, then generates predictions for each subset, and averages them. Then it repeats this for a counterfactual version where x2 is set to 2 for all observations.

In R, prediction::prediction() will give you equivalents of the at() using the at argument. And it will also give you equivalents of over() by passing subsets of data to the data argument.

So, for your Problem 2:

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

And for your Problem 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

In neither of these cases can you replicate Stata's output by just doing a predict(x) and aggregating the predictions because the predictions are occurring on counterfactual datasets.

And, again, variances are not currently implemented (as of August 2018).

Answer 2

I had the same problem, and have found the following workaround. The thread is of course an old one. But I thought my solution would be easier to find when added to this thread.

I have simulated data of a dependent variable dv that is explained by variables level and treat, as well as their interaction.

1. Simulation of data

   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. Estimation

   As the dependent variable is binary, I estimate a Logit model. As is well understood, interaction terms in Logit (as in any non-linear model) cannot be directly interpreted.

   This is why I want marginal effects of "level" over "treat":

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

3. Marginal Effects

   R can actually recover marginal effects with confidence intervals when subsetting the data, as in:

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

   The following is a (somewhat involved) way of doing this for all treatments:

   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. Plotting the result (as in Stata's 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")