Introduction

This R Markdown file demonstrates how to estimate the Survivor Average Causal Effect (SACE) and the Conditional Survivor Average Causal Effects (CSACEs) using Bayesian Additive Regression Trees (BART). The method applies to scenarios where a binary intermediate variable truncates the continuous outcome under the Rubin potential outcome framework [@Rubin1974]. Specifically, let \(Y(a)\in\mathbb{R}\) and \(D(a)\in\{0,1\}\) respectively denote the potential final and intermediate outcomes under treatment \(a\in\{0,1\}\). In the setting of interest, the final outcome is unobserved, i.e., \(Y(a)=\ast\), when \(D(a)=1\). An example is the truncation by death, where \(D(a)\) is the survival status of the patient, and the final outcome will be missing if the patient dies before it is recorded.

This missing not at random pattern makes it difficult to draw causal inference, because certain patients only survive under treatment or control, creating a population imbalance for causal comparisons. Therefore, we adopt the principal stratification framework [@Frangakis2002biostats] to identify patients who will survive under both treatment and control, because their potential outcomes can be observed under either conditions (treatment or control). We use \(S=D(1)D(0)\) to denote the principal stratum. For example, those who survive under treatment but die under control belong to the stratum \(S=10\), since \(D(1)=1\) and \(D(0)=0\) for these patients. In theory, there are four strata in total, i.e., \(\{10,01,00,11\}\). However, the strata can only be identified under certain assumptions, which are documented in detail in @Chen2024. One of the assumptions is monotonicity, where it is assumed that the treatment is not harming the patients, i.e., \(D(1)\geq D(0)\), thus leading to the existence of three principal strata, \(\{10,00,11\}\). The estimand of interest, i.e., the SACE, is defined as \[\Delta=\mathbb{E}\{Y(1)-Y(0)|S=11\}.\] The principal stratum \(S=11\) is called the always-survivors. The estimand that captures the heterogeneity of the treatment effect due to different patient characteristics is the CSACE, \[\Delta(X)=\mathbb{E}\{Y(1)-Y(0)|S=11,X\},\] where \(X\) is the vector of baseline covariates profiling a subgroup of patients.

This R markdown file provides the workflow for implementing the approach in @Chen2024, which includes two main steps:

Estimation of SACE and CSACE using BART.
“Fit-the-fit” approach with the classification and regression trees (CART) to identify covariate-defined subgroups’ differential causal effects.

We will provide a demonstration using the Acute Respiratory Distress Syndrome Network (ARDSNet) ARDS respiratory management (ARMA) trial data [@ARDS2000], where the outcome of interest is the days to return home (DTRH), which is truncated by the intermediate variable of death. That is, the outcome of interest is missing for patients who die before the DTRH is recorded. The target estimand is the average causal effect of switching from the standard (control) condition, high tidal volume ventilation (12 mL/kg), to the new (treatment) condition, low tidal volume ventilation (6 mL/kg).

Functions

Step 1: Estimating the SACE and CSACE using BART

Function: `sace_bart`

This function estimates the SACE and CSACE using a combined Gibbs sampler. At each iteration, we update the outcome and strata models using BART, and sample latent stratum memberships via data augmentation.

Parameters

Y:
Type: Vector
Description: A numeric vector representing the outcome for each observation.
D:
Type: Vector (binary)
Description: A binary treatment assignment vector. Values should be either TRUE (treated) or FALSE (control). The vector length should match Y.
G:
Type: Vector (categorical)
Description: A categorical treatment assignment vector. Values should be either 0 (never-survivor), 1 (protected), or 2 (always-survivor). The vector length should match Y.
id:
Type: Vector
Description: A numeric vector representing the unique for each observation. The vector length should match Y.
X:
Type: Matrix
Description: A matrix of covariates for the model. Each row represents an observation, and each column is a covariate. The number of rows should match the length of Y.
seed:
Type: Integer (optional)
Default: 12345
Description: A random seed to ensure reproducibility of the estimation results.
n_trees:
Type: Integer
Description: The number of trees to be used in the BART model.
t_var:
Type: Numerical
Description: The prior variance of the terminal nodes.
n_burn:
Type: Integer Description: The number of burn-in iterations to be discarded.
n_iter:
Type: Integer Description: The number of MCMC samples kept for posterior analyses.

Returns

The function returns a list with the following elements:

SIGMA:
Type: Vector
Description: A vector containing the posterior samples of the variance parameter in the outcome model.
SACE:
Type: Vector
Description: A vector containing the posterior samples of the SACE.
CSACE:
Type: Matrix
Description: A matrix containing the posterior samples of the CSACEs.
GV:
Type: Matrix
Description: A matrix containing the posterior samples of the values of the principal stratum indicators.
GF:
Type: Matrix
Description: A matrix containing the posterior samples of the fractions of the principal strata.
M.oc.11.1:
Type: Matrix
Description: A matrix containing the posterior samples of the outcome model means of the always survivors under treatment.
M.oc.11.0:
Type: Matrix
Description: A matrix containing the posterior samples of the outcome model means of the always survivors under control.
M.oc.10.1:
Type: Matrix
Description: A matrix containing the posterior samples of the outcome model means of the protected under treatment.
m.oc.11.1.sampler:
Type: dbarts sampler
Description: The dbarts sampler of the outcome model of the always survivors under treatment.
m.oc.11.0.sampler:
Type: dbarts sampler
Description: The dbarts sampler of the outcome model of the always survivors under control.

