Keep in mind what Stata's margins command is doing by default (calculating "predictive margins"). The specific command you've typed is equivalent to margins, over(cyl) which is also a subsetting operation. This package calculates "marginal effects" (the equivalent of Stata's margins, dydx(*) command). This means the effect on the outcome of changing from the baseline/reference/intercept to a given value of the independent variable.

The "prediction" package provides functionality akin to margins, at(...):

> library("prediction")
> summary(prediction(fit))
 Prediction      SE     z         p  lower  upper
     0.4062 0.07402 5.488 4.056e-08 0.2612 0.5513
> summary(prediction(fit, at = list(cyl = c("4", "6", "8"))))
 at(cyl) Prediction     SE     z         p   lower  upper
       4     0.7273 0.1187 6.127 8.942e-10  0.4946 0.9599
       6     0.4286 0.1339 3.201 1.372e-03  0.1661 0.6910
       8     0.1429 0.1339 1.067 2.860e-01 -0.1196 0.4053

which is not a subsetting operation but provides similar inference to your example Stata code.

Thanks for clarifying, this is a big help. I'm trying to get predicted SEs from a glmer model on the probability scale. I noticed that "prediction" will compute the fit for merMod object but not the SEs. Any suggestions on where to turn for the SEs?

Thank you

There are some suggestions here: leeper/prediction#9

Thanks, appreciate it.

The "emmeans" package is also working