TwoStg_Mopt {quantoptr}

Estimate the Two-stage Mean-Optimal Treatment Regime

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

d_2(x_{stage2})=I\left(\beta_{20} +\beta_{21} x_{21} +...+ \beta_{2j} x_{2j} > 0\right)

Usage

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.

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\},

   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.

Value

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

Author(s)

Yu Zhou, [e-mail]

References

Zhang B, Tsiatis AA, Laber EB and Davidian M (2013). “Robust estimation of optimal dynamic treatment regimes for sequential treatment decisions.” Biometrika, 100(3).

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