sace_bart <- function(Y, D, G, id, X, seed = 12345, n_trees, t_var, n_burn, n_iter) {
  
  # model information
  N <- nrow(X)
  Nst <- length(which(!is.na(Y)))
  df_tmp <- cbind(G = G, id = id, X)
  S <- n_burn + n_iter
  
  set.seed(seed)
  ## initial values
  # outcome models
  library(BayesTree, quietly = T)
  m.oc.11.1 <- rep(NA, N)
  # initial estimation using BART
  m.oc.fit <- bart(x.train = as.matrix(X[!is.na(Y),]), y.train = Y[!is.na(Y)], 
                   x.test = as.matrix(X[is.na(Y),]), ndpost = 200, verbose = F)
  m.oc.11.1[!is.na(Y)] <- m.oc.fit$yhat.train.mean
  m.oc.11.1[is.na(Y)] <- m.oc.fit$yhat.test.mean
  m.oc.11.0 <- m.oc.11.1
  m.oc.10.1 <- m.oc.11.1
  sigma2 <- m.oc.fit$sigest
  
  # strata models
  library(geepack, quietly = T)
  # initial estimation using generalized linear models
  lmm1 <- geeglm((G == 0) ~ . - id - 1, family = binomial, id = id, data = df_tmp)
  m.Z <- lmm1$fitted.values[, 1]
  
  lmm2 <- geeglm((G == 1) ~ . - id - 1, family = binomial, id = id, data = df_tmp[G != 0,])
  m.W <- (as.matrix(X) %*% lmm2$coefficients)[, 1]
  
  ## priors 
  # outcome model
  cc <- 0.001  # outcome variance price terms
  dd <- 0.001
  
  # initial values of latent variables
  library(truncnorm, quietly = T)
  # every individual has a Z
  Z <- ifelse (G < 1, rtruncnorm(1, a = 0, mean = 0, sd = sqrt(1)), 
               rtruncnorm(1, b = 0, mean = 0, sd = sqrt(1)))
  # every individual in strata 10 & 11 has a W
  W <- rep(NA, N)
  W[G == 1] <- rtruncnorm(sum(G == 1), a = 0, mean = 0, sd = sqrt(1))
  W[G == 2] <- rtruncnorm(sum(G == 2), b = 0, mean = 0, sd = sqrt(1))
  W[G == 0] <- NA
  
  ## store the SACE and CSACEs at each iteration of the Gibbs sampler
  SIGMA <- rep(0, S) # outcome variance
  SACE <- rep(0, S) # SACE
  CSACE <- array(NA, dim = c(N, S)) # CSACEs
  
  # store the outcome model means at each iteration
  # BART outcome model means for G = 2 & D = 1
  M.oc.11.1 <- matrix(NA, nrow = N, ncol = S)
  # BART outcome model means for G = 2 & D = 0
  M.oc.11.0 <- matrix(NA, nrow = N, ncol = S)
  # BART outcome model means for G = 1 & D = 1
  M.oc.10.1 <- matrix(NA, nrow = N, ncol = S)
  
  # store fractions of principal strata at each iteration
  GF <- matrix(NA, nrow = 3, ncol = S)
  # store principal stratum indicators at each iteration
  GV <- matrix(NA, nrow = N, ncol = S)
  
  # store the strata model means at each iteration
  # BART strata model means for Z
  M.Z <- matrix(NA, nrow = N, ncol = S)
  # BART strata model means for W
  M.W <- matrix(NA, nrow = N, ncol = S) 
  
  # specify dparts controls for the Gibbs sampler
  library(dbarts, quietly = T)
  db.control.1 <- dbartsControl(updateState = F, verbose = F, n.burn = 0L, 
                                n.samples = 1L, n.trees = n_trees, n.thin = 1L, 
                                n.chains = 1L)
  db.control.2 <- dbartsControl(updateState = F, verbose = F, n.burn = 0L, 
                                n.samples = 1L, n.trees = n_trees, n.thin = 1L, 
                                n.chains = 1L)
  
  # start the Gibbs sampler
  for (i in 1:S) {
    # update BART outcome model means, m.oc.11.1, m.oc.11.0, m.oc.10.1
    # use df[(G == g) & (D == d),] to train, and df[!((G == g) & (D == d)),] to test and predict
    # have to save all BART means, because G may change after each iteration
    if (sum((G == 2) & (D == 1)) > ncol(X) + 1) {
      m.oc.11.1.train <- cbind(Y = Y[(G == 2) & (D == 1)], X[(G == 2) & (D == 1),])
      m.oc.11.1.test <- cbind(Y = Y[!((G == 2) & (D == 1))], X[!((G == 2) & (D == 1)),])
      m.oc.11.1.sampler <- dbarts(Y ~ ., m.oc.11.1.train, m.oc.11.1.test, 
                                  node.prior = normal(t_var), control = db.control.1)
      m.oc.11.1.samples <- m.oc.11.1.sampler$run()
      m.oc.11.1 <- rep(NA, N)
      m.oc.11.1[(G == 2) & (D == 1)] <- m.oc.11.1.samples$train[, 1]
      m.oc.11.1[!((G == 2) & (D == 1))] <- m.oc.11.1.samples$test[, 1]
    }
    M.oc.11.1[, i] <- m.oc.11.1
    
    m.oc.11.0.train <- cbind(Y = Y[(G == 2) & (D == 0)], X[(G == 2) & (D == 0),])
    m.oc.11.0.test <- cbind(Y = Y[!((G == 2) & (D == 0))], X[!((G == 2) & (D == 0)),])
    m.oc.11.0.sampler <- dbarts(Y ~ ., m.oc.11.0.train, m.oc.11.0.test, 
                                node.prior = normal(t_var), control = db.control.1)
    m.oc.11.0.samples <- m.oc.11.0.sampler$run()
    m.oc.11.0 <- rep(NA, N)
    m.oc.11.0[(G == 2) & (D == 0)] <- m.oc.11.0.samples$train[, 1]
    m.oc.11.0[!((G == 2) & (D == 0))] <- m.oc.11.0.samples$test[, 1]
    M.oc.11.0[, i] <- m.oc.11.0
    
    if (sum((G == 1) & (D == 1)) > ncol(X) + 1) {
      m.oc.10.1.train <- cbind(Y = Y[(G == 1) & (D == 1)], X[(G == 1) & (D == 1),])
      m.oc.10.1.test <- cbind(Y = Y[!((G == 1) & (D == 1))], X[!((G == 1) & (D == 1)),])
      m.oc.10.1.sampler <- dbarts(Y ~ ., m.oc.10.1.train, m.oc.10.1.test, 
                                  node.prior = normal(t_var), control = db.control.1)
      m.oc.10.1.samples <- m.oc.10.1.sampler$run()
      m.oc.10.1 <- rep(NA, N)
      m.oc.10.1[(G == 1) & (D == 1)] <- m.oc.10.1.samples$train[, 1]
      m.oc.10.1[!((G == 1) & (D == 1))] <- m.oc.10.1.samples$test[, 1]
    }
    M.oc.10.1[, i] <- m.oc.10.1
    
    # update sigma2, the variance in the outcome model
    rate <- dd + 0.5 * (sum((Y[(G == 2) & (D == 1)] - m.oc.11.1[(G == 2) & (D == 1)])**2)) + 
      0.5 * (sum((Y[(G == 1) & (D == 1)] - m.oc.10.1[(G == 1) & (D == 1)])**2)) +
      0.5 * (sum((Y[(G == 2) & (D == 0)] - m.oc.11.0[(G == 2) & (D == 0)])**2))
    sigma2 <- rgamma(1, shape = (cc + Nst/2), rate = rate)**(-1)
    SIGMA[i] <- sigma2
    
    # calculate the SACE and the CSACEs at each iteration
    # the mean potential outcomes under treatment for the always-survivor group
    Y1 <- m.oc.11.1[G == 2]  
    # the mean potential outcomes under control for the always-survivor group
    Y0 <- m.oc.11.0[G == 2]
    # calculate the SACE
    sace <- mean(Y1 - Y0)  
    SACE[i] <- sace
    # calculate the CSACEs
    CSACE[, i] <- m.oc.11.1 - m.oc.11.0
    
    # update BART strata model means, m.Z, m.W
    m.Z.train <- cbind(Z = Z, X)
    m.Z.sampler <- dbarts(Z ~ ., m.Z.train, node.prior = normal(t_var), 
                          control = db.control.2)
    m.Z.samples <- m.Z.sampler$run()
    m.Z <- m.Z.samples$train[, 1]
    M.Z[, i] <- m.Z

    m.W.train <- cbind(W = W[G != 0], X[G != 0,])
    m.W.test <- X[G == 0,]
    m.W.sampler <- dbarts(W ~ ., m.W.train, m.W.test, node.prior = normal(t_var), 
                          control = db.control.2)
    m.W.samples <- m.W.sampler$run()
    m.W = rep(NA, N)
    m.W[G != 0] <- m.W.samples$train[, 1]
    m.W[G == 0] <- m.W.samples$test[, 1]
    M.W[, i] <- m.W

    # update principal stratum indicators G via data augmentation
    G <- rep(NA, N)
    # D = 1, S(1) = 0, identified never survivors
    G[is.na(Y) & (D == 1)] <- 0
    # D = 0, S(0) = 1, identified always survivors
    G[(!is.na(Y)) & (D == 0)] <- 2 
    # D = 0, S(0) = 0, protected or never survivors
    tgt.D0S0 <- intersect(which(is.na(Y)), which(D == 0))
    p00.Z <- pnorm(m.Z[tgt.D0S0])
    p10.Z <- (1 - pnorm(m.Z[tgt.D0S0])) * pnorm(m.W[tgt.D0S0])
    G[tgt.D0S0] <- rep(1, length(tgt.D0S0)) - rbinom(length(tgt.D0S0), 1, p00.Z / (p00.Z + p10.Z))
    # D = 1, S(1) = 1, protected or always survivors
    tgt.D1S1 <- intersect(which(!is.na(Y)), which(D == 1))
    p10.W <- pnorm(m.W[tgt.D1S1]) * dnorm(Y[tgt.D1S1], m.oc.10.1[tgt.D1S1], sqrt(sigma2))
    p11.W <- (1 - pnorm(m.W[tgt.D1S1])) * dnorm(Y[tgt.D1S1], m.oc.11.1[tgt.D1S1], sqrt(sigma2))
    G[tgt.D1S1] <- rep(2, length(tgt.D1S1)) - rbinom(length(tgt.D1S1), 1, p10.W / (p10.W + p11.W))
    GV[, i] <- G
    GF[, i] <- c(sum(G == 0), sum(G == 1), sum(G == 2)) / N

    # update Z & W
    Z <- rep(NA, N)
    Z[G == 0] <- rtruncnorm(sum(G == 0), a = 0, mean = (m.Z[G == 0]), sd = 1)
    Z[G != 0] <- rtruncnorm(sum(G != 0), b = 0, mean = (m.Z[G != 0]), sd = 1)
    W <- rep(NA, N)
    if (sum(G == 1) > 0) {
      W[G == 1] <- rtruncnorm(sum(G == 1), a = 0, mean = (m.W[G == 1]), sd = 1)
    }
    if (sum(G == 2) > 0) {
      W[G == 2] <- rtruncnorm(sum(G == 2), b = 0, mean = (m.W[G == 2]), sd = 1)
    }
    W[G == 0] <- NA

    if (i %% 100 == 0) {
      cat("\r", paste("sampling iteration", i, sep=" "))
      flush.console()
    }
  }   
  return(list(SIGMA = SIGMA, SACE = SACE, CSACE = CSACE, GV = GV, GF = GF,
              M.oc.11.1 = M.oc.11.1, M.oc.11.0 = M.oc.11.0, M.oc.10.1 = M.oc.10.1,
              M.Z = M.Z, M.W = M.W, m.oc.11.1.sampler = m.oc.11.1.sampler, 
              m.oc.11.0.sampler = m.oc.11.0.sampler))
}

Step 2: “Fit-the-Fit” Approach with the Classification and Regression Trees (CART) Model for CSACE Subgroup Identification

In the second stage, the “fit-the-fit” approach is used to find covariate-defined subgroups exhibiting heterogeneity of causal effect. In particular, a CART model is fit with the estimated CSACE values from the first stage as the outcome and the covariates as possible predictors. The model is first fit under default CART hyperparameter settings.

# Load necessary libraries and source the custom functions
library(caret)
library(rpart)
library(rpart.plot)
library(flextable)
library(officer)
library(ggh4x)

Function: `cart_fit`

This function fits a CART model on the estimated CSACE to identify differential treatment effects across covariate-defined subgroups.

Parameters

csace.est:
Type: Numeric vector
Description: The estimated conditional survivor average causal effect (CSACE) values for each observation, obtained from the first-stage BART model.
X:
Type: Matrix or Data Frame
Description: A matrix or data frame of covariates to be used as predictors in the CART model.
maxdepth:
Type: Integer (optional)
Default: 3
Description: Maximum depth of the CART tree. A lower value prunes the tree for better interpretability.

Returns

This function produces a fitted CART model that identifies subgroups with different CSACE values and plots the resulting decision tree. Additionally, for each subgroup identified by the CART model, the function prints the following:

Mean CSACE: The average CSACE for the subgroup.
Lower 2.5% Quantile: The lower bound of the 95% credible interval for the subgroup CSACE.
Upper 97.5% Quantile: The upper bound of the 95% credible interval for the subgroup CSACE.

cart_fit <- function(csace.est, X, maxdepth = 3){
  
  # Combine covariates and CSACE estimates into one data frame
  data <- data.frame(X, csace.est = csace.est)
  
  # Fit the CART model using rpart, limiting the tree depth for interpretability
  cartmod <- rpart(csace.est ~ ., data = data, method = "anova", maxdepth = maxdepth)
  
  # Plot the CART model
  rpart.plot(cartmod, yesno = 2)
  
  # Loop through each unique terminal node in the fitted tree
  terminal_nodes <- unique(cartmod$where)  # Unique terminal node IDs
  for (i in terminal_nodes) {
    # Get the indices of observations in each node
    node_indices <- which(cartmod$where == i)
    
    # Calculate mean and credible intervals (2.5% and 97.5%) for CSACE in the node
    node_mean <- mean(data$csace.est[node_indices])
    node_lower <- quantile(data$csace.est[node_indices], 0.025)
    node_upper <- quantile(data$csace.est[node_indices], 0.975)
    
    # Print results for the current node
    cat("\nNode:", i, "\n")
    cat("Mean CSACE:", node_mean, "\n")
    cat("95% CI:", node_lower, "-", node_upper, "\n")
  }
}

ARDS ARMA data analysis

First, we load the ARMA data.

df <- readRDS("ARMA_ARDS.RDS")

Then, we format the data for the ensuing analysis. \(d\) here represents the treatment indicator variable, which equals 1 for the treatment (6 mL/kg) condition and 0 for the control (12 mL/kg) condition. \(G\) is the principal stratum indicator, which takes values of 0, 1, and 2. There are three principal strata because of the monotonicity assumption, such that the new condition does not harm the health of the patients. Under this assumption, the entire patient population can be divided into three subpopulations: 1. the always survivors (\(G=2\)), who would survive under either the treatment or control condition, 2. the protected (\(G=1\)), who would survive under treatment but die under control, and 3. the never survivors (\(G=0\)), who would die under either the treatment or control condition.

We can identify the principal strata for certain patients because of the monotonicity assumption. This includes the following two groups: 1. those who survived under control (always survivors), and 2. those who died under treatment (never survivors). For the other two combinations of treatment and survival statuses, the principal stratum memberships cannot be directly determined. Specifically, those who died under control could be never survivors or protected, and those who survived under treatment could be always survivors or protected. Therefore, we randomly assign the principal stratum memberships to those patients as initial values.

df$d <- 1*(df$trt1 == "Randomized: 6 ml/kg")
df$G <- NA
for (i in 1:nrow(df)) {
  # control but survived
  if (df$d[i] == 0 & df$day_death[i] == 0) {df$G[i] <- 2} 
  # treatment but died
  if (df$d[i] == 1 & df$day_discharge[i] == 0) {df$G[i] <- 0} 
  # control and died
  if (df$d[i] == 0 & df$day_discharge[i] == 0) {df$G[i] <- rbinom(1, 1, 0.5)} 
  # treatment and survived
  if (df$d[i] == 1 & df$day_death[i] == 0) {df$G[i] <- 1 + rbinom(1, 1, 0.5)} 
}

We standardize continuous covariates and convert categorical covariates into binary indicator variables.

df$day_discharge[df$day_discharge == 0] <- NA
df$day_death[df$day_death == 0] <- NA
std <- function(x) {(x - mean(x)) / sd(x)}
df2 <- data.frame(id = 1:nrow(df), Y = df$day_discharge, d = df$d, G = df$G, 
                  TIDAL = std(df$TIDAL), PEEP = std(df$PEEP), FIO2 = std(df$FIO2), 
                  PACO2 = std(df$PACO2), PAO2 = std(df$PAO2), ARTPH = std(df$ARTPH), 
                  AGE = std(df$AGE), APACHE = std(df$APACHE), glasgow = std(df$glasgow), 
                  SYSBP = std(df$SYSBP), PAFI = std(df$PAFI), PLATE = std(df$PLATE), 
                  CREAT = std(df$CREAT), BILI = std(df$BILI), aado2 = std(df$aado2), 
                  ptof = std(df$ptof), GENDER = df$GENDER - 1, VASO = df$VASO - 1, 
                  Non_White = 1*(df$ETHNIC != 1), PNEUM1 = 1*(df$PNEUM == 1), 
                  PNEUM2 = 1*(df$PNEUM == 2), SEPSIS1 = 1*(df$SEPSIS == 1), 
                  SEPSIS2 = 1*(df$SEPSIS == 2), ASPIR1 = 1*(df$ASPIR == 1), 
                  ASPIR2 = 1*(df$ASPIR == 2), TRAUMA1 = 1*(df$TRAUMA == 1), 
                  TRAUMA2 = 1*(df$TRAUMA == 2), OTHER1 = 1*(df$OTHER == 1), 
                  OTHER2 = 1*(df$OTHER == 2), MULTIRAN1 = 1*(df$MULTRAN == 1), 
                  MULTIRAN2 = 1*(df$MULTRAN == 2))

