# The Doubly Robust Estimator of the Quantile-Optimal Treatment Regime `DR_Qopt` implements the doubly robust estimation method to estimate the quantile-optimal treatment regime. The double robustness property means that it is consistent when either the propensity score model is correctly specified, or the conditional quantile function is correctly specified. Both linear and nonlinear conditional quantile models are considered. See 'Examples' for an illustrative example. ```r DR_Qopt(data, regimeClass, tau, moPropen = "BinaryRandom", nlCondQuant_0 = FALSE, nlCondQuant_1 = FALSE, moCondQuant_0, moCondQuant_1, max = TRUE, length.out = 200, s.tol, it.num = 8, cl.setup = 1, p_level = 1, pop.size = 3000, hard_limit = FALSE, start_0 = NULL, start_1 = NULL) ``` ## Arguments - `data`: a data frame, must contain all the variables that appear in `moPropen`, `RegimeClass`, `moCondQuant_0`, `moCondQuant_1`, and a column named `y` as the observed 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).d(x)=I(\beta_0 +\beta_1 * x1 + \beta_2 * x2 > 0). $$ 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`. - `nlCondQuant_0`: Logical. When `nlCondQuant_0=TRUE`, this means the prespecified model for the conditional quantile function given a=0 is nonlinear, so the provided `moCondQuant_0` should be nonlinear. - `nlCondQuant_1`: Logical. When `nlCondQuant_1=TRUE`, this means the prespecified model for the conditional quantile function given a=1 is nonlinear, so the provided `moCondQuant_1` should be nonlinear. - `moCondQuant_0`: Either a formula or a string representing the parametric form of the conditional quantile function given that treatment=0. - `moCondQuant_1`: Either a formula or a string representing the parametric form of the conditional quantile function given that treatment=1. - `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`. - `length.out`: an integer greater than 1. If one of the conditional quantile model is set to be nonlinear, this argument will be triggered and we will fit `length.out` models across quantiles equally spaced between 0.001 and 0.999. Default is 200. - `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. - `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`). - `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`. - `start_0`: a named list or named numeric vector of starting estimates for the conditional quantile function when `treatment = 0`. This is required when `nlCondQuant_0=TRUE`. - `start_1`: a named list or named numeric vector of starting estimates for the conditional quantile function when `treatment = 1`. This is required when `nlCondQuant_1=TRUE`. ## Returns This function returns an object with 9 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 - **`moCondQuant_0`**: the user specified conditional quantile model for treatment 0 - **`moCondQuant_1`**: the user specified conditional quantile model for treatment 1 ## Details * Standardization on covariates AND explanation on the differences between the two returned regime parameters. 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 2016 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\}.d(x)= I{\beta_0 + \beta_1*x_1 + ... + \beta_k*x_k \\> 0}. $$ The estimated $\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 each 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. * Property of the doubly robust(DR) estimator. The DR estimator `DR_Qopt` is consistent if either the propensity score model or the conditional quantile regression model is correctly specified. (Wang et. al. 2016) ## Examples ```r ilogit <- function(x) exp(x)/(1 + exp(x)) GenerateData.DR <- function(n) { x1 <- runif(n,min=-1.5,max=1.5) x2 <- runif(n,min=-1.5,max=1.5) tp <- ilogit( 1 - 1*x1^2 - 1* x2^2) a <-rbinom(n,1,tp) y <- a * exp(0.11 - x1- x2) + x1^2 + x2^2 + a*rgamma(n, shape=2*x1+3, scale = 1) + (1-a)*rnorm(n, mean = 2*x1 + 3, sd = 0.5) return(data.frame(x1=x1,x2=x2,a=a,y=y)) } regimeClass <- as.formula(a ~ x1+x2) moCondQuant_0 <- as.formula(y ~ x1+x2+I(x1^2)+I(x2^2)) moCondQuant_1 <- as.formula(y ~ exp( 0.11 - x1 - x2)+ x1^2 + p0 + p1*x1 + p2*x1^2 + p3*x1^3 +p4*x1^4 ) start_1 = list(p0=0, p1=1.5, p2=1, p3 =0,p4=0) n <- 400 testdata <- GenerateData.DR(n) ## Examples below correctly specified both the propensity model and ## the conditional quantile model. system.time( fit1 <- DR_Qopt(data=testdata, regimeClass = regimeClass, tau = 0.25, moPropen = a~I(x1^2)+I(x2^2), moCondQuant_0 = moCondQuant_0, moCondQuant_1 = moCondQuant_1, nlCondQuant_1 = TRUE, start_1=start_1, pop.size = 1000)) fit1 ## Go parallel for the same fit. It would save a lot of time. ### Could even change the cl.setup to larger values ### if more cores are available. system.time(fit2 <- DR_Qopt(data=testdata, regimeClass = regimeClass, tau = 0.25, moPropen = a~I(x1^2)+I(x2^2), moCondQuant_0 = moCondQuant_0, moCondQuant_1 = moCondQuant_1, nlCondQuant_1 = TRUE, start_1=start_1, pop.size = 1000, cl.setup=2)) fit2 ``` ## References Rdpack::insert_ref(key="wang2017quantile",package="quantoptr") ## See Also `dr_quant_est`, `augX` ## Author(s) Yu Zhou, zhou0269@umn.edu

DR_Qopt function