# Pairwise marginal effects at specific quintiles STATA-like in R's marginaleffects

I would like to estimate marginal effects of two different variables. In STATA, I have a two-way fixed effects model of type:

xtset region year
reg y x1 x2 c.x1#c.x2 x3 x4 i.region i.year [aw=pop], robust


In turn, I create the margins based on a specific set of quintiles from the variables, obtaining a 25-rows matrix of pairwise estimates.

margins, at(x1 =(0.02033865 0.0737 0.25173 1.04338 40.41884) x2 =(0.69 1.60 2.39 3.36 6.91))


What would be the correct translation of this code into R, using the marginaleffects package?

Here is the output of the STATA prediction applied to my dataset.

Below, instead, I created a reproducible example in both STATA and R, in order to make testing easier.

STATA:

set more off

cd "/your/dir"
use "https://dss.princeton.edu/training/Panel101.dta"

* Run regression with the interaction
xtset country year

reg y x1 x2 c.x1#c.x2 x3 i.country i.year, robust

centile x1 , centile(5 25 50 75 95)
centile x2 , centile(5 25 50 75 95)

* Create the margins based on distribution above
margins, at(x1 =(-.16 .32 .64 1.1 1.4) x2 =(-1.54 -1.22 -.46 1.61 1.80)) saving(predictions_test, replace)


Same code for R, where I am not sure the marginaleffects function is called correctly in order to replicate STATA's:

R:

library(tidyverse)
library(sandwich)
library(lmtest)
library(plm)
library(multiwayvcov)
library(margins)
library(MASS)
library(marginaleffects)
library(data.table)
library(MASS)
library(jtools) # for clustered SE
library(dplyr)

setwd('/your/dir/')

# Load Data -----------------------------------------------------------

library(foreign)
library(plm)

setkeyv(data,
c("country", "year")
)

m = "x1 + x2 + I(x1 * x2) + x3 + factor(country) + factor(year)"
f = as.formula(paste("y ~", m, collapse=" + "))
ols = lm(f,
# weights = pop,
data = data
)

x1quants = c(-.16, .32, .64, 1.1, 1.4)
x2quants = c(-1.54, -1.22, -.46, 1.61, 1.80)

pred = predictions(ols,
by = c("x1", "x2"),
variables = c("x1", "x2"),
newdata=datagridcf(x1=x1quants, x2=x2quants),
# cross=T
)

write.csv(pred, file='prediction_tests_R.csv', row.names = F)

• Can you include the stata output? I'm not familliar with stata and would like to make sure I can give you what you're looking for Commented Mar 29, 2023 at 12:00
• Yes, Stata output would be very useful. In the meantime, you can also read this page, where I do side-by-side comparisons between marginaleffects and Stata: vincentarelbundock.github.io/marginaleffects/articles/… Commented Mar 29, 2023 at 13:48
• You don't need the variables argument in the predictions() call. I can't replicate your results because I don't have access to Stata and because your two models don't match (no weights in R model). Have you read the link above and been able to exactly replicate the Stata examples? In particular, these: vincentarelbundock.github.io/marginaleffects/articles/… Commented Mar 29, 2023 at 17:30
• Thanks a lot @Vincent . I had misunderstood the role of "variables" in together with "by" argument for predictions. Now the estimates are much more similar to STATA's one, although still computationally more expensive. Thanks for the tip (and the great package!) Commented Mar 30, 2023 at 8:57

From your description of the output, I'm lead to believe the right approach is

library(tidyverse)
library(marginaleffects)

x1 <- rnorm(1000, 2, 1)
x2 <- rnorm(1000, -2, 1)
y <- 2*x1 - x2 + 0.5*x1*x2 + rnorm(1000)

d <- tibble(x1, x2, y)

fit <- lm(y~x1*x2, data=d)

comparisons(
fit,
variables = c('x1', 'x2'),
newdata = datagrid(x1=quantile(x1), x2=quantile(x2)),
cross = T
)


Here, I have specified that the marginal effects should be computed at the quantile of each variable (you're free to specify your own values). This returns 25 marginal effects estimated at each combination of the variables specified in newdata = datagrid(x1=quantile(x1), x2=quantile(x2)).

If you edit your question to include stata output and a reproducible example in R, then I can verify this is the right appraoch.

• I think that in Stata, margins will nearly always return an average estimate, so the equivalent would be avg_comparisons(). Also, the at argument in Stata builds what marginaleffects calls "counterfactual" dataset, by replicating the full observed data for each combination of values in at. So I think the analogous call in marginaleffects would use the datagridcf() function (note the "cf"). Commented Mar 29, 2023 at 13:50
• Demetri and @Vincent , thanks a lot for the feedback. I have updated the question with both the STATA output and a reproducible example for testing, I am not sure the call is correct. It doesn't matter now to programmatically include the computation of the quantiles, could be other user-generated inputs, I mostly wanna make sure the translation between STATA and R is correct. Commented Mar 29, 2023 at 15:14
• @Vincent side note: I have tried avg_comparisons in marginaleffects with datagridcf() on a 8GB ram laptop with a dataset of roughly 4500 observations. It seems to saturate the RAM, and very slow computation. Seems to be very computationally inefficient with respect to margins, at() counterpart in STATA, which takes only a few seconds? Commented Mar 29, 2023 at 16:53

As per @Vincent 's suggestion, the correct translation for the pairwise STATA behavior would be:

pred = predictions(ols,
by = c("x1", "x2"),
newdata=datagridcf(x1=x1quants, x2=x2quants),
)