IPWE_Qopt_IndCen {QTOCen}

Function to estimate the quantile-optimal treatment regime: the independent censoring Case

Description

This function implements the estimation method proposed in Chapter 2 of (Zhou 2018). It estimates the quantile-optimal treatment regime for a given quantile level of interest from a single-stage clinical randomized experiment or a single-stage observational study under the independent censoring assumption. In other words, we estimate the parameters indexing the quantile-optimal treatment regime.

Our assumption of independent censoring means the distribution of the censoring time is the same conditional on baseline covariates, treatment group and the two potential survival times.

Usage

IPWE_Qopt_IndCen(data, regimeClass, tau, moPropen = "BinaryRandom",
  Domains = NULL, cluster = FALSE, p_level = 1, s.tol = 1e-04,
  it.num = 8, pop.size = 6000, sign_beta1 = NULL,
  Penalty.level = 0)

Arguments

data

a data.frame, containing variables in the moPropen and RegimeClass and also the response variables, namely censor_y as the censored response, and delta as the censoring indicator.

regimeClass

a formula specifying the class of treatment regimes to search, e.g. if regimeClass = a~x1+x2, and then this function will search the class of treatment regimes of the form

d(x) = I \left(\beta_0 +\beta_1 x_1 + \beta_2 x_2 > 0\right).

Polynomial arguments are also supported.

tau

a value between 0 and 1. This is the quantile of interest.

moPropen

The propensity score model for the probability of receiving treatment level 1. When moPropen equals the string "BinaryRandom", the proportion of observations receiving treatment level 1 in the sample will be plugged in as an estimate of the propensity. Otherwise, this argument should be a formula/string, based on which this function will fit a logistic regression on the treatment level. e.g. a1~x1.

Domains

default is NULL. Otherwise, the object should be a nvars *2 matrix used as the space of parameters, which will be supplied to rgenoud::genoud.

cluster

default is FALSE, meaning do not use parallel computing for the genetic algorithm(GA).

p_level

choose between 0,1,2,3 to indicate different levels of output from the genetic function. Specifically, 0 (minimal printing), 1 (normal), 2 (detailed), and 3 (debug).

s.tol

tolerance level for the GA algorithm. This is input for parameter solution.tolerance in function rgenoud::genoud.

it.num

the maximum GA iteration number

pop.size

an integer with the default set to be 3000. This is roughly the number individuals for the first generation in the genetic algorithm (rgenoud::genoud).

sign_beta1

logical. Default is NULL. FALSE if the coefficient for the first continuous variable is fixed to be negative one; TRUE if positive one.

Penalty.level

0: stop if the marginal quantiles cannot be further optimized; 1: continue the search among treatment regimes with with same value for the TR with the smallest intended proportion of treatment.

Details

The input argument data is the dataframe that contains:

1. a observed treatment assignment

2. censor_y the censored response variable

3. delta the censoring indicator

The naming of these three columns should be strict.

Note that this function currently only works for scenarios in which treatment is binary.

References

Zhou Y (2018). Quantile-Optimal Treatment Regimes with Censored Data. Ph.D. thesis, University of Minnesota.

Horowitz JL (1992). “A smoothed maximum score estimator for the binary response model.” Econometrica: journal of the Econometric Society, 505–531.

Examples

     
GenerateData <- function(n)
{
  x1 <- runif(n, min=-0.5,max=0.5)
  x2 <- runif(n, min=-0.5,max=0.5)
  error <- rnorm(n, sd= 1)
  ph <- exp(-0.5+1*(x1+x2))/(1+exp(-0.5 + 1*(x1+x2)))
  a <- rbinom(n = n, size = 1, prob=ph)
  c <- 1 + 1*a + runif(n = n, min=0, max=2)
  cmplt_y <-  pmin(2+x1+x2 +  a*(1 - x1 - x2) +  (0.2 + a*(1+x1+x2)) * error, 4.4)
  censor_y <- pmin(cmplt_y, c)
  delta <- as.numeric(c > cmplt_y)
  return(data.frame(x1=x1,x2=x2,a=a, censor_y = censor_y, delta=delta))
}
n <- 400

data <- GenerateData(n)
fit1 <- IPWE_Qopt_IndCen(data = data, regimeClass = a~x1+x2, tau=0.25)

# We can used the returned model to visualize the Kaplan-meier
# estimate of survival function of the censoring time variable,
# justified by the independent censoring assumption.
library(survminer)
ggsurvplot(fit1$survfitCensorTime, data=fit1$data_aug, risk.table = TRUE)