We run the combined Gibbs sampler.

Y <- df2$Y
D <- df2$d
G <- df2$G
id <- df2$id
X <- df2[,c('TIDAL', 'PEEP', 'FIO2', 'PACO2', 'PAO2', 'ARTPH', 'AGE',
            'GENDER', 'Non_White', 'APACHE', 'glasgow', 'SYSBP', 'PLATE',
            'CREAT', 'BILI', 'aado2', 'ptof', 'VASO')]
n_trees <- 50 
t_var <- 4 
n_burn <- 5000
n_iter <- 5000

Res <- sace_bart(Y, D, G, id, X, seed = 12345, n_trees, t_var, n_burn, n_iter)

##  sampling iteration 100 sampling iteration 200 sampling iteration 300 sampling iteration 400 sampling iteration 500 sampling iteration 600 sampling iteration 700 sampling iteration 800 sampling iteration 900 sampling iteration 1000 sampling iteration 1100 sampling iteration 1200 sampling iteration 1300 sampling iteration 1400 sampling iteration 1500 sampling iteration 1600 sampling iteration 1700 sampling iteration 1800 sampling iteration 1900 sampling iteration 2000 sampling iteration 2100 sampling iteration 2200 sampling iteration 2300 sampling iteration 2400 sampling iteration 2500 sampling iteration 2600 sampling iteration 2700 sampling iteration 2800 sampling iteration 2900 sampling iteration 3000 sampling iteration 3100 sampling iteration 3200 sampling iteration 3300 sampling iteration 3400 sampling iteration 3500 sampling iteration 3600 sampling iteration 3700 sampling iteration 3800 sampling iteration 3900 sampling iteration 4000 sampling iteration 4100 sampling iteration 4200 sampling iteration 4300 sampling iteration 4400 sampling iteration 4500 sampling iteration 4600 sampling iteration 4700 sampling iteration 4800 sampling iteration 4900 sampling iteration 5000 sampling iteration 5100 sampling iteration 5200 sampling iteration 5300 sampling iteration 5400 sampling iteration 5500 sampling iteration 5600 sampling iteration 5700 sampling iteration 5800 sampling iteration 5900 sampling iteration 6000 sampling iteration 6100 sampling iteration 6200 sampling iteration 6300 sampling iteration 6400 sampling iteration 6500 sampling iteration 6600 sampling iteration 6700 sampling iteration 6800 sampling iteration 6900 sampling iteration 7000 sampling iteration 7100 sampling iteration 7200 sampling iteration 7300 sampling iteration 7400 sampling iteration 7500 sampling iteration 7600 sampling iteration 7700 sampling iteration 7800 sampling iteration 7900 sampling iteration 8000 sampling iteration 8100 sampling iteration 8200 sampling iteration 8300 sampling iteration 8400 sampling iteration 8500 sampling iteration 8600 sampling iteration 8700 sampling iteration 8800 sampling iteration 8900 sampling iteration 9000 sampling iteration 9100 sampling iteration 9200 sampling iteration 9300 sampling iteration 9400 sampling iteration 9500 sampling iteration 9600 sampling iteration 9700 sampling iteration 9800 sampling iteration 9900 sampling iteration 10000

Load the required R packages for the results visualization.

library(ggplot2)
library(reshape2)
library(plyr)
library(gridExtra)
library(grid)
library(lattice)
library(Matrix)
library(hrbrthemes)

We first produce a caterpillar plot for the CSACEs.

# We first find the identified always survivors and save their indexes as S1.
# S1 has 260 patients.
S1 <- intersect(which(!is.na(df2$Y)), which(df2$d == 0))

# We then find the likely always survivors which are with posterior probabilities
# greater than 0.8 being an always-survivor.
G11_fr <- sapply(1:nrow(df2), function(x) sum(Res$GV[x, -(1:n_burn)] == 2) / n_iter)

# We save the indexes as S1Sp.
# S1Sp has 522 patients.
S1Sp <- which(G11_fr >= 0.8)

# Calculate the estimated CSACE for each patient in S1Sp.
IE_mean_S1Sp <- apply(Res$CSACE[S1Sp, -(1:n_burn)], 1, mean)
IE_S1Sp <- Res$CSACE[S1Sp, -(1:n_burn)]
# Order these CSACEs increasingly.
IE_S1Sp_ordered <- IE_S1Sp[order(IE_mean_S1Sp, decreasing = F),]

# Calculate the lower and upper bounds for the confidence bands.
ub <- apply(IE_S1Sp_ordered, 1, function(x) quantile(x, 0.975))
est <- apply(IE_S1Sp_ordered, 1, mean)
lb <- apply(IE_S1Sp_ordered, 1, function(x) quantile(x, 0.025))
dat_IE_S1Sp <- data.frame(index = 1:nrow(IE_S1Sp_ordered), est = est, ub = ub, lb = lb)

# Produce the caterpillar plot.
fig_IE_S1Sp <- ggplot(dat_IE_S1Sp, aes(x = index,y = est)) +
  geom_smooth(aes(ymin = lb, ymax = ub), stat = "identity", fill = alpha(5,0.3), colour = 4) +
  xlab("") + ylab("") + ggtitle("Estimated CSACE for Likely Always-Survivors") + theme_minimal() +
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12)) + 
  scale_x_continuous(breaks = round(seq(0, 530, by = 100), 1)) + 
  geom_hline(yintercept = 0, color = 10, size = 1.2, linetype = 2)

yleft <- textGrob("CSACE", rot = 90, gp = gpar(fontsize = 15))
bottom <- textGrob("Participants Index", gp = gpar(fontsize = 15))

grid.arrange(arrangeGrob(fig_IE_S1Sp, nrow = 1, ncol = 1), left = yleft, bottom = bottom)

We then produce mirror histograms of the SACE. Specifically, the first histogram of the SACE is produced using its posterior values. The second histogram of the SACEs is produced using the posterior mean CSACEs of all patients in S1Sp.

SACE_data <- data.frame(ind = 1:length(colMeans(IE_S1Sp)), sace = Res$SACE[-(1:n_burn)], csace = colMeans(IE_S1Sp))

mirror <- ggplot(SACE_data, aes(x = ind) ) +
  geom_histogram(aes(x = sace, y = ..density..), fill = alpha(5, 0.3), colour = 4) +
  geom_label(aes(x = -12, y = 0.15, label = "SACE 1 Posterior"), color = 4) +
  geom_histogram(aes(x = csace, y = -..density..), fill = alpha(2, 0.3), colour = 10) +
  geom_label(aes(x = -12, y = -0.15, label = "SACE 2 Posterior"), color = 10) +
  theme_ipsum() +
  xlab("") + ylab("") + 
  geom_segment(aes(x = mean(Res$SACE[-(1:n_burn)]), y = 0, xend = mean(Res$SACE[-(1:n_burn)]), yend = 0.15), 
               color = 4, size = 1.2, linetype = 2) + 
  geom_segment(aes(x = mean(colMeans(IE_S1Sp)), y = 0, xend = mean(colMeans(IE_S1Sp)), yend = -0.15), 
               color = 10, size = 1.2, linetype = 2) + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12))

