# example single run of a 2-stage g-estimation analysis
set.seed(1)
# expit function
expit <- function(x) { 1.0 / (1.0 + exp(-x)) }
# sample size
n <- 10000
# variables (X = patient information, A = treatment)
X1 <- rnorm(n)
A1 <- rbinom(n, 1, expit(X1))
X2 <- rnorm(n)
A2 <- rbinom(n, 1, expit(X2))
# blip functions
gamma1 <- A1 * (1 + X1)
gamma2 <- A2 * (1 + X2)
# observed outcome: treatment-free outcome plus blip functions
Y <- exp(X1) + exp(X2) + gamma1 + gamma2 + rnorm(n)
# models to be passed to DTRreg
# blip model
blip.mod <- list(~ X1, ~ X2)
# treatment model (correctly specified)
treat.mod <- list(A1 ~ X1, A2 ~ 1)
# treatment-free model (incorrectly specified)
tf.mod <- list(~ X1, ~ X2)
# perform G-estimation
mod1 <- DTRreg(twoStageCont$Y, blip.mod, treat.mod, tf.mod,
data = twoStageCont, method = "gest")
# predicted Y for optimal treatment
dat <- data.frame(X1, X2, A1, A2)
predict(mod1, newdata = dat)
Run the code above in your browser using DataLab