如何用R/marginaleffects或Stata/margins转换边际OR至不同基线风险的风险差？

我在R中拟合了logistic回归模型，关注的协变量为离散型治疗效应（数据来自对比两种治疗的随机对照试验）。已掌握使用R包marginaleffects的G公式估计总体平均（边际）比值比（OR）的方法，但实际需求是估计不同基线风险下的比例差（即风险差、绝对风险降低/增加）。不确定该转换是否可行，有观点认为可通过Stata的margins命令实现，但查阅文档后未找到具体方法。现提出两个问题：

• 该转换是否可行？
• 若可行，如何通过R的marginaleffects或Stata的margins实现？

示例代码 #

# 加载包
library(marginaleffects)
library(rms)
remotes::install_github("jamesotto852/ggdensity")

# 加载数据
load(url(
  "https://github.com/vincentarelbundock/modelarchive/raw/main/data-raw/gusto.rda"
))

gusto <- subset(gusto, tx %in% c("tPA", "SK"))
gusto$tx <- factor(gusto$tx, levels = c("tPA", "SK"))

# 拟合模型
mod <- glm(
  day30 ~ tx + rcs(age, 4) + Killip + pmin(sysbp, 120) + lsp(pulse, 50) +
    pmi + miloc + sex, family = binomial(),
  data = gusto)

# 基于G估计的总体平均（边际）比值比
mOR = 
  comparisons(
  mod,
  variables = "tx",
  transform_pre = "lnoravg",
  transform_post = exp) |>
  summary()

mOR
#>   Term                 Contrast Effect Std. Error z value   Pr(>|z|) 2.5 %
#> 1   tx ln(odds(SK) / odds(tPA))  1.192    0.04489   26.54 < 2.22e-16 1.104
#>   97.5 %
#> 1   1.28
#> 
#> Model type:  glm 
#> Prediction type:  response 
#> Post-transformation:  exp

# 基于G估计的总体平均（边际）风险差
mRD = 
  comparisons(
    mod,
    variables = "tx") |>
  summary()

mRD
#>   Term Contrast  Effect Std. Error z value   Pr(>|z|)    2.5 % 97.5 %
#> 1   tx SK - tPA 0.01108   0.002766   4.005 6.1955e-05 0.005658 0.0165
#> 
#> Model type:  glm 
#> Prediction type:  response

# 个体水平风险差
uRD = 
  comparisons(
    mod,
    variables = "tx")

uRD |> 
  ggplot() +
  aes(x = predicted_hi, y = comparison) +
  ggdensity::geom_hdr_points() +
  scale_color_viridis_d(option = "magma", begin = .8, end = 0, direction = -1) +
  labs(x = "基线风险", y = "风险差") +
  theme_minimal() 

图表说明 #

现有超过3万个数据点，"probs"标识了多数点的分布位置（依次为50%、80%等），图表展示了基线风险与个体水平风险差的关系，数据点按密度分层呈现。

会话信息 #

sessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Catalina 10.15.7
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] rms_6.3-0                  SparseM_1.81              
#> [3] Hmisc_4.7-0                ggplot2_3.3.6             
#> [5] Formula_1.2-4              survival_3.3-1            
#> [7] lattice_0.20-45            marginaleffects_0.7.0.9000
#> 
#> loaded via a namespace (and not attached):
#>  [1] viridisLite_0.4.0   splines_4.1.2       assertthat_0.2.1   
#>  [4] highr_0.9           latticeExtra_0.6-29 yaml_2.3.5         
#>  [7] remotes_2.4.2       pillar_1.8.0        backports_1.4.1    
#> [10] quantreg_5.93       glue_1.6.2          digest_0.6.29      
#> [13] RColorBrewer_1.1-3  checkmate_2.1.0     colorspace_2.0-3   
#> [16] sandwich_3.0-2      htmltools_0.5.2     Matrix_1.4-1       
#> [19] pkgconfig_2.0.3     purrr_0.3.4         mvtnorm_1.1-3      
#> [22] scales_1.2.0        jpeg_0.1-9          MatrixModels_0.5-0 
#> [25] htmlTable_2.4.0     tibble_3.1.8        farver_2.1.1       
#> [28] generics_0.1.3      TH.data_1.1-1       withr_2.5.0        
#> [31] nnet_7.3-17         cli_3.3.0           magrittr_2.0.3     
#> [34] polspline_1.1.20    evaluate_0.15       fs_1.5.2           
#> [37] fansi_1.0.3         nlme_3.1-158        MASS_7.3-57        
#> [40] foreign_0.8-82      tools_4.1.2         data.table_1.14.3  
#> [43] lifecycle_1.0.1     multcomp_1.4-19     stringr_1.4.0      
#> [46] munsell_0.5.0       reprex_2.0.1        cluster_2.1.3      
#> [49] ggdensity_0.1.0.900 compiler_4.1.2      rlang_1.0.4        
#> [52] grid_4.1.2          rstudioapi_0.13     htmlwidgets_1.5.4  
#> [55] labeling_0.4.2      base64enc_0.1-3     rmarkdown_2.14     
#> [58] gtable_0.3.0        codetools_0.2-18    DBI_1.1.3          
#> [61] curl_4.3.2          R6_2.5.1            gridExtra_2.3      
#> [64] zoo_1.8-10          knitr_1.39          dplyr_1.0.9        
#> [67] fastmap_1.1.0       utf8_1.2.2          insight_0.18.0.4   
#> [70] stringi_1.7.6       vctrs_0.4.1         rpart_4.1.16       
#> [73] png_0.1-7           tidyselect_1.1.2    xfun_0.31

^{Created on 2022-08-15 by the reprex package (v2.0.1)}

内容的提问来源于stack exchange，提问作者Arthur Albuquerque