yleft <- textGrob("Density", rot = 90, gp = gpar(fontsize = 15))
bottom <- textGrob("SACE", gp = gpar(fontsize = 15))

grid.arrange(arrangeGrob(mirror, nrow = 1, ncol = 1), left = yleft, bottom = bottom)

We produce histograms of the posterior means of the CSACEs and their posterior mean density.

dat_3.3 <- data.frame(x = apply(Res$CSACE[S1Sp, -(1:n_burn)], 1, mean))
dat_3.3.2 <- data.frame(x = c(Res$CSACE[S1Sp, -(1:n_burn)]))

fig_3.3.1 <- ggplot() + 
  geom_histogram(data = dat_3.3, aes(x = x, y = ..density..),
                 colour = 4, fill = alpha(5, 0.3)) +
  geom_vline(xintercept = 0, color = 10, size = 1.2, linetype = 2) +
  geom_hline(yintercept = 0, color = 4, size = 0.5, linetype = 1) +
  theme_minimal() +
  xlab("") + ylab("") + ggtitle("Histogram of Posterior Means of CSACE") + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12)) + 
  scale_x_continuous(breaks = round(seq(-60, 20, by = 20), 1), limits = c(-80, 30))

fig_3.3.2 <- ggplot() + 
  geom_density(data = dat_3.3.2, aes(x = x),
               lwd = 1.2, linetype = 1, colour = 4) + 
  geom_vline(xintercept = 0, color = 10, size = 1.2, linetype = 2) +
  geom_hline(yintercept = 0, color = 4, size = 0.5, linetype = 1) +
  theme_minimal() +
  xlab("") + ylab("") + ggtitle("Posterior Mean Density of CSACE") + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size=12)) + 
  scale_x_continuous(breaks = round(seq(-60, 20, by = 20), 1), limits = c(-80, 30))

grid.arrange(fig_3.3.1, fig_3.3.2, nrow = 1)

We perform the “fit-the-fit” analysis.

# Form the dataset by combining the CSACE and selecting covariates 
dat_cart <- cbind(apply(IE_S1Sp, 1, mean), 
                  df[S1Sp, names(df2) %in% c("TIDAL", "PEEP", "FIO2", "PACO2", 
                                             "PAO2", "ARTPH", "AGE", "APACHE",
                                             "glasgow", "SYSBP", "PLATE", "CREAT",
                                             "BILI", "aado2", "ptof", "GENDER", 
                                             "VASO", "Non_White")])
names(dat_cart)[1] = "CSACE"

cart_fit(dat_cart[,1], dat_cart[,-1], maxdepth = 3)

## 
## Node: 14 
## Mean CSACE: -19.44485 
## 95% CI: -27.9618 - -12.48688 
## 
## Node: 11 
## Mean CSACE: -25.39102 
## 95% CI: -31.83523 - -19.21139 
## 
## Node: 15 
## Mean CSACE: -14.74963 
## 95% CI: -21.99009 - -8.656816 
## 
## Node: 12 
## Mean CSACE: -21.607 
## 95% CI: -27.42787 - -16.70351 
## 
## Node: 5 
## Mean CSACE: -31.01948 
## 95% CI: -37.59328 - -24.07658 
## 
## Node: 4 
## Mean CSACE: -36.23433 
## 95% CI: -42.89518 - -29.85862 
## 
## Node: 7 
## Mean CSACE: -26.73198 
## 95% CI: -31.32094 - -20.67475 
## 
## Node: 8 
## Mean CSACE: -21.55447 
## 95% CI: -27.51481 - -17.10354

Additional plots

We create additional plots to better visualize the CSACEs. Specifically, we categorize the CSACEs into the lower, middle, and upper thirds and examine the effects of various covariates on these causal effects [@Buell2024]. First, we create the histogram of the means of the CSACEs.

# Categorize the CSACEs into the lower, middle, and upper thirds.
est <- IE_mean_S1Sp
qt_grp <- rep(0, length(est))
for (i in 1:length(est)) {
  if (est[i] < quantile(est, 1/3)) {
    qt_grp[i] <- "Lower third"
  } else {
    if (est[i] < quantile(est, 2/3)) {
      qt_grp[i] <- "Middle third"
    } else {
      qt_grp[i] <- "Upper third"
    }
  }
}

dat_IE_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = est, 
                         qt_grp = as.factor(qt_grp), predictor = "ITE")

hist_S1Sp <- ggplot() + 
  geom_histogram(data = dat_IE_S1Sp, aes(x = est, fill = qt_grp, colour = qt_grp), 
                 alpha = 1, position = "identity") +
  theme_bw() +
  xlab("") + ylab("") + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12)) + 
  scale_fill_manual(values = c("red3", "grey80", "blue3")) + 
  scale_color_manual(values = c("red3", "grey80", "blue3")) + 
  guides(fill = guide_legend(title = ""), colour = guide_legend(title = ""))

yleft <- textGrob("No. of patients", rot = 90, gp = gpar(fontsize = 15))
bottom <- textGrob("CSACE", gp = gpar(fontsize = 15))

grid.arrange(arrangeGrob(hist_S1Sp, nrow = 1, ncol = 1), left = yleft, bottom = bottom)

Next, we create the plot visualizing the effects of different covariates. This is conducted by estimating the counterfactual CSACEs when replacing the specific covariate values by its median.

df2_S1Sp <- df2[S1Sp,]

m111.sampler <- Res$m.oc.11.1.sampler
m110.sampler <- Res$m.oc.11.0.sampler

# The following covariates are considered
# GLASGOW, PACO2, BILI, TIDAL, SYSBP, AGE, ARTPH, PtoF, PLATE, VASO, AADO2, APACHE, FIO2
# Create a function to generate counterfactual CSACEs
counter_gen <- function(df, sampler1, sampler0, variable) {
  df_1 <- df
  df_1[,variable] <- median(df_1[,variable])
  
  sampler1$setTestPredictor(df_1)
  sampler0$setTestPredictor(df_1)
  
  new_samples_1 <- sampler1$run(numBurnIn = 0, numSamples = 100, updateState = FALSE)
  new_samples_0 <- sampler0$run(numBurnIn = 0, numSamples = 100, updateState = FALSE)
  
  predictions_1 <- new_samples_1$test
  mean_preds_1 <- rowMeans(predictions_1)
  
  predictions_0 <- new_samples_0$test
  mean_preds_0 <- rowMeans(predictions_0)
  
  counter <- mean_preds_1 - mean_preds_0
  
  return(counter)
}

dat_SACE_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = mean(IE_mean_S1Sp), 
                           qt_grp = as.factor(qt_grp), predictor = "SACE")

## glasgow
counter_glasgow <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'glasgow')
dat_glasgow_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_glasgow, 
                               qt_grp = as.factor(qt_grp), predictor = "GLASGOW")

## paco2
counter_paco2 <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'PACO2')
dat_paco2_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_paco2, 
                             qt_grp = as.factor(qt_grp), predictor = "PACO2")

## bili
counter_bili <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'BILI')
dat_bili_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_bili, 
                            qt_grp = as.factor(qt_grp), predictor = "BILI")

## tidal
counter_tidal <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'TIDAL')
dat_tidal_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_tidal, 
                             qt_grp = as.factor(qt_grp), predictor = "TIDAL")

## sysbp
counter_sysbp <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'SYSBP')
dat_sysbp_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_sysbp, 
                             qt_grp = as.factor(qt_grp), predictor = "SYSBP")

