# how to replicate the stata "margins, atmeans" command in R with the margins library [closed]

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$$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)

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

• I think the quotation that you give outlines exactly the steps that you need to follow. Suppose your data lives in data and var is the variable that you're interested in. You'd replace the raw values of var with the average of the values of var. That's what the line data$var <- mean(data$var, na.rm = TRUE) accomplishes. Then you just supply data to margins. Can you clarify what part you're having trouble with?
– Sycorax
Commented Feb 6, 2019 at 1:32
• hi, i'm not understanding. could you take the minimal example i've provided and match the ucla-published numbers using R? Commented Feb 6, 2019 at 19:00

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)
1 0.1174438 0.06010156 0.2013822   male 52.23
2 0.2832163 0.19213931 0.3900793 female 52.23

• Sidenote... the UCLA website is using the term 'predicted probability'. To me this sounds a bit as a confusing term. It confuses me when I start to think what a 'predicted probability' could literally mean... Is it a prediction for the mean-parameter (ie probability) in a Bernouilli/binomial distributed random variable?. It is at least not a prediction interval, it is a (confidence) interval expressing the error of the estimates of the mean. Commented Feb 7, 2019 at 12:19

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

library(foreign)

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