IPWE_Qopt {quantoptr}

Estimate the Quantile-optimal Treatment Regime

Description

Estimate the Quantile-optimal Treatment Regime by inverse probability of weighting

Usage

IPWE_Qopt(data, regimeClass, tau, moPropen = "BinaryRandom", max = TRUE,
  s.tol, it.num = 8, hard_limit = FALSE, cl.setup = 1, p_level = 1,
  pop.size = 3000)

Arguments

data

a data frame, containing variables in the moPropen and RegimeClass and a component y as the response.

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. See also 'Details'.

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 employed as a good estimate of the probability for each observation. 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.

max

logical. If max=TRUE, it indicates we wish to maximize the marginal quantile; if max=FALSE, we wish to minimize the marginal quantile. The default is TRUE.

s.tol

This is the tolerance level used by genoud. Default is 10^{-5} times the difference between the largest and the smallest value in the observed responses. This is particularly important when it comes to evaluating it.num.

it.num

integer > 1. This argument will be used in rgeound::geound function. If there is no improvement in the objective function in this number of generations, rgenoud::genoud will think that it has found the optimum.

hard_limit

logical. When it is true the maximum number of generations in rgeound::geound cannot exceed 100. Otherwise, in this function, only it.num softly controls when genoud stops. Default is FALSE.

cl.setup

the number of nodes. >1 indicates choosing parallel computing option in rgenoud::genoud. Default is 1.

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.)

pop.size

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

Details

Note that all estimation functions in this package use the same type of standardization on covariates. Doing so would allow us to provide a bounded domain of parameters for searching in the genetic algorithm.

This estimated parameters indexing the quantile-optimal treatment regime are returned in two scales:

1. The returned coefficients is the set of parameters after covariates X are standardized to be in the interval [0, 1]. To be exact, every covariate is subtracted by the smallest observed value and divided by the difference between the largest and the smallest value. Next, we carried out the algorithm in Wang et al. 2017 to get the estimated regime parameters, coefficients, based on the standardized data. For the identifiability issue, we force the Euclidean norm of coefficients to be 1.

2. In contrast, coef.orgn.scale corresponds to the original covariates, so the associated decision rule can be applied directly to novel observations. In other words, let \beta denote the estimated parameter in the original scale, then the estimated treatment regime is:

   d(x)= I\{\hat{\beta}_0 + \hat{\beta}_1 x_1 + ... + \hat{\beta}_k x_k > 0\}.

   The estimated \bm{\hat{\beta}} is returned as coef.orgn.scale. The same as coefficients, we force the Euclidean norm of coef.orgn.scale to be 1.

If, for every input covariate, the smallest observed value is exactly 0 and the range (i.e. the largest number minus the smallest number) is exactly 1, then the estimated coefficients and coef.orgn.scale will render identical.

Value

This function returns an object with 7 objects. Both coefficients and coef.orgn.scale were normalized to have unit euclidean norm.

coefficients

the parameters indexing the estimated quantile-optimal treatment regime for standardized covariates.

coef.orgn.scale

the parameter indexing the estimated quantile-optimal treatment regime for the original input covariates.

tau

the quantile of interest

hatQ

the estimated marginal tau-th quantile when the treatment regime indexed by coef.orgn.scale is applied on everyone. See the 'details' for connection between coef.orgn.scale and coefficient.

call

the user's call.

moPropen

the user specified propensity score model

regimeClass

the user specified class of treatment regimes

Author(s)

Yu Zhou, [e-mail] with substantial contribution from Ben Sherwood.

References

Wang L, Zhou Y, Song R and Sherwood B (2017). “Quantile-Optimal Treatment Regimes.” Journal of the American Statistical Association.

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= 0.5)
  tp <- exp(-1+1*(x1+x2))/(1+exp(-1+1*(x1+x2)))
  a <- rbinom(n = n, size = 1, prob=tp)
  y <-  1+x1+x2 +  a*(3 - 2.5*x1 - 2.5*x2) +  (0.5 + a*(1+x1+x2)) * error
  return(data.frame(x1=x1,x2=x2,a=a,y=y))
}
n <- 300
testData <- GenerateData(n)

# 1. Estimate the 0.25th-quantile optimal treatment regime. ###

fit1 <- IPWE_Qopt(data = testData, regimeClass = "a~x1+x2",
           tau = 0.25, moPropen="a~x1+x2")
fit1


# 2. Go parallel. This saves time in calculation. ###

fit2 <- IPWE_Qopt(data = testData, regimeClass = "a~x1+x2",
           tau = 0.25, moPropen="a~x1+x2", cl.setup=2)
fit2




# 3. Set a quardratic term in the class #######################

fit3 <- IPWE_Qopt(data = testData, regimeClass = "a~x1+x2+I(x1^2)",
                  tau = 0.25, moPropen="a~x1+x2", pop.size=1000)
fit3


# 4. Set screen prints level. #######################
# Set the p_level to be 0, 
# then all screen prints from the genetic algorithm will be suppressed.

fit4 <- IPWE_Qopt(data = testData, regimeClass = "a~x1+x2",
           tau = 0.25, moPropen="a~x1+x2", cl.setup=2, p_level=0)
fit4