how to replicate the stata "margins, atmeans" command in R with the margins library

Question

i'm trying to replicate the output of margins female, atmeans in R shown here:

https://stats.idre.ucla.edu/stata/dae/using-margins-for-predicted-probabilities/

i can re-create the same setup shown on the ucla page with

library(foreign)
library(margins)
x <- read.dta( "https://stats.idre.ucla.edu/stat/data/hsbdemo.dta" )
x$honors <- as.numeric( x$honors == 'enrolled' )

and this code matches the logit honors i.female read line on ucla's page

this_model <- glm( honors ~ female + read , data = x , family = binomial() )
summary(this_model)

from here, i'm confused about how to modify

summary(margins(this_model))

so that i'm replicating the stata atmeans parameter. my goal is to reproduce the output shown on the ucla page, but in R instead of stata:

margins female, atmeans
Adjusted predictions                              Number of obs   =        200
Model VCE    : OIM
Expression   : Pr(honors), predict()
at           : 0.female        =        .455 (mean)
               1.female        =        .545 (mean)
               read            =       52.23 (mean)

         |            Delta-method
         |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]

-------------+----------------------------------------------------------------
      female |
          0  |   .1127311   .0350115     3.22   0.001     .0441097    .1813524
          1  |   .2804526   .0509114     5.51   0.000     .1806681    .3802371

the R margins library help pages discusses the atmeans command, but it's not obvious to me how to implement this:

atmeans: calculate marginal effects at the mean (MEMs) of a dataset rather than the default behavior of calculating average marginal effects (AMEs)

Stata’s atmeans argument is not implemented in margins() for various reasons, including because it is possible to achieve the effect manually through an operation like data$var <- mean(data$var, na.rm = TRUE) and passing the modified data frame to margins(x, data = data).

Answer 1

The function margin only reports the effect.

> mar <-  margins(this_model, 
+                          data = x,
+                          type = "response",
+                          var = c("female"),
+                          at = list(read = mean(x$read))
+ )
> summary(mar)
       factor    read    AME     SE      z      p  lower  upper
 femalefemale 52.2300 0.1677 0.0585 2.8673 0.0041 0.0531 0.2824

You can see that this is sort of equivalent to the stata code, for instance $0.2804526-0.1127311 = 0.1677215 \approx 0.1677$

---

This is mentioned here

https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html

The functionality of Stata's command to produce predictive margins is not ported, as this is easily obtained from the prediction package

---

So you can use the prediction function. For example:

newdata <- data.frame(female = c("male","female"),
                        read = rep(mean(x$read),2))
p <- predict(object = this_model,
             newdata = newdata,
             type = "response",
             se.fit = TRUE)
mult <- qnorm(0.5*(1-0.95))
out <- cbind(p$fit,
             p$se.fit,
             p$fit+p$se.fit*mult,
             p$fit-p$se.fit*mult)
rownames(out) <- levels(newdata$female)[newdata$female]
colnames(out) <- c("margin", "Std.Err", "lower 95% conf", "upper 95% conf")

with output

> out
          margin    Std.Err lower 95% conf upper 95% conf
male   0.1127310 0.03501142     0.04410991      0.1813522
female 0.2804526 0.05091129     0.18066832      0.3802369

In the above the confidence intervals are computed manually since the predict function does not do this for glm objects. I believe it is a bit tricky to compute the confidence interval for glm objects since the estimate of the variance (and maybe also the effective degrees of freedom?) is not straightforward. The method above, which matches the stata output uses a computation of 95% confidence intervals based on $\pm 1.96$ times the standard error and is based on a normal distribution (for the estimate of the mean/margin with known variance). This works because the data is modeled with a binomial distribution (where estimate of mean and variance are related). But it does not work in general for GLM and sometimes you should use a t-distribution which gives a slightly wider interval (or much larger interval when the sample size is small).

There is alternatively the glm.predict function which computes all this stuff based on a simulation (what kind of simulation I do not know).

> library(glm.predict)
> predicts(this_model,"F;mean",sim.count = 1000000, set.seed = 1)
       mean      lower     upper female  read
1 0.1174438 0.06010156 0.2013822   male 52.23
2 0.2832163 0.19213931 0.3900793 female 52.23

Answer 2

following up with info to replicate the second half of this ucla page:

library(foreign)

x <- read.dta( "https://stats.idre.ucla.edu/stat/data/hsbdemo.dta" )
x$honors <- as.numeric( x$honors == 'enrolled' )
x$female <- as.numeric( x$female == 'female' )

this_model <- glm( honors ~ female + read , data = x , family = binomial() )


newdata <- expand.grid(female = mean(x$female),
       read = seq(20,70,10))
 p <- predict(object = this_model,
     newdata = newdata,
     type = "response",
     se.fit = TRUE)
 mult <- qnorm(0.5*(1-0.95))
 out <- cbind(p$fit,
             p$se.fit,
     p$fit+p$se.fit*mult,
     p$fit-p$se.fit*mult)
rownames(out) <- seq(20,70,10)
colnames(out) <- c("margin", "Std.Err", "lower 95% conf", "upper 95% conf")

matches stata margins, at(read=(20(10)70)) atmeans vsquish post output:

   margin      Std.Err lower 95% conf upper 95% conf
20 0.002226388 0.001962525   -0.001620091    0.006072867
30 0.009363869 0.006100713   -0.002593308    0.021321046
40 0.038500176 0.016282790    0.006586493    0.070413859
50 0.145023969 0.031199410    0.083874250    0.206173689
60 0.418114778 0.049886309    0.320339409    0.515890147
70 0.752714028 0.067025506    0.621346451    0.884081606