Back to Index

Simulate Data

• Assume a severity variable increases the treatment selection by an odds ratio of 1.5
• The baseline probability of getting treatment is 0.2
• The severity variable also increase the outcome risk by an odds ratio of 2
• The treatment reduces the outcome risk by by an odds ratio of 0.5
• The baseline probaiblity of outcome (death) is 0.2

set.seed(123456)
n <- 2000
dt <- data.table(severity = rnorm(n), covar1 = rnorm(n), covar2 = rnorm(n), covar3 = rnorm(n), covar_trt = rnorm(n))
p0 <- 0.2
or_trt <- 1.5
dt <- dt[, odds_trt_0 := p0 / (1 - p0)
         ][, odds2 := 1.5
           ][, odds3 := 2.5
             ][, log_odds_trt := log(odds_trt_0) + severity * log(or_trt) + covar_trt * log(2.5)
           ][, p_trt := exp(log_odds_trt) / (1 + exp(log_odds_trt))]

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

vsample <- Vectorize(vsample)

dt <- dt[, trt := vsample(p_trt)]

p <- 0.2
dt <- dt[, odds_outcome_0 := p / (1 - p)
         ][, log_odds_outcome := log(odds_outcome_0) + trt * log(0.5) + severity * log(2) + covar1 * log(2.5) + covar2 * log(1.5)
           ][, p_outcome := exp(log_odds_outcome) / (1 + exp(log_odds_outcome))]


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


Vars <- c("severity"
        , "covar1"
        , "covar2"
        , "covar3"
        , "covar_trt")

chr <- tbl1n(data = dt
           , vars = Vars
           , factorVars = NULL
           , strata = "trt"
           , test = TRUE
             )

chr[,] %>% prt()

Propensity Score Method

m_trt <- glm(trt ~ severity + covar1 + covar2 + covar3 + covar_trt, data = dt, family = binomial(logit))

library(sjPlot)

tab_model(m_trt)

dt$pred <- m_trt$fitted.values

p1 <- dt[trt == 1]$pred
p0 <- dt[trt == 0]$pred

plot_ly(alpha = 0.5, xbins = list(start = 0, end = 1, size = 0.02)) %>%
  add_histogram(x = ~ p1
              , name = "Treated"
              , inherit = TRUE
              ## , xbins = seq(0, 1, 0.05)
                ) %>%
  add_histogram(x = ~ p0
              , name = "Not Treated"
              , inherit = TRUE
              ## , xbins = seq(0, 1, 0.05)
                ) %>%
    layout(barmode = "overlay"
         , xaxis = list(title = paste0("Predicted Probabilities of Being Treated"),
                      zeroline = FALSE),
         yaxis = list(title = "Count",
                      zeroline = FALSE))

ggplot(data=dt, aes(x=pred, colour=factor(trt))) +
  geom_density(alpha=0.5) +
  labs(y="Density",x="pred") +
  scale_y_continuous(breaks=NULL,label=c("")) +
    theme(legend.position=c(0.8,0.2))

font_tick_size <- 10
font_axis_size <- 10
size_i <- 1
ggplot(data = dt, aes(x = pred, fill = factor(trt))) +
  geom_histogram(position = "dodge", breaks = seq(0,1,0.1), closed = "left", alpha = 1) +
  scale_fill_manual(values=c("#333333", "#999999")) +
  labs(x = "Predicted Probability of Being Treated"
     , y = "Count"
     , fill = "Treatment") +
  scale_x_continuous(breaks = seq(0,1,0.1), limits = c(0,1.05), expand = c(0, 0)) +
  scale_y_continuous(breaks = seq(0,100,20), limits = c(0,100), expand = c(0, 0)) +
        theme(
            axis.text.x = element_text(size = font_tick_size)
          , axis.text.y = element_text(size = font_tick_size)
          , axis.title.x = element_text(size = font_axis_size)
          , axis.title.y = element_text(size = font_axis_size)
         ## ,  plot.margin = unit(0.2, "cm")
           , axis.line = element_line(colour = "#111111")         
         ## , axis.ticks = element_blank()
         ## , panel.grid.major = element_line(colour = "#DDDDDD", size = size_i)
           ## , panel.grid.minor = element_blank()
        , panel.background = element_blank(),
          legend.position = c(.85, .85),
          )

m_outcome_mv <- glm(outcome ~ severity + trt + covar1 + covar2 + covar3 + covar_trt, data = dt, family = binomial(logit))

tab_model(m_outcome_mv)

library(mgcv)

m <- bam(outcome ~ severity + trt + covar1 + covar2 + covar3 + covar_trt + pred
       , data = dt
       , family = binomial(logit)
       , discrete = TRUE
       , nthreads = 12
         )

t <- summary(m)

pt <- data.table(
    name=names(t$p.coeff)
  , coef=t$p.coeff
  , se=t$se
  , p=t$p.pv
)

