How to get the marginal effects after lm_robust() with clustered standard errors?

Question

I am running a regression with clustered standard errors by year. This is easy to do with Stata but I have to do it with R, so I run it using the lm_robust() function from the estimatr package. The problem is that I must now get the marginal effects of some variables but I cannot do it and I guess it is because of the cluster standard error. I followed what is on the manual for lm_robust() and I've seen they only used the margins command from the margins package for other functions without clustered standard errors... Does anyone have a clue on how can I get and plot the marginal effects?

set.seed(42)
library(fabricatr)
library(randomizr)
dat <- fabricate(
  N = 100,                        # sample size
  x = runif(N, 0, 1),             # pre-treatment covariate
  y0 = rnorm(N, mean = x),        # control potential outcome
  y1 = y0 + 0.35,                 # treatment potential outcome
  z = complete_ra(N),             # complete random assignment to treatment
  y = ifelse(z, y1, y0),          # observed outcome

  # We will also consider clustered data
  clust = sample(rep(letters[1:20], each = 5)),
  z_clust = cluster_ra(clust),
  y_clust = ifelse(z_clust, y1, y0)
)

Then when I run the regression with the lm_robust() function:

library(estimatr)
lmout_cl <- lm_robust(
  y_clust ~ z_clust + x,
  data = dat,
  clusters = clust
)

And finally, I try to get the margins...

library(margins)
mar_cl <- margins(lmout_cl)

But this results in an error:

Error in attributes(.Data) <- c(attributes(.Data), attrib) :'names' attribute 
[1] must be the same length as the vector [0]

Answer 1

Apologies for this bug which prevents margins() from working with lm_robust() objects with non-numeric clusters in estimatr versions 0.10 and earlier. This was created by the internal way both estimatr::lm_robust() and margins::margins() handle which variables are in the model.

The bug has since been solved and so you have two solutions within estimatr.

Let me first generate the data.

library(fabricatr)
library(randomizr)
dat <- fabricate(
  N = 100,
  x = runif(N),
  clust = sample(rep(letters[1:20], each = 5)),
  y_clust = rnorm(N),
  z_clust = cluster_ra(clust),
)

---

Get the latest version of estimatr (v0.11.0)

The dev version on https://declaredesign.org/r/estimatr has a fix for this bug, and it will be on CRAN in the next month or so.

install.packages("estimatr", dependencies = TRUE,
                 repos = c("http://r.declaredesign.org", "https://cloud.r-project.org"))
library(estimatr)
lmout_cl <- lm_robust(
  y_clust ~ z_clust + x,
  data = dat,
  clusters = clust
)
library(margins)
mar_cl <- margins(lmout_cl)

---

Use numeric clusters with CRAN version of estimatr (v0.10.0)

A workaround with the existing version of estimatr on CRAN is to use numeric clusters instead of character clusters

dat <- fabricate(
  N = 100,
  x = runif(N),
  clust = sample(rep(1:20, each = 5)),
  y_clust = rnorm(N),
  z_clust = cluster_ra(clust),
)
install.packages("estimatr")
library(estimatr)
lmout_cl <- lm_robust(
  y_clust ~ z_clust + x,
  data = dat,
  clusters = clust
)
mar_cl <- margins(lmout_cl)