DR_Qopt {quantoptr} | R Documentation |
The Doubly Robust Estimator of the Quantile-Optimal Treatment Regime
Description
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.
Usage
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 |
regimeClass |
a formula specifying the class of treatment regimes to search,
e.g. if
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 |
nlCondQuant_0 |
Logical. When |
nlCondQuant_1 |
Logical. When |
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 |
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
|
s.tol |
This is the tolerance level used by |
it.num |
integer > 1. This argument will be used in |
cl.setup |
the number of nodes. >1 indicates choosing parallel computing option in
|
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 ( |
hard_limit |
logical. When it is true the maximum number of generations
in |
start_0 |
a named list or named numeric vector of starting estimates for
the conditional quantile function when |
start_1 |
a named list or named numeric vector of starting estimates for
the conditional quantile function when |
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:
The returned
coefficients
is the set of parameters after covariatesX
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 ofcoefficients
to be 1.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 ascoef.orgn.scale
. The same ascoefficients
, we force the Euclidean norm ofcoef.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
andcoef.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)
Value
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 betweencoef.orgn.scale
andcoefficient
.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
Author(s)
Yu Zhou, zhou0269@umn.edu
References
Wang L, Zhou Y, Song R and Sherwood B (2017). “Quantile-Optimal Treatment Regimes.” Journal of the American Statistical Association.
See Also
Examples
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