Back to Index

Simulate Data

set.seed(123456)
n <- 1000
dt <- data.table(
    p0 = rep(0.2, n)
  , or1 = rep(1, n)
  , var1 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.3, 0.7))
  , or2 = rep(1.1, n)
  , var2 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.4, 0.6))
  , or3 = rep(1.2, n)
  , var3 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.2, 0.8))
  , or4 = rep(1.5, n)
  , var4 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.3, 0.7))
  , or5 = rep(1.7, n)
  , var5 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.5, 0.5))
  , or6 = rep(2, n)
  , var6 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.4, 0.6))
  , or7 = rep(5, n)
  , var7 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.1, 0.9))
)


dt <- dt[, odds0 := p0 / (1 - p0)
         ][, log_odds := log(odds0) +
                 var1 * log(or1) +
                 var2 * log(or2) +
                 var3 * log(or3) +
                 var4 * log(or4) +
                 var5 * log(or5) +
                 var6 * log(or6) +
                 var7 * log(or7)
           ][, p := exp(log_odds)/ (1 + exp(log_odds))]

vsample <- function(p){
    sample(c(1, 0), size = 1, replace = TRUE, prob = c(p, 1 - p))
}

vsample <- Vectorize(vsample)

dt <- dt[, outcome := vsample(p)]

brms Package

• The function to report brms model needs to be modified to reflect the logistic model

library(brms)
library(bayestestR)
library(sjstats)

m <- brm(formula = outcome ~ var1 + var2 + var3 + var4 + var5 + var6 + var7,  
                   data = dt, 
                   family = bernoulli(link = "logit"),
                   warmup = 500, 
                   iter = 2000, 
                   chains = 2, 
                   inits = "0", 
                   cores = 2,
                   seed = 123)



report_brms <- function(obj, digits=2, equiv=0.1){
    pe <- unique(point_estimate(obj))
    colnames(pe) <- tolower(colnames(pe))
    pe <- pe[, -1]
    he <- unique(hdi(obj, ci=0.89))
    colnames(he) <- paste0("hdi_", tolower(colnames(he)))
    et <- unique(eti(obj, ci=0.95))
    colnames(et) <- paste0("eti_", tolower(colnames(et)))
    pd <- unique(pd(obj))
    colnames(pd) <- tolower(colnames(pd))
    pd <- pd[, 1:2]
    rtn <- cbind(pe, he, et, pd)
    rtn <- as.data.table(rtn)
    rtn$equivalence_l <- -equiv
    rtn$equivalence_r <- equiv
    rtn <- as.data.table(rtn)
    rtn <- rtn[, hdi := paste0(
                     as.character(round(median, digits))
                   , " ["
                   , as.character(round(hdi_ci_low, digits))
                   , ", "
                   , as.character(round(hdi_ci_high, digits))
                   , "]"
                 )]
    rtn <- rtn[, eti := paste0(
                     as.character(round(mean, digits))
                   , "["
                   , as.character(round(eti_ci_low, digits))
                   , ", "
                   , as.character(round(eti_ci_high, digits))
                   , "]"
                 )]
    rtn <- rtn[, pd_n := pd][, pd := Wu::percent(pd, digits=digits)]
    for(i in 1:nrow(rtn)){
        rp <- equivalence_test(
            obj
          , range = as.numeric(rtn[i, list(equivalence_l, equivalence_r)])
          , ci=0.89
          , parameters=rtn$parameter[i])
        rtn$rope_n[i] <- rp$ROPE_Percentage
        rtn$rope[i] <- Wu::percent(rp$ROPE_Percentage, digits=digits)
    }
    label(rtn$parameter) <- "Parameter"
    label(rtn$hdi) <- "Median [89% HDI]"
    label(rtn$eti) <- "Mean [95% ETI]"
    label(rtn$pd) <- "Probability of Direction"
    label(rtn$rope) <- "ROPE"
    invisible(rtn)
}

t <- report_brms(m, digits=2, equiv = 0.06)
t2 <- t[-1, ]


