Back to Index

Simulate Data

library(data.table)

set.seed(123456)
n <- 5000
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))
  , var1n = rnorm(n, 0, 1)
  , or2 = rep(1.1, n)
  , var2 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.4, 0.6))
  , var2n = rnorm(n, 0, 2)
  , or3 = rep(1.2, n)
  , var3 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.2, 0.8))
  , var3n = rnorm(n, 0, 2)
  , or4 = rep(1.5, n)
  , var4 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.3, 0.7))
  , var4n = rnorm(n, 0, 2)
  , or5 = rep(1.7, n)
  , var5 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.5, 0.5))
  , var5n = rnorm(n, 0, 2)
  , or6 = rep(2, n)
  , var6 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.4, 0.6))
  , var6n = rnorm(n, 0, 2)
  , or7 = rep(5, n)
  , var7 = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.1, 0.9))
  , var7n = rnorm(n, 0, 2)
)


dt <- dt[, odds0 := p0 / (1 - p0)
         ][, log_odds := log(odds0) +
                 var1 * log(or1) + var1n * log(or1) + 
                 var2 * log(or2) + var2n * log(or2) + 
                 var3 * log(or3) + var3n * log(or3) + 
                 var4 * log(or4) + var4n * log(or4) + 
                 var5 * log(or5) + var5n * log(or5) + 
                 var6 * log(or6) + var6n * log(or6) + 
                 var7 * log(or7) + var7n * 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)]

unique(dt[, .(or1, or2, or3, or4, or5, or6, or7)]) %>% prt(caption = "Variables with Odds Ratios")

GLM

m <- glm(outcome ~ var1 + var2 + var3 + var4 + var5 + var6 + var7 +
             var1n + var2n + var3n + var4n + var5n + var6n + var7n
       , data = dt
       , family = binomial
         )

library(sjPlot)

tab_model(m)

CART: the Full Model

library(rpart)
library(rpart.plot)
library(RColorBrewer)
library(rattle)

predictors <- c(
    "var1"
  , "var2"
  , "var3"
  , "var4"
  , "var5"
  , "var6"
  , "var7"
  , "var1n"
  , "var2n"
  , "var3n"
  , "var4n"
  , "var5n"
  , "var6n"
  , "var7n"
)

frml <- Wu::wu_formula(outcome = "outcome", predictors = predictors)



set.seed(123456)
tr <- rpart(
    frml
  , data = dt
  , method = "class"
  , model = TRUE
  , x = TRUE
  , y = TRUE
  , parms = list(split = "information")
  , control = rpart.control(cp = 0
                          , xval = 20
                          , maxdepth = 30
                          , minsplit = 10
                          , minbucket = 5
                            )
)

rpart.plot(tr, type = 2, extra = 106, tweak = 2, under = TRUE)

Prune the CART Model

tr$cptable %>% prt(caption = "CP Table")

plotcp(tr)

bestcp <- tr$cptable[which.min(tr$cptable[,"xerror"]),"CP"]
## bestcp <- tr$cptable[12, "CP"]
trp <- prune(tr, cp = bestcp)
rpart.plot(trp, type = 2, extra = 106, tweak = 2, under = TRUE)

## fancyRpartPlot(trp, uniform=TRUE, main="Pruned Classification Tree")

bestcp2 <- tr$cptable[9,"CP"]


fancyRpartPlot(trp, caption = "fancyRpartPlot")

AUC on the Training Data

library(pROC)
pred <- predict(object = trp, newdata = dt, type = "prob")
r <- roc(dt$outcome, pred[, 2], ci = TRUE, direction = "<")
plot(r)

plt_roc(r) %>% ann("Pruned CART")

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] pROC_1.18.0 rattle_5.5.1 bitops_1.0-7
[4] tibble_3.1.7 RColorBrewer_1.1-3 rpart.plot_3.1.1
[7] rpart_4.1-15 sjPlot_2.8.10 Wu_0.0.0.9000
[10] flexdashboard_0.5.2 lme4_1.1-29 Matrix_1.4-0
[13] mgcv_1.8-38 nlme_3.1-152 png_0.1-7
[16] scales_1.2.0 nnet_7.3-16 labelled_2.9.1
[19] kableExtra_1.3.4 plotly_4.10.0 gridExtra_2.3
[22] ggplot2_3.3.6 DT_0.23 tableone_0.13.2
[25] magrittr_2.0.3 lubridate_1.8.0 dplyr_1.0.9
[28] plyr_1.8.7 data.table_1.14.2 rmdformats_1.0.4
[31] knitr_1.39

loaded via a namespace (and not attached): [1] TH.data_1.1-1 minqa_1.2.4 colorspace_2.0-3 ellipsis_0.3.2
[5] sjlabelled_1.2.0 estimability_1.4 parameters_0.18.1 rstudioapi_0.13
[9] fansi_1.0.3 mvtnorm_1.1-3 xml2_1.3.3 codetools_0.2-18 [13] splines_4.2.0 sjmisc_2.8.9 jsonlite_1.8.0 nloptr_2.0.3
[17] ggeffects_1.1.2 broom_0.8.0 effectsize_0.7.0 compiler_4.2.0
[21] httr_1.4.3 sjstats_0.18.1 emmeans_1.7.5 backports_1.4.1
[25] assertthat_0.2.1 fastmap_1.1.0 lazyeval_0.2.2 survey_4.1-1
[29] cli_3.3.0 htmltools_0.5.3 tools_4.2.0 gtable_0.3.0
[33] glue_1.6.2 Rcpp_1.0.8.3 jquerylib_0.1.4 vctrs_0.4.1
[37] svglite_2.1.0 crosstalk_1.2.0 insight_0.18.0 xfun_0.31
[41] stringr_1.4.0 rvest_1.0.2 lifecycle_1.0.1 klippy_0.0.0.9500 [45] MASS_7.3-54 zoo_1.8-10 hms_1.1.1 sandwich_3.0-2
[49] yaml_2.3.5 sass_0.4.1 stringi_1.7.8 highr_0.9
[53] bayestestR_0.12.1 boot_1.3-28 rlang_1.0.4 pkgconfig_2.0.3
[57] systemfonts_1.0.4 evaluate_0.15 lattice_0.20-45 purrr_0.3.4
[61] htmlwidgets_1.5.4 tidyselect_1.1.2 bookdown_0.27 R6_2.5.1
[65] generics_0.1.2 multcomp_1.4-19 DBI_1.1.2 pillar_1.7.0
[69] haven_2.5.0 withr_2.5.0 survival_3.2-13 datawizard_0.4.1 [73] performance_0.9.1 modelr_0.1.8 crayon_1.5.1 utf8_1.2.2
[77] rmarkdown_2.14 grid_4.2.0 forcats_0.5.1 digest_0.6.29
[81] webshot_0.5.3 xtable_1.8-4 tidyr_1.2.0 munsell_0.5.0
[85] viridisLite_0.4.0 bslib_0.3.1 mitools_2.4