Estimate the marginal quantile response of a specific dynamic TR
Assume we have binary treatment options for two sequential stages with a fixed time duration between them. This means for each subject in the target population if the censored survival time or the time-to-event is beyond the timepoint of the second treatment.
This function evaluates a given dynamic treatment regime and returns the estimated marginal quantile response.
We assume the space of two-stage treatment regimes is a cartesian product of two single-stage linear treatment regime space.
The user facing function that applies this function is IPWE_Qopt_DTR_IndCen.
est_quant_TwoStg_ipwe(n, beta, sign_beta1.stg1, sign_beta1.stg2, txVec1, txVec2_na_omit, s_Diff_Time, nvars.stg1, nvars.stg2, p.data1, p.data2, censor_y, delta, ELG, w_di_vec, tau, check_complete = TRUE, Penalty.level = 0)
n: the sample size
beta: the vector of coefficients indexing a two-stage treatment regime
sign_beta1.stg1: Is sign of the coefficient for the first non-intercept variable for the first stage known? Default is NULL, meaning user does not have contraint on the sign; FALSE if the coefficient for the first continuous variable is fixed to be -1; TRUE if 1. We can make the search space discrete because we employ scale normalizaion.
sign_beta1.stg2: Default is NULL. Similar to sign_beta1.stg1.
txVec1: the vector of treatment received at the first stage
txVec2_na_omit: the vector of second stage treatment for patients who indeed second stage treatment
s_Diff_Time: the length of time between the first stage treatment and the second stage treatment
nvars.stg1: number of coeffients for the decision rule of the first stage
nvars.stg2: number of coeffients for the decision rule of the second stage
p.data1: the design matrix to be used for decision in stage one
p.data2: the design matrix to be used for decision in stage two
censor_y: Numeric vector. The censored survival times from all observed data, i.e. censor_y = min(Y, C)
delta: Numeric vector. The censoring indicators from all observed data. We use 1 for uncensored, 0 for censored.
ELG: the boolean vector of whether patients get the second stage treatment
w_di_vec: the inverse probability weight for two stage experiments
tau: a value between 0 and 1. This is the quantile of interest.
check_complete: logical. Since this value estimation method is purely nonparametric, we need at least one unit in collected data such that the observed treatment assignment is the same what the regime parameter suggests. If check_complete
is TRUE. It will check if any observation satisfies this criterion. When none observation satisfies, a message is printed to console to raise users awareness that the input regime parameter beta does not agree with any observed treatment level assignment. Then a sufficiently small number is returned from this function, to keep the genetic algorithm running smoothly.
Penalty.level: the level that determines which objective function to use. Penalty.level = 0 indicates no regularization; Penalty.level = 1 indicates the value function estimation minus the means absolute average coefficient is the output, which is useful trick to achieve uniqueness of estimated optimal TR when resolution of input response is low.
########################################################################## # 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)
Useful links