t2[, list(parameter, hdi, eti, pd, rope)] %>%
    prt(col.names=c("Parameter", "Median [89% HDI]", "Mean [95% ETI]", "Probability of Direction", "ROPE")) %>%
    footnote(symbol=c("HDI: High Density Credible Interval"
                      , "ETI: Equal-Tailed Credible Interval")
        , number=c(
                 "Probability of Direction (pd) describes effect existence. It represents the certainty associated with the most probable direction (positive or negative) of the effect."
                 , "pd <= 95%: uncertain"
                 , "pd > 95%: possibly existing"
                 , "pd > 97%: likely existing"
                 , "pd > 99%: probably existing"
                 , "pd > 99.9%: certainly existing"
          )
        , alphabet=c("ROPE: Region of Practical Equivalence with range $[-0.1*SD_{Y},0.1*SD_Y]$"
              , "> 99%: negligible (we can accept the null hypothesis)"
              , "> 97.5%: probably negligible"
              , "<= 97.5% & >= 2.5%: undecided significance"
              , "< 2.5%: probably significant"
              , "< 1%: significant (we can reject the null hypothesis)"
          )
             )

rmsb Package

• Multiple imputation

library(rmsb)

dt <- dt[, missing_indicator := sample(c(0, 1), size = n, replace = TRUE, prob = c(0.9, 0.1))
         ][, var3 := case_when(
                 missing_indicator == 1 ~ as.numeric(NA)
               , TRUE ~ var3)
           ]


dd <- datadist(dt)
options(datadist = "dd")
set.seed(123456)

mi <- aregImpute( ~ var1 + var2 + var3 + var4 + var5 + var6 + var7
               , data = dt
               , n.impute = 5
               , pr = FALSE)

m <- stackMI(outcome ~ var1 + var2 + var3 + var4 + var5 + var6 + var7
           , fitter=blrm
           , iter=2000
           ## , warmup=1000
           , chain=2
           ## , cores=2
           , data=dt
           , xtrans=mi
           , loo=TRUE
             )


get_post_m <- function(mod){
    mu_draws <- as.data.table(mod$draws)
    coef1 <- lapply(mu_draws, median)
    coef2 <- lapply(mu_draws, function(x){HPDint(x, prob = 0.95)[1]})
    coef3 <- lapply(mu_draws, function(x){HPDint(x, prob = 0.95)[2]})
    coef4 <- lapply(mu_draws, function(x){sum(x > 0) / length(x)})
    rtn <- cbind(names(coef1)
               , as.data.table(unlist(coef1))
               , as.data.table(unlist(coef2))
               , as.data.table(unlist(coef3))
               , as.data.table(unlist(coef4))
                 )
    colnames(rtn) <- c("coef", "median", "lower", "upper", "probgt0")
    rtn <- rtn[, hdi := ci_to_str(fit=exp(median), lower=exp(lower), upper=exp(upper), digits=2)
               ]
    invisible(rtn[-1, .(coef, hdi, probgt0)])
}

hdi <- get_post_m(m)

hdi %>% prt(caption="Multivariable Model"
                , col.names=c("Variable", "Odds Ratio (HDI)", "Probability > 1"))

R sessionInfo

sessionInfo()

R version 4.2.0 (2022-04-22) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 20.04.3 LTS

Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale: [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] rmsb_0.1.0 rms_6.3-0 SparseM_1.81
[4] Hmisc_4.7-0 Formula_1.2-4 survival_3.2-13
[7] lattice_0.20-45 sjstats_0.18.1 bayestestR_0.12.1
[10] brms_2.17.0 Rcpp_1.0.9 Wu_0.0.0.9000
[13] flexdashboard_0.5.2 lme4_1.1-29 Matrix_1.4-0
[16] mgcv_1.8-38 nlme_3.1-152 png_0.1-7
[19] scales_1.2.0 nnet_7.3-16 labelled_2.9.1
[22] kableExtra_1.3.4 plotly_4.10.0 gridExtra_2.3
[25] ggplot2_3.3.6 DT_0.23 tableone_0.13.2
[28] magrittr_2.0.3 lubridate_1.8.0 dplyr_1.0.9
[31] plyr_1.8.7 data.table_1.14.2 rmdformats_1.0.4
[34] knitr_1.39