pt <- pt[, lower := coef - qnorm(0.975) * se
         ][, upper := coef + qnorm(0.975) * se
           ][, ci := ci_to_str(exp(coef), exp(lower), exp(upper), 2)]
pt[2:7, .(name, ci, p)] %>% prt() %>% add_footnote(label=Reduce(paste,(deparse(m$formula))))

dt$prob <- m_trt$fitted.values

dt <- dt[, weight_pre := case_when(
             trt %in% c(1) ~ 1 / prob
           , TRUE ~ 1 / (1 - prob)
           )][, weight := weight_pre / mean(weight_pre)]

library(mgcv)


m0 <- bam(outcome ~ trt
       , data = dt
       , weights = weight
       , family = binomial(logit)
       , discrete = TRUE
       , nthreads = 12
         )

m <- bam(outcome ~ severity + trt + covar1 + covar2 + covar3 + covar_trt
       , data = dt
       , weights = weight
       , family = binomial(logit)
       , discrete = TRUE
       , nthreads = 12
         )

t <- summary(m)

pt <- data.table(
    name=names(t$p.coeff)
  , coef=t$p.coeff
  , se=t$se
  , p=t$p.pv
)

pt <- pt[, lower := coef - qnorm(0.975) * se
         ][, upper := coef + qnorm(0.975) * se
           ][, ci := ci_to_str(exp(coef), exp(lower), exp(upper), 2)]
pt[2:7, .(name, ci, p)] %>% prt() %>% add_footnote(label=Reduce(paste,(deparse(m$formula))))

library(optmatch)
match <- match_on(
    trt ~ prob
  , data = dt
  , method = "mahalanobis"
  )

matched_1to1 <- pairmatch(match, controls = 1, data = dt)
dt$matched_id_11 <- ifelse(is.na(matched_1to1), NA, paste0("", matched_1to1, sep = ""))

dt2 <- dt[!is.na(matched_id_11)]

library(RItools)
bt <- xBalance(trt ~ severity + covar1 + covar2 + covar3 + covar_trt
             , data = dt2
               , strata = list(unstrat = NULL, pt = ~ matched_id_11)
               , report = "all")

t <- bt$results[,, 1]

t %>% prt(digits = 4, caption = "Balance Table")

t1 <- bt$results[,, 1][, c(1, 2, 3, 4, 7)]
t <- t1
t %>% prt(digits=4
        , caption = "Balance Table after 1:1 Matching"
        , col.names=c("Mean (0)", "Mean (1)", "Difference", "Standardized Difference", "P-value"))

dt2 <- dt2[, matched_id_11 := factor(matched_id_11)]
m <- bam(outcome ~ severity + trt + covar1 + covar2 + covar3 + covar_trt + s(matched_id_11, bs = "re")
       , data = dt2
       , family = binomial(logit)
       , discrete = TRUE
       , nthreads = 12
         )

t <- summary(m)

pt <- data.table(
    name=names(t$p.coeff)
  , coef=t$p.coeff
  , se=t$se
  , p=t$p.pv
)

pt <- pt[, lower := coef - qnorm(0.975) * se
         ][, upper := coef + qnorm(0.975) * se
           ][, ci := ci_to_str(exp(coef), exp(lower), exp(upper), 2)]
pt[2:7, .(name, ci, p)] %>% prt() %>% add_footnote(label=Reduce(paste,(deparse(m$formula))))

library(optmatch)
match <- match_on(
    trt ~ prob
  , data = dt
  , method = "mahalanobis"
  )

matched_1ton <- fullmatch(match, data = dt)
dt$matched_id_1n <- ifelse(is.na(matched_1ton), NA, paste0("", matched_1ton, sep = ""))

dt2 <- dt[!is.na(matched_id_1n)]

library(RItools)
bt <- xBalance(trt ~ severity + covar1 + covar2 + covar3 + covar_trt
             , data = dt2
               , strata = list(unstrat = NULL, pt = ~ matched_id_1n)
               , report = "all")

t <- bt$results[,, 1]

t %>% prt(digits = 4, caption = "Balance Table")

t1 <- bt$results[,, 1][, c(1, 2, 3, 4, 7)]
t <- t1
t %>% prt(digits=4
        , caption = "Balance Table after 1:n Matching"
        , col.names=c("Mean (0)", "Mean (1)", "Difference", "Standardized Difference", "P-value"))

dt2 <- dt2[, matched_id_1n := factor(matched_id_1n)]
m <- bam(outcome ~ severity + trt + covar1 + covar2 + covar3 + covar_trt + s(matched_id_1n, bs = "re")
       , data = dt2
       , family = binomial(logit)
       , discrete = TRUE
       , nthreads = 12
         )

t <- summary(m)

pt <- data.table(
    name=names(t$p.coeff)
  , coef=t$p.coeff
  , se=t$se
  , p=t$p.pv
)

pt <- pt[, lower := coef - qnorm(0.975) * se
         ][, upper := coef + qnorm(0.975) * se
           ][, ci := ci_to_str(exp(coef), exp(lower), exp(upper), 2)]
