TwoStg_Mopt: Estimate the Two-stage Mean-Optimal Treatment Regime
In quantoptr: Algorithms for Quantile- And Mean-Optimal Treatment Regimes

Description Usage Arguments Details Value Author(s) References Examples

This function implements the estimator of two-stage mean-optimal treatment regime by inverse probability of weighting proposed by Baqun Zhang. 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:

d1(x_stage1)=I(β_10 +β_11 * x11 +... + β_1k * x1k > 0),

d2(x_stage2)=I(β_20 + β_21 * x21 +... + β_2j * x2j > 0)

TwoStg_Mopt(data, regimeClass.stg1, regimeClass.stg2,
  moPropen1 = "BinaryRandom", moPropen2 = "BinaryRandom", max = TRUE,
  s.tol, cl.setup = 1, p_level = 1, it.num = 10, pop.size = 3000,
  hard_limit = FALSE)

`data`	a data frame, containing variables in the `moPropen` and `RegimeClass` and a component `y` as the response.
`regimeClass.stg1`	a formula or a string specifying the Class of treatment regimes at stage 1, e.g. `a1~x1+x2`
`regimeClass.stg2`	a formula or a string specifying the Class of treatment regimes at stage 2, e.g. `a2~x1+a1+x2`
`moPropen1`	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`.
`moPropen2`	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`.
`max`	logical. If `max=TRUE`, it indicates we wish to maximize the marginal mean; if `max=FALSE`, we wish to minimize the marginal mean. 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`.
`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.)
`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.
`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`.

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 mean-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 2016 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 β denote the estimated parameter in the original scale, then the estimated treatment regime is:

d(x)= I{β_0 + β_1*x_1 + ... + β_k*x_k > 0},

where the β 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.

This function returns an object with 6 objects. Both coef.1, coef.2 and coef.orgn.scale.1, coef.orgn.scale.2 were normalized to have unit euclidean norm.

coef.1, coef.2: the set of parameters indexing the estimated mean-optimal treatment regime for standardized covariates.
coef.orgn.scale.1, coef.orgn.scale.2: the set of parameter indexing the estimated mean-optimal treatment regime for the original input covariates.
hatM: the estimated marginal mean 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.
call: the user's call.
moPropen1, moPropen2: the user specified propensity score models for the first and the second stage respectively
regimeClass.stg1, regimeClass.stg2: the user specified class of treatment regimes for the first and the second stage respectively

Yu Zhou, zhou0269@umn.edu

\insertRef

zhang2013robustquantoptr

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

n <- 400
testdata <- GenerateData.2stg(n)

fit <- TwoStg_Mopt(data=testdata, 
                   regimeClass.stg1="a1~x1", regimeClass.stg2="a2~x1+a1+x2",
                   moPropen1="a1~x1", moPropen2="a2~x2",
                   cl.setup=2)
fit

fit2 <- TwoStg_Mopt(data=testdata, 
                   regimeClass.stg1="a1~x1", regimeClass.stg2="a2~a1+x1*x2",
                   moPropen1="a1~x1", moPropen2="a2~x2",
                   cl.setup=2)
fit2