## age
counter_age <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'AGE')
dat_age_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_age, 
                           qt_grp = as.factor(qt_grp), predictor = "AGE")

## artph
counter_artph <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'ARTPH')
dat_artph_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_artph, 
                             qt_grp = as.factor(qt_grp), predictor = "ARTPH")

## ptof
counter_ptof <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'ptof')
dat_ptof_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_ptof, 
                            qt_grp = as.factor(qt_grp), predictor = "PTOF")

## plate
counter_plate <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'PLATE')
dat_plate_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_plate, 
                             qt_grp = as.factor(qt_grp), predictor = "PLATE")

## vaso
counter_vaso <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'VASO')
dat_vaso_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_vaso, 
                            qt_grp = as.factor(qt_grp), predictor = "VASO")

## aado2
counter_aado2 <- counter_gen(df2_S1Sp, m111.sampler, m110.sampler, 'aado2')
dat_aado2_S1Sp <- data.frame(index = 1:length(IE_mean_S1Sp), est = counter_aado2, 
                             qt_grp = as.factor(qt_grp), predictor = "AADO2")

dat_figure_1a = rbind(dat_SACE_S1Sp, dat_glasgow_S1Sp, dat_paco2_S1Sp, dat_bili_S1Sp, 
                      dat_tidal_S1Sp, dat_sysbp_S1Sp, dat_age_S1Sp, dat_artph_S1Sp, 
                      dat_ptof_S1Sp, dat_plate_S1Sp, dat_vaso_S1Sp, dat_aado2_S1Sp)
dat_figure_1a$id = rep(1:length(IE_mean_S1Sp), 12)

dat_figure_1a$predictor <- 
  factor(dat_figure_1a$predictor, levels = c("SACE", "GLASGOW", "PACO2", "BILI", 
                                             "TIDAL", "SYSBP", "AGE", "ARTPH", 
                                             "PTOF", "PLATE", "VASO", "AADO2"))

lines_S1Sp <- ggplot(data = dat_figure_1a, aes(x = predictor, y = est, group = id)) + 
  geom_line(aes(color = qt_grp)) +
  geom_point(aes(color = qt_grp)) +
  theme_bw() +
  xlab("") + ylab("") + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12)) + 
  scale_color_manual(values = c("red3", "grey80", "blue3")) + 
  guides(fill = guide_legend(title = ""), colour = guide_legend(title = "")) + 
  scale_x_discrete(limits = c("SACE", "GLASGOW", "PACO2", "BILI", "TIDAL", "SYSBP", 
                              "AGE", "ARTPH", "PTOF", "PLATE", "VASO", "AADO2")) 

yleft = textGrob("CSACEs", rot = 90, gp = gpar(fontsize = 15))
bottom = textGrob("Baseline covariates", gp = gpar(fontsize = 15))

grid.arrange(arrangeGrob(lines_S1Sp, nrow = 1, ncol = 1), left = yleft, bottom = bottom)

Last, we plot the covariate effects on three example individuals from the lower, middle, and upper thirds.

upper_third_example = unique(dat_figure_1a$id[dat_figure_1a$qt_grp == "Upper third"])[1]
middle_third_example = unique(dat_figure_1a$id[dat_figure_1a$qt_grp == "Middle third"])[1]
lower_third_example = unique(dat_figure_1a$id[dat_figure_1a$qt_grp == "Lower third"])[1]

dat_efigure_3_u = dat_figure_1a[dat_figure_1a$id == upper_third_example,]
dat_efigure_3_u_percent = (dat_efigure_3_u$est[-nrow(dat_efigure_3_u)] - dat_efigure_3_u$est[-1])/dat_efigure_3_u$est[-nrow(dat_efigure_3_u)]

