TwoStg_Mopt function

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]$

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 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.

Examples

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