DynTxRegime (version 3.2)

optimalSeq: Regression Based Value-Search Estimation of Optimal Dynamic Treatment Regimes

Description

A doubly robust Augmented Inverse Propensity Weighted Estimator (AIPWE) or Inverse Propensity Weighted Estimator (IPWE) for population mean outcome is optimized over a restricted class of regimes. Methods are available for both single-decision-point and multiple-decision-point regimes. This method requires the rgenoud package.

Usage

optimalSeq(..., moPropen, moMain, moCont, data, response, txName, regimes,
           fSet = NULL, refit = FALSE, iter = 0, verbose = TRUE)

Arguments

...

additional arguments required by rgenoud. At a minimum, this should include Domains, pop.size and starting.values. See ?rgenoud for more information.

moPropen

An object of class "modelObj," an object of class "list" containing objects of class "modelObj," or an object of class "list" containing objects of class "modelObjSubset." This object specifies the model(s) of the propensity score regression and the R methods to be used to obtain parameter estimates and predictions. The method(s) specified to obtain predictions must return the prediction on the scale of the probability, i.e., returned values must be in the interval (0,1). If performing a single-decision-point analysis, moPropen is an object of class "modelObj" or a list of objects of class "modelObjSubset." If performing a multiple-decision-point analysis, moPropen is a list of objects of class "modelObj" or a list of objects of class "modelObjSubset." See moPropen for further information.

moMain

An object of class "modelObj," an object of class "list" containing objects of class "modelObj," or an object of class "list" containing objects of class "modelObjSubset." This object specifies the model(s) of the main effects component of the outcome regression(s) and the R methods to be used to obtain parameter estimates and predictions. The method chosen to obtain predictions must return the prediction on the scale of the response variable. If performing a single-decision-point analysis, moMain is an object of class "modelObj" or a list of objects of class "modelObjSubset." If performing a multiple-decision-point analysis, moMain is a list of objects of class "modelObj" or a list of objects of class "modelObjSubset."

moCont

An object of class "modelObj," an object of class "list" containing objects of class "modelObj," or an object of class "list" containing objects of class "modelObjSubset." This object specifies the model(s) of the contrasts component of the outcome regression(s) and the R methods to be used to obtain parameter estimates and predictions. The method chosen to obtain predictions must return the prediction on the scale of the response variable. If performing a single-decision-point analysis, moCont is an object of class "modelObj" or a list of objects of class "modelObjSubset." If performing a multiple-decision-point analysis, moCont is a list of objects of class "modelObj" or a list of objects of class "modelObjSubset."

data

An object of class "data.frame." The covariates and treatment histories.

response

An object of class "vector." A vector of the outcome of interest.

txName

An object of class "character." For a single-decision-point analysis, the column header of the stage treatment variable as given in input data. For multiple-decision-point analyses, a vector, the ith element of which gives the column header of data containing the treatment variable for the ith stage.

regimes

An object of class "function" or an object of class "list" containing objects of class "function." Function(s) defining the class of decision rule to be considered. See Details section for further information

fSet

An object of class "function," an object of class "list" containing objects of class "function," or NULL. See fSet section for further information.

refit

No longer used.

iter

An object of class "numeric." If >0, the iterative method will be used to obtain parameter estimates in the outcome regression step. See iter for further information.

verbose

An object of class "logical." If FALSE, screen prints will be suppressed.

Value

Returns an object of class "OptimalSeq" that inherits directly from class "DynTxRegime."

Methods

coef

signature(object = "OptimalSeq"): Retrieve parameter estimates for all regression steps.

DTRstep

signature(object = "OptimalSeq"): Retrieve description of method used to create object.

estimator

signature(x = "OptimalSeq"): Retrieve the estimated value of the estimated optimal regime for the training data set.

fitObject

signature(object = "OptimalSeq"): Retrieve value object returned by regression methods.

genetic

signature(object = "OptimalSeq"): Retrieve value object returned by genoud().

optTx

signature(x = "OptimalSeq", newdata = "missing"): Retrieve the estimated optimal treatment regime for training data set.

optTx

signature(x = "OptimalSeq", newdata = "data.frame"): Estimate the optimal treatment regime for newdata.

outcome

signature(x = "OptimalSeq"): Retrieve value object returned by outcome regression methods.

plot

signature(x = "OptimalSeq"): Generate plots for regression analyses.

print

signature(object = "OptimalSeq"): Print main results of analysis.

propen

signature(x = "OptimalSeq"): Retrieve value object returned by propensity score regression methods.