loaded via a namespace (and not attached): [1] backports_1.4.1 systemfonts_1.0.4 igraph_1.3.4
[4] lazyeval_0.2.2 splines_4.2.0 crosstalk_1.2.0
[7] TH.data_1.1-1 rstantools_2.2.0 inline_0.3.19
[10] digest_0.6.29 htmltools_0.5.3 fansi_1.0.3
[13] checkmate_2.1.0 cluster_2.1.2 modelr_0.1.8
[16] RcppParallel_5.1.5 matrixStats_0.62.0 sandwich_3.0-2
[19] xts_0.12.1 svglite_2.1.0 prettyunits_1.1.1
[22] jpeg_0.1-9 colorspace_2.0-3 rvest_1.0.2
[25] mitools_2.4 haven_2.5.0 xfun_0.31
[28] callr_3.7.0 crayon_1.5.1 jsonlite_1.8.0
[31] zoo_1.8-10 glue_1.6.2 gtable_0.3.0
[34] emmeans_1.7.5 MatrixModels_0.5-0 webshot_0.5.3
[37] sjmisc_2.8.9 distributional_0.3.0 pkgbuild_1.3.0
[40] rstan_2.21.5 abind_1.4-5 mvtnorm_1.1-3
[43] DBI_1.1.2 miniUI_0.1.1.1 htmlTable_2.4.0
[46] performance_0.9.1 viridisLite_0.4.0 xtable_1.8-4
[49] klippy_0.0.0.9500 foreign_0.8-81 stats4_4.2.0
[52] StanHeaders_2.21.0-7 survey_4.1-1 datawizard_0.4.1
[55] htmlwidgets_1.5.4 httr_1.4.3 threejs_0.3.3
[58] RColorBrewer_1.1-3 posterior_1.2.2 ellipsis_0.3.2
[61] pkgconfig_2.0.3 loo_2.5.1 farver_2.1.1
[64] sass_0.4.1 utf8_1.2.2 effectsize_0.7.0
[67] tidyselect_1.1.2 rlang_1.0.4 reshape2_1.4.4
[70] later_1.3.0 munsell_0.5.0 tools_4.2.0
[73] cli_3.3.0 generics_0.1.3 sjlabelled_1.2.0
[76] broom_0.8.0 ggridges_0.5.3 evaluate_0.15
[79] stringr_1.4.0 fastmap_1.1.0 yaml_2.3.5
[82] processx_3.6.0 purrr_0.3.4 quantreg_5.93
[85] mime_0.12 xml2_1.3.3 compiler_4.2.0
[88] bayesplot_1.9.0 shinythemes_1.2.0 rstudioapi_0.13
[91] tibble_3.1.8 bslib_0.3.1 stringi_1.7.8
[94] highr_0.9 parameters_0.18.1 ps_1.7.0
[97] Brobdingnag_1.2-7 forcats_0.5.1 nloptr_2.0.3
[100] markdown_1.1 shinyjs_2.1.0 tensorA_0.36.2
[103] vctrs_0.4.1 pillar_1.8.0 lifecycle_1.0.1
[106] jquerylib_0.1.4 bridgesampling_1.1-2 estimability_1.4
[109] insight_0.18.0 httpuv_1.6.5 latticeExtra_0.6-29 [112] R6_2.5.1 bookdown_0.27 promises_1.2.0.1
[115] codetools_0.2-18 polspline_1.1.20 boot_1.3-28
[118] colourpicker_1.1.1 MASS_7.3-54 gtools_3.9.2.2
[121] assertthat_0.2.1 withr_2.5.0 shinystan_2.6.0
[124] multcomp_1.4-19 parallel_4.2.0 hms_1.1.1
[127] rpart_4.1-15 grid_4.2.0 tidyr_1.2.0
[130] coda_0.19-4 minqa_1.2.4 rmarkdown_2.14
[133] shiny_1.7.1 base64enc_0.1-3 dygraphs_1.1.1.6