I estimated a probit model by using the glm function.
From the package lmtest is very easy to compute clustered standard errors by using the coeftest function.
To compute the average marginal effect, I type: AME=margins(model1, variables = c("x", "y")).
How can I compute the clustered errors of AME?
The margins package is now in “maintenance mode”, which means that it is not being actively developed. I recommend you use the marginaleffects package instead, which was designed as a successor (disclaimer: I am the author): https://marginaleffects.com
You can use clustered standard errors by specifying a one-sided formula in the vcov argument:
library(marginaleffects)
mod <- glm(vs ~ mpg + hp + factor(cyl), data = mtcars, family = binomial)
avg_slopes(mod,
variables = c("mpg", "hp"),
vcov = ~cyl)
#
# Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
# hp -0.00282 0.000476 -5.937 <0.001 28.4 -0.00376 -0.00189
# mpg -0.01242 0.015213 -0.817 0.414 1.3 -0.04224 0.01739
#
# Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Compare to classical IID standard errors:
avg_slopes(mod, variables = c("mpg", "hp"))
#
# Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
# hp -0.00282 0.00228 -1.236 0.216 2.2 -0.0073 0.00165
# mpg -0.01242 0.01681 -0.739 0.460 1.1 -0.0454 0.02053
#
# Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
See the documentation for vcov in ?avg_slopes for details.