show

signature(object = "OptimalSeq"): Show main results of analysis.

summary

signature(object = "OptimalSeq"): Retrieve summary information from regression analyses.

Details

The IPWE estimator can be chosen by specifying moMain and moCont as NULL.

The method uses a genetic algorithm (R package rgenoud) to maximize the AIPWE/IPWE estimator over the class of treatment regimes specified by the treatment rules.

A regimes function defines the class of the decision rule to be considered. The formal arguments of each function must include the parameters to be estimated and the data.frame. NOTE: THE LAST ARGUMENT OF THE FUNCTION MUST BE THE DATA FRAME. For example, for $$d_i = I(eta_1 > x1),$$

regimes <- function(a,data){
  as.numeric(a > data\$x1)
} 

For a single-decision-point analysis, regimes is a single function. For multiple-decision-point analyses, regimes is a list of functions where each element of the list corresponds to the decision point (1st element <- 1st decision point, etc.)]

References

Zhang, B., Tsiatis, A. A., Laber, E. B., and Davidian, M. (2012). A Robust Method for Estimating Optimal Treatment Regimes. Biometrics, 68, 1010--1018.

Zhang, B., Tsiatis, A. A., Laber, E. B., and Davidian, M. (2013) Robust Estimation of Optimal Dynamic Treatment Regimes for Sequential Treatment Decisions. Biometrika, 100, 681--694.

Mebane, W. and Sekhon, J. S. (2011). Genetic Optimization Using Derivatives : The rgenoud package for R. Journal of Statistical Software, 42, 1--26.

Examples

# NOT RUN {
# Load and process data set
  data(bmiData)

  # define response y to be the negative 12 month
  # change in BMI from baseline
  bmiData$y <- -100*(bmiData[,6] - bmiData[,4])/bmiData[,4]

# Define the propensity for treatment model and methods.
  moPropen <- buildModelObj(model =  ~ 1, 
                            solver.method = 'glm', 
                            solver.args = list('family'='binomial'),
                            predict.method = 'predict.glm',
                            predict.args = list(type='response'))


# Create modelObj object for main effect component
  moMain <- buildModelObj(model = ~ gender + parentBMI + month4BMI,
                          solver.method = 'lm')

# Create modelObj object for contrast component
  moCont <- buildModelObj(model = ~ parentBMI + month4BMI,
                          solver.method = 'lm')

# treatment regime rules at each decision point.
  regimes <- function(a,b,c,data){
               tst <- a + b*data$parentBMI + c*data$month4BMI > 0 
               res <- character(nrow(data))
               res[tst] <- "MR"
               res[!tst] <- "CD"
               return(res)
             }

# genoud requires some additional information 
  c1 <- c(-1,-1,-1)
  c2 <- c( 1, 1, 1)
  Domains <- cbind(c1,c2)
  starts <- c(0,0,0)

#!! A LARGER VALUE FOR POP.SIZE IS RECOMMENDED
#!! THIS VALUE WAS CHOSEN TO MINIMIZE RUN TIME OF EXAMPLES
  pop.size <- 50
# }
# NOT RUN {
  ft <- optimalSeq(moPropen = moPropen, moMain = moMain, moCont = moCont, 
                   data = bmiData,  response = bmiData$y,  txName = "A2", 
                   regimes = regimes, iter = 0L,
                   pop.size = pop.size, starting.values = starts, 
                   Domains = Domains, solution.tolerance = 0.0001)

## Available methods

  # Coefficients of the propensity score and outcome regression
  coef(ft)

  # Description of method used to obtain object
  DTRstep(ft)

  # Estimated value of estimated optimal treatment for training data
  estimator(ft)

  # Value object returned by outcome regression method
  fitObject(ft)

  # Value object returned by genetic algorithm
  genetic(ft)

  # Estimated optimal treatment for training data
  optTx(ft)

  # Estimated optimal treatment for new data
  optTx(ft, newdata = bmiData)

  # Value object returned by outcome regression method
  outcome(ft)

  # Plots if defined by outcome regression method
  dev.new()
  par(mfrow = c(2,4))
  plot(ft)

  dev.new()
  par(mfrow = c(2,4))
  plot(ft, suppress = TRUE)

  # Value object returned by propensity score regression method
  propen(ft)

  # Estimated decision function parameters
  regimeCoef(ft)

  # Show main results of method
  show(ft)

  # Show summary results of method
  summary(ft)

# }
# NOT RUN {
# }