# NOT RUN {
##########################################################################
# Note: the preprocessing steps prior to calling est_quant_TwoStg_ipwe() #
# are wrapped up in IPWE_Qopt_DTR_IndCen(). #
# w_di_vec is the inverse probability weight for two stage experiments #
# We recommend users to use function IPWE_Qopt_DTR_IndCen() directly. #
# Below is a simple customized calculation of the weight that only works #
# for this example #
##########################################################################
library(survival)
# Simulate data
n=200
s_Diff_Time = 1
D <- simJLSDdata(n, case="a")
# give regime classes
regimeClass.stg1 <- as.formula(a0~x0)
regimeClass.stg2 <- as.formula(a1~x1)
# extract columns that matches each stage's treatment regime formula
p.data1 <- model.matrix(regimeClass.stg1, D)
# p.data2 would only contain observations with non-null value.
p.data2 <- model.matrix(regimeClass.stg2, D)
txVec1 <- D[, "a0"]
# get none-na second stage treatment levels in data
txVec2 <- D[, "a1"]
txVec2_na_omit <- txVec2[which(!is.na(txVec2))]
# Eligibility flag
ELG <- (D$censor_y > s_Diff_Time)
# Build weights
D$deltaC <- 1 - D$delta
survfit_all <- survfit(Surv(censor_y, event = deltaC)~1, data=D)
survest <- stepfun(survfit_all$time, c(1, survfit_all$surv))
D$ghat <- survest(D$censor_y)
g_s_Diff_Time <- survest(s_Diff_Time)
D$w_di_vec <- rep(-999, n)
for(i in 1:n){
if (!ELG[i]) {
D$w_di_vec[i] <- 0.5 * D$ghat[i]} else {
D$w_di_vec[i] <- 0.5* D$ghat[i] * 0.5
}
}
qhat <- est_quant_TwoStg_ipwe(n=n, beta=c(2.5,2.8),
sign_beta1.stg1 = FALSE, sign_beta1.stg2=FALSE,
txVec1=txVec1, txVec2_na_omit=txVec2_na_omit, s_Diff_Time=1,
nvars.stg1=2, nvars.stg2=2,
p.data1=p.data1,
p.data2=p.data2,
censor_y=D$censor_y,
delta=D$delta,
ELG=ELG, w_di_vec=D$w_di_vec,
tau=0.3)
# }
Run the code above in your browser using DataLab