This function implements the estimator of two-stage quantile-optimal treatment regime by inverse probability of weighting proposed by Lan Wang, et al. As there are more than one stage, the second stage treatment regime could take into account the evolving status of an individual after the first stage and the treatment level received in the first stage. We assume the options at the two stages are both binary and take the form: $$d_1(x)=I\left(\beta_{10} +\beta_{11} x_{11} +...+ \beta_{1k} x_{1k} > 0\right), $$ $$d_2(x)=I\left(\beta_{20} +\beta_{21} x_{21} +...+ \beta_{2p} x_{2p} > 0\right)$$
TwoStg_Qopt(data, tau, regimeClass.stg1, regimeClass.stg2,
moPropen1 = "BinaryRandom", moPropen2 = "BinaryRandom", s.tol = 1e-04,
it.num = 8, max = TRUE, cl.setup = 1, p_level = 1, pop.size = 1000,
hard_limit = FALSE)a data frame, containing variables in the moPropen and RegimeClass and
a component y as the response.
a value between 0 and 1. This is the quantile of interest.
a formula or a string specifying the Class of treatment regimes
at stage 1, e.g. a1~x1+x2
a formula or a string specifying the Class of treatment regimes
at stage 2, e.g. a2~x1+a1+x2
The propensity score model for the probability of receiving
treatment level 1 at the first stage .
When moPropen1 equals the string "BinaryRandom", the proportion of observations
receiving treatment level 1 in the sample at the first stage 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.
The propensity score model for the probability of receiving
treatment level 1 at the second stage .
When moPropen2 equals the string "BinaryRandom", the proportion of observations
receiving treatment level 1 in the sample at the second stage 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. a2~x1+a1+x2.
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.
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.
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.
the number of nodes. >1 indicates choosing parallel computing option in
rgenoud::genoud. Default is 1.
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.)
an integer with the default set to be 3000. This is the population number for the first generation
in the genetic algorithm (rgenoud::genoud).
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.
This function returns an object with 7 objects. Both coefficients
and coef.orgn.scale were normalized to have unit euclidean norm.
coef.1, coef.2the set of parameters indexing the estimated quantile-optimal treatment regime for standardized covariates.
coef.orgn.scale.1, coef.orgn.scale.2the set of parameter indexing the estimated quantile-optimal treatment regime for the original input covariates.
tauthe quantile of interest
hatQthe estimated marginal quantile when the treatment
regime indexed by coef.orgn.scale.1 and coef.orgn.scale.2
is applied on the entire population.
See the 'details' for connection between coef.orgn.scale.k and
coef.k.
callthe user's call.
moPropen1, moPropen2the user specified propensity score models for the first and the second stage respectively
regimeClass.stg1, regimeClass.stg2the user specified class of treatment regimes for the first and the second stage respectively
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.
For every stage k, \(k=1,2\), this estimated parameters indexing the
two-stage quantile-optimal treatment regime are returned in two scales:
, the returned coef.k
is the set of parameters that we estimated after standardizing
every covariate available for decision-making
at stage k 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, coef.k, based on the standardized data.
For the identifiability issue, we force the Euclidean norm of coef.k
to be 1.
The difference between coef.k and coef.orgn.scale.k is that the latter
set of parameters correspond 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\{\beta_0 + \beta_1 x_1 + ... + \beta_k x_k > 0\},$$
where the \(\beta\) values are returned as coef.orgn.scale.k, and the
the vector \((1, x_1,...,x_k)\) corresponds to the specified class of treatment
regimes in the kth stage.
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
coef.k and coef.orgn.scale.k will render identical.
wang2017quantilequantoptr
# NOT RUN {
ilogit <- function(x) exp(x)/(1 + exp(x))
GenerateData.2stg <- function(n){
x1 <- runif(n)
p1 <- ilogit(-0.5+x1)
a1 <- rbinom(n, size=1, prob=p1)
x2 <- runif(n,x1,x1+1)
p2 <- ilogit(-1 + x2)
a2 <- rbinom(n, size=1, prob=p2)
mean <- 1+x1+a1*(1-3*(x1-0.2)^2) +x2 + a2*(1-x2-x1)
y <- mean + (1+a1*(x1-0.5)+0.5*a2*(x2-1))*rnorm(n,0,sd = 1)
return(data.frame(x1,a1,x2,a2,y))
}
# }
# NOT RUN {
n <- 400
testdata <- GenerateData.2stg(n)
fit <- TwoStg_Qopt(data=testdata, tau=0.2,
regimeClass.stg1=a1~x1, regimeClass.stg2=a2~x1+a1+x2,
moPropen1=a1~x1, moPropen2=a2 ~ x2,
cl.setup=2)
fit
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab