我不太熟悉Stata中的over()选项。如果你描述了计算出的具体数学量，并提供了一个带有Stata输出和公共数据的可重现的例子，我也许能够确认这是否可以用R中的marginaleffects包来实现。
作为参考，您可以在此小插曲中看到marginaleffects和Stata的并行示例，用于许多等效调用： www.example.com
根据我从Stata文档中收集的信息，你发布的代码估计了结果的预测概率，在年龄值的反事实网格上估计，并在一个预测器的类别中取平均值。在marginaleffects中的类似例子是：

library(marginaleffects)
mod <- glm(vs ~ mpg * factor(cyl), data = mtcars, family = binomial)

# equivalent to ,at() option in Stata
avg_predictions(mod, by = "cyl", newdata = datagridcf(mpg = seq(10, 33, by = 5)))
 
#  cyl Estimate Pr(>|z|)    2.5 % 97.5 %
#    4 8.85e-01    0.269 1.70e-01  0.997
#    6 5.47e-01    0.835 1.70e-01  0.877
#    8 1.17e-09    0.998 2.22e-16  1.000
# 
# Columns: cyl, estimate, p.value, conf.low, conf.high

# equivalent to ,dydx at() option in Stata
avg_slopes(mod, by = "cyl", newdata = datagridcf(mpg = seq(10, 33, by = 5)))

#  Term          Contrast cyl  Estimate Std. Error         z Pr(>|z|)     2.5 %    97.5 %
#   cyl mean(6) - mean(4)   4 -3.72e-01   2.30e-01 -1.62e+00   0.1050 -8.22e-01  7.78e-02
#   cyl mean(6) - mean(4)   6 -3.72e-01   2.30e-01 -1.62e+00   0.1050 -8.22e-01  7.78e-02
#   cyl mean(6) - mean(4)   8 -3.72e-01   2.30e-01 -1.62e+00   0.1050 -8.22e-01  7.78e-02
#   cyl mean(8) - mean(4)   4 -8.82e-01   2.25e-01 -3.93e+00   <0.001 -1.32e+00 -4.42e-01
#   cyl mean(8) - mean(4)   6 -8.82e-01   2.25e-01 -3.93e+00   <0.001 -1.32e+00 -4.42e-01
#   cyl mean(8) - mean(4)   8 -8.82e-01   2.25e-01 -3.93e+00   <0.001 -1.32e+00 -4.42e-01
#   mpg mean(dY/dX)         4  4.18e-03   3.17e-02  1.32e-01   0.8949 -5.79e-02  6.62e-02
#   mpg mean(dY/dX)         6 -5.14e-02   2.92e-02 -1.76e+00   0.0785 -1.09e-01  5.85e-03
#   mpg mean(dY/dX)         8 -3.19e-21   2.25e-06 -1.42e-15   1.0000 -4.41e-06  4.41e-06
# 
# Columns: term, contrast, cyl, estimate, std.error, statistic, p.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo

编辑：请注意，由于逆变换和平均的顺序，Stata和R之间可能存在一些数值差异。您可能想尝试type参数的不同值（“link”或“response”）。