lines_efigure_S1Sp_u = ggplot(data = dat_efigure_3_u[-1.], aes(x = predictor, y = est)) + 
  geom_point(color = "grey80", size = 1) +
  geom_point(x = "GLASGOW", y = dat_efigure_3_u$est[1], shape = 21, size = 5) +
  geom_point(x = "AADO2", y = dat_efigure_3_u$est[12], shape = 4, size = 5) +
  theme_bw() +
  scale_x_discrete(limits = c("GLASGOW", "PACO2", "BILI", "TIDAL", "SYSBP", "AGE", "ARTPH", "PTOF", "PLATE", "VASO", "AADO2"))  + 
  xlab("") + ylab("") + # ggtitle("ISCE Sp") + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12)) + 
  geom_segment(x = "GLASGOW", xend = "PACO2", y = dat_efigure_3_u$est[2], yend = dat_efigure_3_u$est[2], colour = "grey80", size = 1) +
  geom_segment(x = "PACO2", xend = "BILI", y = dat_efigure_3_u$est[3], yend = dat_efigure_3_u$est[3], colour = "grey80", size = 1) +
  geom_segment(x = "BILI", xend = "TIDAL", y = dat_efigure_3_u$est[4], yend = dat_efigure_3_u$est[4], colour = "grey80", size = 1) +
  geom_segment(x = "TIDAL", xend = "SYSBP", y = dat_efigure_3_u$est[5], yend = dat_efigure_3_u$est[5], colour = "grey80", size = 1) +
  geom_segment(x = "SYSBP", xend = "AGE", y = dat_efigure_3_u$est[6], yend = dat_efigure_3_u$est[6], colour = "grey80", size = 1) +
  geom_segment(x = "AGE", xend = "ARTPH", y = dat_efigure_3_u$est[7], yend = dat_efigure_3_u$est[7], colour = "grey80", size = 1) +
  geom_segment(x = "ARTPH", xend = "PTOF", y = dat_efigure_3_u$est[8], yend = dat_efigure_3_u$est[8], colour = "grey80", size = 1) +
  geom_segment(x = "PTOF", xend = "PLATE", y = dat_efigure_3_u$est[9], yend = dat_efigure_3_u$est[9], colour = "grey80", size = 1) +
  geom_segment(x = "PLATE", xend = "VASO", y = dat_efigure_3_u$est[10], yend = dat_efigure_3_u$est[10], colour = "grey80", size = 1) +
  geom_segment(x = "VASO", xend = "AADO2", y = dat_efigure_3_u$est[11], yend = dat_efigure_3_u$est[11], colour = "grey80", size = 1) + 
  geom_segment(x = "GLASGOW", y = dat_efigure_3_u$est[1], xend = "GLASGOW", yend = dat_efigure_3_u$est[2], lineend = "round", linejoin = "round",
    size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "PACO2", y = dat_efigure_3_u$est[2], xend = "PACO2", yend = dat_efigure_3_u$est[3], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "BILI", y = dat_efigure_3_u$est[3], xend = "BILI", yend = dat_efigure_3_u$est[4], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "TIDAL", y = dat_efigure_3_u$est[4], xend = "TIDAL", yend = dat_efigure_3_u$est[5], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "SYSBP", y = dat_efigure_3_u$est[5], xend = "SYSBP", yend = dat_efigure_3_u$est[6], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "AGE", y = dat_efigure_3_u$est[6], xend = "AGE", yend = dat_efigure_3_u$est[7], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "ARTPH", y = dat_efigure_3_u$est[7], xend = "ARTPH", yend = dat_efigure_3_u$est[8], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "PTOF", y = dat_efigure_3_u$est[8], xend = "PTOF", yend = dat_efigure_3_u$est[9], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "PLATE", y = dat_efigure_3_u$est[9], xend = "PLATE", yend = dat_efigure_3_u$est[10], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "VASO", y = dat_efigure_3_u$est[10], xend = "VASO", yend = dat_efigure_3_u$est[11], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "AADO2", y = dat_efigure_3_u$est[11], xend = "AADO2", yend = dat_efigure_3_u$est[12], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  annotate(geom = "text", x = "GLASGOW", y = dat_efigure_3_u$est[2] + 0.2, label = "+36.9%", color = "red3") + 
  annotate(geom = "text", x = "PACO2", y = dat_efigure_3_u$est[3] - 0.2, label = "-18.3%", color = "blue3") + 
  annotate(geom = "text", x = "BILI", y = dat_efigure_3_u$est[4] - 0.2, label = "-14.9%", color = "blue3") + 
  annotate(geom = "text", x = "TIDAL", y = dat_efigure_3_u$est[5] + 0.2, label = "+10.0%", color = "red3") + 
  annotate(geom = "text", x = "SYSBP", y = dat_efigure_3_u$est[6] + 0.2, label = "+27.4%", color = "red3") + 
  annotate(geom = "text", x = "AGE", y = dat_efigure_3_u$est[7] - 0.2, label = "-31.0%", color = "blue3") + 
  annotate(geom = "text", x = "ARTPH", y = dat_efigure_3_u$est[8] + 0.2, label = "+13.9%", color = "red3") + 
  annotate(geom = "text", x = "PTOF", y = dat_efigure_3_u$est[9] - 0.2, label = "-17.6%", color = "blue3") + 
  annotate(geom = "text", x = "PLATE", y = dat_efigure_3_u$est[10] + 0.2, label = "+20.2%", color = "red3") + 
  annotate(geom = "text", x = "VASO", y = dat_efigure_3_u$est[11] - 0.2, label = "-15.0%", color = "blue3") +
  annotate(geom = "text", x = "AADO2", y = dat_efigure_3_u$est[12] - 0.2, label = "-11.2%", color = "blue3") + 
  geom_text(x = -Inf, y = Inf, label = "Patient from upper third", hjust = 0, vjust = 1)


dat_efigure_3_m = dat_figure_1a[dat_figure_1a$id == middle_third_example,]
dat_efigure_3_m_percent = (dat_efigure_3_m$est[-nrow(dat_efigure_3_m)] - dat_efigure_3_m$est[-1])/dat_efigure_3_m$est[-nrow(dat_efigure_3_m)]

lines_efigure_S1Sp_m = ggplot(data = dat_efigure_3_m[-1.], aes(x = predictor, y = est)) + 
  geom_point(color = "grey80", size = 1) +
  geom_point(x = "GLASGOW", y = dat_efigure_3_m$est[1], shape = 21, size = 5) +
  geom_point(x = "AADO2", y = dat_efigure_3_m$est[12], shape = 4, size = 5) +
  theme_bw() +
  scale_x_discrete(limits = c("GLASGOW", "PACO2", "BILI", "TIDAL", "SYSBP", "AGE", "ARTPH", "PTOF", "PLATE", "VASO", "AADO2"))  + 
  xlab("") + ylab("") + # ggtitle("ISCE Sp") + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12)) + 
  geom_segment(x = "GLASGOW", xend = "PACO2", y = dat_efigure_3_m$est[2], yend = dat_efigure_3_m$est[2], colour = "grey80", size = 1) +
  geom_segment(x = "PACO2", xend = "BILI", y = dat_efigure_3_m$est[3], yend = dat_efigure_3_m$est[3], colour = "grey80", size = 1) +
  geom_segment(x = "BILI", xend = "TIDAL", y = dat_efigure_3_m$est[4], yend = dat_efigure_3_m$est[4], colour = "grey80", size = 1) +
  geom_segment(x = "TIDAL", xend = "SYSBP", y = dat_efigure_3_m$est[5], yend = dat_efigure_3_m$est[5], colour = "grey80", size = 1) +
  geom_segment(x = "SYSBP", xend = "AGE", y = dat_efigure_3_m$est[6], yend = dat_efigure_3_m$est[6], colour = "grey80", size = 1) +
  geom_segment(x = "AGE", xend = "ARTPH", y = dat_efigure_3_m$est[7], yend = dat_efigure_3_m$est[7], colour = "grey80", size = 1) +
  geom_segment(x = "ARTPH", xend = "PTOF", y = dat_efigure_3_m$est[8], yend = dat_efigure_3_m$est[8], colour = "grey80", size = 1) +
  geom_segment(x = "PTOF", xend = "PLATE", y = dat_efigure_3_m$est[9], yend = dat_efigure_3_m$est[9], colour = "grey80", size = 1) +
  geom_segment(x = "PLATE", xend = "VASO", y = dat_efigure_3_m$est[10], yend = dat_efigure_3_m$est[10], colour = "grey80", size = 1) +
  geom_segment(x = "VASO", xend = "AADO2", y = dat_efigure_3_m$est[11], yend = dat_efigure_3_m$est[11], colour = "grey80", size = 1) + 
  geom_segment(x = "GLASGOW", y = dat_efigure_3_m$est[1], xend = "GLASGOW", yend = dat_efigure_3_m$est[2], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "PACO2", y = dat_efigure_3_m$est[2], xend = "PACO2", yend = dat_efigure_3_m$est[3], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "BILI", y = dat_efigure_3_m$est[3], xend = "BILI", yend = dat_efigure_3_m$est[4], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "TIDAL", y = dat_efigure_3_m$est[4], xend = "TIDAL", yend = dat_efigure_3_m$est[5], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "SYSBP", y = dat_efigure_3_m$est[5], xend = "SYSBP", yend = dat_efigure_3_m$est[6], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "AGE", y = dat_efigure_3_m$est[6], xend = "AGE", yend = dat_efigure_3_m$est[7], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "ARTPH", y = dat_efigure_3_m$est[7], xend = "ARTPH", yend = dat_efigure_3_m$est[8], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "PTOF", y = dat_efigure_3_m$est[8], xend = "PTOF", yend = dat_efigure_3_m$est[9], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "PLATE", y = dat_efigure_3_m$est[9], xend = "PLATE", yend = dat_efigure_3_m$est[10], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "VASO", y = dat_efigure_3_m$est[10], xend = "VASO", yend = dat_efigure_3_m$est[11], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "AADO2", y = dat_efigure_3_m$est[11], xend = "AADO2", yend = dat_efigure_3_m$est[12], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  annotate(geom = "text", x = "GLASGOW", y = dat_efigure_3_m$est[2] - 0.2, label = "-7.1%", color = "blue3") + 
  annotate(geom = "text", x = "PACO2", y = dat_efigure_3_m$est[3] - 0.2, label = "-16.5%", color = "blue3") + 
  annotate(geom = "text", x = "BILI", y = dat_efigure_3_m$est[4] + 0.2, label = "+2.9%", color = "red3") + 
  annotate(geom = "text", x = "TIDAL", y = dat_efigure_3_m$est[5] + 0.2, label = "+20.8%", color = "red3") + 
  annotate(geom = "text", x = "SYSBP", y = dat_efigure_3_m$est[6] + 0.2, label = "+15.7%", color = "red3") + 
  annotate(geom = "text", x = "AGE", y = dat_efigure_3_m$est[7] - 0.2, label = "-4.5%", color = "blue3") + 
  annotate(geom = "text", x = "ARTPH", y = dat_efigure_3_m$est[8] - 0.2, label = "-38.9%", color = "blue3") + 
  annotate(geom = "text", x = "PTOF", y = dat_efigure_3_m$est[9] + 0.2, label = "+6.7%", color = "red3") + 
  annotate(geom = "text", x = "PLATE", y = dat_efigure_3_m$est[10] + 0.2, label = "+13.0%", color = "red3") + 
  annotate(geom = "text", x = "VASO", y = dat_efigure_3_m$est[11] + 0.2, label = "+7.7%", color = "red3") +
  annotate(geom = "text", x = "AADO2", y = dat_efigure_3_m$est[12] - 0.2, label = "-19.7%", color = "blue3") + 
  geom_text(x = -Inf, y = Inf, label = "Patient from middle third", hjust = 0, vjust = 1)


dat_efigure_3_l = dat_figure_1a[dat_figure_1a$id == lower_third_example,]
dat_efigure_3_l_percent = (dat_efigure_3_l$est[-nrow(dat_efigure_3_l)] - dat_efigure_3_l$est[-1])/dat_efigure_3_l$est[-nrow(dat_efigure_3_l)]

lines_efigure_S1Sp_l = ggplot(data = dat_efigure_3_l[-1.], aes(x = predictor, y = est)) + 
  geom_point(color = "grey80", size = 1) +
  geom_point(x = "GLASGOW", y = dat_efigure_3_l$est[1], shape = 21, size = 5) +
  geom_point(x = "AADO2", y = dat_efigure_3_l$est[12], shape = 4, size = 5) +
  theme_bw() +
  scale_x_discrete(limits = c("GLASGOW", "PACO2", "BILI", "TIDAL", "SYSBP", "AGE", "ARTPH", "PTOF", "PLATE", "VASO", "AADO2"))  + 
  xlab("") + ylab("") + # ggtitle("ISCE Sp") + 
  theme(plot.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12)) + 
  geom_segment(x = "GLASGOW", xend = "PACO2", y = dat_efigure_3_l$est[2], yend = dat_efigure_3_l$est[2], colour = "grey80", size = 1) +
  geom_segment(x = "PACO2", xend = "BILI", y = dat_efigure_3_l$est[3], yend = dat_efigure_3_l$est[3], colour = "grey80", size = 1) +
  geom_segment(x = "BILI", xend = "TIDAL", y = dat_efigure_3_l$est[4], yend = dat_efigure_3_l$est[4], colour = "grey80", size = 1) +
  geom_segment(x = "TIDAL", xend = "SYSBP", y = dat_efigure_3_l$est[5], yend = dat_efigure_3_l$est[5], colour = "grey80", size = 1) +
  geom_segment(x = "SYSBP", xend = "AGE", y = dat_efigure_3_l$est[6], yend = dat_efigure_3_l$est[6], colour = "grey80", size = 1) +
  geom_segment(x = "AGE", xend = "ARTPH", y = dat_efigure_3_l$est[7], yend = dat_efigure_3_l$est[7], colour = "grey80", size = 1) +
  geom_segment(x = "ARTPH", xend = "PTOF", y = dat_efigure_3_l$est[8], yend = dat_efigure_3_l$est[8], colour = "grey80", size = 1) +
  geom_segment(x = "PTOF", xend = "PLATE", y = dat_efigure_3_l$est[9], yend = dat_efigure_3_l$est[9], colour = "grey80", size = 1) +
  geom_segment(x = "PLATE", xend = "VASO", y = dat_efigure_3_l$est[10], yend = dat_efigure_3_l$est[10], colour = "grey80", size = 1) +
  geom_segment(x = "VASO", xend = "AADO2", y = dat_efigure_3_l$est[11], yend = dat_efigure_3_l$est[11], colour = "grey80", size = 1) + 
  geom_segment(x = "GLASGOW", y = dat_efigure_3_l$est[1], xend = "GLASGOW", yend = dat_efigure_3_l$est[2], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "PACO2", y = dat_efigure_3_l$est[2], xend = "PACO2", yend = dat_efigure_3_l$est[3], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "BILI", y = dat_efigure_3_l$est[3], xend = "BILI", yend = dat_efigure_3_l$est[4], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "TIDAL", y = dat_efigure_3_l$est[4], xend = "TIDAL", yend = dat_efigure_3_l$est[5], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "SYSBP", y = dat_efigure_3_l$est[5], xend = "SYSBP", yend = dat_efigure_3_l$est[6], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "AGE", y = dat_efigure_3_l$est[6], xend = "AGE", yend = dat_efigure_3_l$est[7], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "ARTPH", y = dat_efigure_3_l$est[7], xend = "ARTPH", yend = dat_efigure_3_l$est[8], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "PTOF", y = dat_efigure_3_l$est[8], xend = "PTOF", yend = dat_efigure_3_l$est[9], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "PLATE", y = dat_efigure_3_l$est[9], xend = "PLATE", yend = dat_efigure_3_l$est[10], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  geom_segment(x = "VASO", y = dat_efigure_3_l$est[10], xend = "VASO", yend = dat_efigure_3_l$est[11], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "red3") + 
  geom_segment(x = "AADO2", y = dat_efigure_3_l$est[11], xend = "AADO2", yend = dat_efigure_3_l$est[12], lineend = "round", linejoin = "round",
               size = 1, arrow = arrow(length = unit(0.2, "inches")), color = "blue3") + 
  annotate(geom = "text", x = "GLASGOW", y = dat_efigure_3_l$est[2] - 0.5, label = "-119.3%", color = "blue3") + 
  annotate(geom = "text", x = "PACO2", y = dat_efigure_3_l$est[3] + 0.5, label = "+1.2%", color = "red3") + 
  annotate(geom = "text", x = "BILI", y = dat_efigure_3_l$est[4] - 0.5, label = "-11.0%", color = "blue3") + 
  annotate(geom = "text", x = "TIDAL", y = dat_efigure_3_l$est[5] + 0.5, label = "+4.5%", color = "red3") + 
  annotate(geom = "text", x = "SYSBP", y = dat_efigure_3_l$est[6] - 0.5, label = "-6.3%", color = "blue3") + 
  annotate(geom = "text", x = "AGE", y = dat_efigure_3_l$est[7] + 0.5, label = "+37.1%", color = "red3") + 
  annotate(geom = "text", x = "ARTPH", y = dat_efigure_3_l$est[8] - 0.5, label = "-50.4%", color = "blue3") + 
  annotate(geom = "text", x = "PTOF", y = dat_efigure_3_l$est[9] - 0.5, label = "-4.0%", color = "blue3") + 
  annotate(geom = "text", x = "PLATE", y = dat_efigure_3_l$est[10] - 0.5, label = "-6.4%", color = "blue3") + 
  annotate(geom = "text", x = "VASO", y = dat_efigure_3_l$est[11] + 0.5, label = "+20.6%", color = "red3") +
  annotate(geom = "text", x = "AADO2", y = dat_efigure_3_l$est[12] - 0.5, label = "-26.2%", color = "blue3") + 
  geom_text(x = -Inf, y = Inf, label = "Patient from lower third", hjust = 0, vjust = 1)

yleft = textGrob("Predicted individual treatment effect", rot = 90, gp = gpar(fontsize = 15))
bottom = textGrob("Baseline predictors", gp = gpar(fontsize = 15))

grid.arrange(arrangeGrob(lines_efigure_S1Sp_u, lines_efigure_S1Sp_m, lines_efigure_S1Sp_l, nrow = 3, ncol = 1), left = yleft, bottom = bottom)

SACE Analysis using BART

Xinyuan Chen

Introduction

Functions

Step 1: Estimating the SACE and CSACE using BART

Function: `sace_bart`

Parameters

Returns

Step 2: “Fit-the-Fit” Approach with the Classification and Regression Trees (CART) Model for CSACE Subgroup Identification

Function: `cart_fit`

Parameters

Returns

ARDS ARMA data analysis

Additional plots

SACE Analysis using BART

Xinyuan Chen

Introduction

Functions

Step 1: Estimating the SACE and CSACE using BART

Function: sace_bart

Parameters

Returns

Step 2: “Fit-the-Fit” Approach with the Classification and Regression Trees (CART) Model for CSACE Subgroup Identification

Function: cart_fit

Parameters

Returns

ARDS ARMA data analysis

Additional plots

Function: `sace_bart`

Function: `cart_fit`