pt[2:7, .(name, ci, p)] %>% prt() %>% add_footnote(label=Reduce(paste,(deparse(m$formula))))

Causal Forest

• grf runs regression forest

library(grf)
frml <- Wu::wu_formula(outcome = "", predictors = c("covar1", "covar2", "covar3", "covar_trt", "trt", "severity"))
mmx <- model.matrix.lm(frml, data = dt)

X <- subset(mmx, select=-c(trt))
X <- subset(X, select = -c(1))
Y <- dt$outcome
W <- mmx[, "trt"]

rf <- causal_forest(X, Y, W
                  ## , clusters = dt[["department"]]
                  , Y.hat = NULL
                  , W.hat = NULL
                  , num.trees = 1000
                  , honesty = TRUE
                  , tune.parameters = "none"
                    ## , tune.parameters = c("sample.fraction", "mtry")
                    ## , tune.num.trees = 4000
                  , compute.oob.predictions = TRUE
                  , seed = 123456
                    )


ate <- average_treatment_effect(rf, method="AIPW", target.sample = c("all"))
ate_treated <- average_treatment_effect(rf, method="AIPW", target.sample = c("treated"))


rst <- rbind(c("Average Effect", ate)
           , c("Average Effect on Treated", ate_treated)
             )

rst <- as.data.table(rst)
colnames(rst) <- c("name", "estimate", "se")
rst <- rst[, estimate := as.numeric(estimate)
           ][, se := as.numeric(se)
             ][, lower := estimate - (qnorm(0.975) * se)
           ][, upper := estimate + qnorm(0.975) * se]



rst %>% prt(digits=3)

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] grf_2.2.0 RItools_0.3-1 optmatch_0.10.4
[4] sjPlot_2.8.11 Wu_0.0.0.9000 flexdashboard_0.6.0 [7] lme4_1.1-30 Matrix_1.4-0 mgcv_1.8-38
[10] nlme_3.1-152 png_0.1-7 scales_1.2.0
[13] nnet_7.3-16 labelled_2.9.1 kableExtra_1.3.4
[16] plotly_4.10.0 gridExtra_2.3 ggplot2_3.3.6
[19] DT_0.24 tableone_0.13.2 magrittr_2.0.3
[22] lubridate_1.8.0 dplyr_1.0.9 plyr_1.8.7
[25] data.table_1.14.2 rmdformats_1.0.4 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] class_7.3-19 sjlabelled_1.2.0 estimability_1.4.1 parameters_0.18.2 [9] rstudioapi_0.13 proxy_0.4-27 farver_2.1.1 fansi_1.0.3
[13] mvtnorm_1.1-3 xml2_1.3.3 codetools_0.2-18 splines_4.2.0
[17] cachem_1.0.6 sjmisc_2.8.9 jsonlite_1.8.0 nloptr_2.0.3
[21] ggeffects_1.1.3 broom_1.0.0 effectsize_0.7.0.5 compiler_4.2.0
[25] httr_1.4.4 sjstats_0.18.1 emmeans_1.8.0 backports_1.4.1
[29] assertthat_0.2.1 fastmap_1.1.0 lazyeval_0.2.2 survey_4.1-1
[33] cli_3.3.0 htmltools_0.5.3 tools_4.2.0 coda_0.19-4
[37] gtable_0.3.0 glue_1.6.2 Rcpp_1.0.9 jquerylib_0.1.4
[41] vctrs_0.4.1 svglite_2.1.0 crosstalk_1.2.0 insight_0.18.2
[45] xfun_0.32 stringr_1.4.0 rvest_1.0.2 lifecycle_1.0.1
[49] klippy_0.0.0.9500 MASS_7.3-54 zoo_1.8-10 rlemon_0.2.0
[53] hms_1.1.1 sandwich_3.0-2 SparseM_1.81 yaml_2.3.5
[57] sass_0.4.2 stringi_1.7.8 highr_0.9 bayestestR_0.12.1 [61] e1071_1.7-11 boot_1.3-28 rlang_1.0.4 pkgconfig_2.0.3
[65] systemfonts_1.0.4 evaluate_0.16 lattice_0.20-45 purrr_0.3.4
[69] htmlwidgets_1.5.4 labeling_0.4.2 tidyselect_1.1.2 bookdown_0.28
[73] R6_2.5.1 generics_0.1.3 multcomp_1.4-20 DBI_1.1.3
[77] pillar_1.8.1 haven_2.5.0 withr_2.5.0 abind_1.4-5
[81] survival_3.2-13 datawizard_0.5.0 tibble_3.1.8 performance_0.9.2 [85] modelr_0.1.8 utf8_1.2.2 rmarkdown_2.14 grid_4.2.0
[89] forcats_0.5.1 svd_0.5.2 digest_0.6.29 webshot_0.5.3
[93] xtable_1.8-4 tidyr_1.2.0 munsell_0.5.0 viridisLite_0.4.0 [97] bslib_0.4.0 mitools_2.4