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

# Load and process data set
data(bmiData)

# define the negative 12 month change in BMI from baseline
y12 <- -100*(bmiData[,6L] - bmiData[,4L])/bmiData[,4L]

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

# outcome model second stage
moMain2 <- buildModelObj(model = ~parentBMI+month4BMI,
                         solver.method = 'lm')

moCont2 <- buildModelObj(model = ~race + parentBMI+month4BMI,
                         solver.method = 'lm')

# outcome model first stage
moMain1 <- buildModelObj(model = ~parentBMI+baselineBMI,
                         solver.method = 'lm')

moCont1 <- buildModelObj(model = ~race + parentBMI+baselineBMI,
                         solver.method = 'lm')

moMain <- list(moMain1, moMain2)
moCont <- list(moCont1, moCont2)

# regime function second stage
regime2 <- function(eta1, eta2, data) {
             tst <- {data$parentBMI > eta1} & {data$month4BMI > eta2}
             rec <- rep('MR', nrow(x = data))
             rec[!tst] <- 'CD'
             return( rec )
           }

# regime function first stage
regime1 <- function(eta1, eta2, data) {
             tst <- {data$parentBMI > eta1} & {data$baselineBMI > eta2}
             rec <- rep('MR', nrow(x = data))
             rec[!tst] <- 'CD'
             return( rec )
           }

regimes <- list(regime1, regime2)

#### Analysis using AIPW
## Not run: 
fit_AIPW <- optimalSeq(moPropen = moPropen, 
                       moMain = moMain, moCont = moCont,
                       regimes = regimes,
                       data = bmiData, response = y12,  txName = c('A1','A2'),
                       Domains = cbind(rep(0,4),rep(100,4)),
                       pop.size = 100, starting.values = rep(25,4))

##Available methods

  # Coefficients of the regression objects
  coef(object = fit_AIPW)

  # Description of method used to obtain object
  DTRstep(object = fit_AIPW)

  # Estimated value of the optimal treatment regime for training set
  estimator(x = fit_AIPW)

  # Value object returned by regression methods
  fitObject(object = fit_AIPW)

  # Retrieve the results of genetic algorithm
  genetic(object = fit_AIPW)

  # Estimated optimal treatment and decision functions for training data
  optTx(x = fit_AIPW)

  # Estimated optimal treatment and decision functions for new data
  optTx(x = fit_AIPW, newdata = bmiData)

  # Value object returned by outcome regression method
  outcome(object = fit_AIPW)

  # Plots if defined by regression methods
  dev.new()
  par(mfrow = c(2,4))

  plot(x = fit_AIPW)
  plot(x = fit_AIPW, suppress = TRUE)

  # Retrieve the value object returned by propensity regression method
  propen(object = fit_AIPW)

  # Show main results of method
  show(object = fit_AIPW)

  # Show summary results of method
  summary(object = fit_AIPW)

## End(Not run)
#### Single Decision Point Analysis using IPW

# 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'))

# regime function second stage
regime <- function(eta1, eta2, data) {
            tst <- {data$parentBMI > eta1} & {data$month4BMI > eta2}
            rec <- rep('MR', nrow(x = data))
            rec[!tst] <- 'CD'
            return( rec )
          }
## Not run: 
fit_IPW <- optimalSeq(moPropen = moPropen, 
                      regimes = regime,
                      data = bmiData, response = y12,  txName = 'A2',
                      Domains = cbind(rep(0,2),rep(100,2)),
                      pop.size = 100, starting.values = rep(25,2))

##Available methods

  # Coefficients of the regression objects
  coef(object = fit_IPW)

  # Description of method used to obtain object
  DTRstep(object = fit_IPW)

  # Estimated value of the optimal treatment regime for training set
  estimator(x = fit_IPW)

  # Value object returned by regression method
  fitObject(object = fit_IPW)

  # Retrieve the results of genetic algorithm
  genetic(object = fit_IPW)

  # Estimated optimal treatment and decision functions for training data
  optTx(x = fit_IPW)

  # Estimated optimal treatment and decision functions for new data
  optTx(x = fit_IPW, newdata = bmiData)

  # Value object returned by outcome regression method
  outcome(object = fit_IPW)

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

  plot(x = fit_IPW)
  plot(x = fit_IPW, suppress = TRUE)

  # Retrieve the value object returned by propensity regression method
  propen(object = fit_IPW)

  # Show main results of method
  show(object = fit_IPW)

  # Show summary results of method
  summary(object = fit_IPW)

## End(Not run)