# Estimate the Two-stage Mean-Optimal Treatment Regime 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: [REMOVE_ME]$$ d_1(x_{stage1})=I\left(\beta_{10} +\beta_{11} x_{11} +...+ \beta_{1k} x_{1k} > 0\right),d1(x_stage1)=I(\beta_10 +\beta_11 * x11 +... + \beta_1k * x1k > 0), [REMOVE_ME_2]$$ [REMOVE_ME]$$ d_2(x_{stage2})=I\left(\beta_{20} +\beta_{21} x_{21} +...+ \beta_{2j} x_{2j} > 0\right)d2(x_stage2)=I(\beta_20 + \beta_21 * x21 +... + \beta_2j * x2j > 0) [REMOVE_ME_2]$$ ```r 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) ``` ## Arguments - `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`. ## Returns 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 ## Description 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: $$ d_1(x_{stage1})=I\left(\beta_{10} +\beta_{11} x_{11} +...+ \beta_{1k} x_{1k} > 0\right),d1(x_stage1)=I(\beta_10 +\beta_11 * x11 +... + \beta_1k * x1k > 0), $$ $$ d_2(x_{stage2})=I\left(\beta_{20} +\beta_{21} x_{21} +...+ \beta_{2j} x_{2j} > 0\right)d2(x_stage2)=I(\beta_20 + \beta_21 * x21 +... + \beta_2j * x2j > 0) $$ ## 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. For every stage `k`, $k=1,2$, this estimated parameters indexing the two-stage mean-optimal treatment regime are returned **in two scales:** 1. , 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. 2. 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\},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 `k`th 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. ## Examples ```r 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 ``` ## References Rdpack::insert_ref(key="zhang2013robust",package="quantoptr") ## Author(s) Yu Zhou, zhou0269@umn.edu

TwoStg_Mopt function