R: DTR estimation and inference for time-to-event data using...

DWSurv {DTRreg}

R Documentation

DTR estimation and inference for time-to-event data using DWSurv

Description

Dynamic treatment regimen estimation and inference via dynamic weighted survival modeling (DWSurv). Inference for the blip estimators with single- and multi-stage data.

The function DWSurv() allows estimating an optimal dynamic treatment regime from multi-stage trials or observational data when the outcome of interest is survival time subject to right-censoringg. The dynamic weighted survival modeling (DWSurv) algorithm is implemented. The method focuses on estimating the parameters of the blip: a model of the difference in expected outcome under the observed treatment and some reference treatment (usually a control) at a given stage, assuming identical histories and optimal treatment thereafter.

The method requires the specification of four models for each stage of the analysis: a treatment model (conditional mean of the treatment variable), a censoring model, a treatment-free model (conditional mean of outcome assuming only reference treatments are used), and a blip model. Only the blip model must be correctly specified (or over-specified), with consistent parameter estimates obtainable if at least one of the treatment-free or the treatment and censoring models are correctly specified. Note that all of these must be specified as lists of formula objects, even if only one stage of treatment is considered.

Note that as is conventional, it is assumed a larger survival time is preferred (which can be easily achieved via transformation of your data if necessary).

Several treatment weight function options have been implemented within the package:

"none": No treatment weights applied. If method = "dWOLS", this selection results in the implementation of Q-learning, modified slightly to use the dWOLS style pseudo-outcome (computed using the observed outcome modified by the estimated treatment effect) rather than the traditional Q-learning outcome (predicted based on model only, rather than observed outcome with treatment effect).
"ipw": weights based on the inverse probability of treatment. For binary treatments, a logistic regression is used. For multinomial, a multinomial log-linear model is fit using multinom. For continuous treatments, a GLM with the specified family and link function provided in the treat.fam argument is used.
"cipw": inverse probability of treatment weights as described for "ipw" and capped at the 99th percentile of the observed weights.
"qpom": weights based on the stabilized inverse probability of treatment applied to the categorized (into n.bins bins) continuous doses or multinomial treatments; probabilities are calculated using a proportional odds model. This weight is appropriate only for continuous and multinomial treatments.
"wo": overlap weights for the categorized continuous doses or multinomial treatments (Li and Li, 2019). This weight is appropriate only for continuous treatments.
"abs": Absolute difference |A - E[A|...]|. appropriate only for binary treatments.
"manual": User provides treatment weights through input treat.wgt.man.
"manual.with.censor": User provides combined treatment * censoring weights through input treat.wgt.man. Note that 'cens.mod' should be specified with the event indicator on the right-hand side of the formula (e.g., ~ status).

Usage

DWSurv(
  time,
  blip.mod,
  treat.mod,
  tf.mod,
  cens.mod,
  data = NULL,
  method = c("dwols", "qlearn"),
  interactive = FALSE,
  treat.type = c("bin", "multi", "cont"),
  treat.fam = gaussian(link = "identity"),
  weight = c("abs", "ipw", "cipw", "qpom", "wo", "none", "manual", "manual.with.censor"),
  n.bins = 3L,
  treat.range = NULL,
  treat.wgt.man = NULL,
  var.estim = c("none", "bootstrap", "sandwich"),
  bootstrap.controls = list(B = 100L, M = 0L, type = "standard", truncate = 0, verbose =
    FALSE, interrupt = FALSE),
  dtr = TRUE,
  full.cov = FALSE
)

Arguments

`time`	A list of formula specifying the survival time variable for each stage in order. The time variable should be specified on the right hand side of the formula. No dependent variable should be specified. The list should be as long as the number of stages.
`blip.mod`	A list of formula objects specifying covariates of the blip function for each stage in order. No dependent variable should be specified. Note that this input should include the treatment variable ONLY if the blip model is quadratic in treatment. See Details for further clarification.
`treat.mod`	A list of formula objects specifying the treatment model for each stage in order. Treatment variable should be included as the dependent variable. If treatment is binary `glm(family = binomial)` will be used to obtain parameter estimates; if treatment is multi-nomial, `multinom()` will be used to obtain parameter estimates; and if treatment is continuous, `lm()` will be used.
`tf.mod`	A list of formula objects specifying covariates of the treatment-free model for each stage in order. No dependent variable should be specified.
`cens.mod`	A list of formula objects specifying the censoring model for each stage in order. The event indicator, which takes value 1 if an event was observed and 0 otherwise, should be included as the dependent variable and should be the same across stages. In the absence of censoring or if censoring weights are provided by the user through 'treat.wgt.man', (i.e., `weight = 'manual.with.censor'`) one still needs to specify an event indicator on the right-hand side of the formula and leave the left-hand side empty (see example below).
`data`	A data frame containing all necessary covariates and treatments contained in the models. Missing data should be coded as `NA`.
`method`	The DTR method to be used, choose "dwols" for dynamic WOLS, or "qlearn" for Q-learning.
`interactive`	If `TRUE` on-screen prompts will guide the user through the specification of blip, treatment, and treatment-free models.
`treat.type`	A character object. Must be one of {"bin", "multi", "cont"} indicating that the treatments at each stage are binary, multinomial, or continuous, respectively. Each stage must have the same treatment type.
`treat.fam`	A character or family object. The description of the dose distribution along with the link function to be used in the treatment model for computing weights; should be specified in a similar format as that used in `glm()`. If character object, must be one of {"gaussian", "Gamma"}, for which `gaussian(link = "identity")` or `Gamma(link = "log")` will be used, respectively. Input is ignored for `treat.type = "bin"` and `treat.type = "multi"`.
`weight`	The form of the treatment weight. See details.
`n.bins`	An integer object. The number of bins (levels) to be used for categorizing continuous doses. This input is required only when `treat.type = "cont"` and `weight = "wo"` or `weight = "qpom"`.
`treat.range`	For continuous treatments. Specify the maximum/minimum value that treatments can be take. If unspecified then the minimum/maximum value of observed treatments is used. If you wish to have unrestricted treatments set this option to `c(-Inf, +Inf)`. If each stage has its own range, provide as a list, the ith element providing the min and max for the ith stage treatment.
`treat.wgt.man`	NULL or a list of vectors of known treatment (or treatment * censoring) weights can be specified to be used instead of hard-coded treatment weight options. The `i^{th}` element of the list contains the multiplicative weights for the `i^{th}` stage. Each vector must be of length `n`, the number of participants. Used only for `method = "dwols"`. If providing the treatment * censoring weights, `cens.mod = NA` must be used.
`var.estim`	Covariance matrix estimation method, either "bootstrap" or "sandwich" for sandwich estimation.
`bootstrap.controls`	A named list specifying control parameters of the bootstrap if `var.estim = "bootstrap"`. Available controls are: B: The number of bootstrap samples. M: The subsample size for m out of n bootstrap. type: The type of bootstrap. Must be one of {"standard", "empirical", "normal"}. The last two are parametric bootstraps. truncate: A number between 0 and 0.5. The lowest and highest specified proportion of parameter estimates will be replaced by the relevant quantiles affording some robustness to extreme values when estimating covariance. verbose: If TRUE, estimated time to completion will be printed to the console every ~30 seconds. interrupt: If TRUE then user will be given the option to abort the bootstrap without error if estimated time to completion exceeds 10 minutes.
`dtr`	A logical object. If `TRUE`, use the DTR estimation approach, which estimates the stage pseudo-outcome by adding a regret function. If `FALSE`, use an 'effect estimation' approach, which treats the observed outcome as being equal to an outcome assuming no treatment is received at any stage, plus a blip component at each stage; each stage pseudo-outcome is generated by subtracting a blip function. Note that most of the DTR-specific output will either be suppressed or irrelevant.
`full.cov`	A logical. If `TRUE`, the full covariance matrix will be returned. If `FALSE`, only the terms pertaining to the blip parameters are returned.

Value

An object of class DWSurv, a list including elements

`K:`	The number of decision points.
`beta:`	A list. The ith element contains the parameter estimates of the ith stage treatment-free model.
`psi:`	A list. The ith element contains the parameter estimates of the ith stage blip model.
`covmat:`	A list. The ith element contains covariance matrix of the ith stage blip parameter estimates.
`nonreg:`	Non-regularity estimates.
`setup:`	A list detailing the input parameter settings used for the analysis

models: A list of the models used for the analysis.
method: The parameter estimation method.
var.estim: The variance estimation method.
cc.modeled: If TRUE, missing data was modeled. If FALSE, cases with missing data were removed from the analysis.
tx.weight: The treatment weighting used for the analysis.
tx.type: Treatment was binary, multinomial, or continuous.
n.bins: The number of bins (levels) used for categorizing continuous doses when tx.weight = "wo" or tx.weight = "qpom".
tx.wgt.man: Any user provided treatment weights.
tx.range: For continuous treatments, the range of allowed treatment values.
tx.family: The description of the dose distribution along with the link function used in the continuous treatment model.
boot.controls: A list of the bootstrap controls.
type: The type of effect. Dynamic treatment regime or treatment effect.

training_data:

A list containing the training data.

data: The covariates and treatment data.
outcome: The outcome of interest.
A: The treatment variables, possibly recoded to adhere to internal code requirements.

analysis:

A list containing the primary results of each stage analysis.

n: The number of participants included in the stage analysis.
last.stage: The last stage each participant was included in the analysis.
prob.cens: The complete case probabilities.
cens.mod.fitted: The regression objects returned for estimating the complete case probabilities.
cens.wgt: The complete case weights.
cts: The treatment type at each stage.
tx.mod.fitted: The regression objects returned for estimating the treatment probabilities.
A.hat: The estimated or provided treatment probabilities.
tx.wgt: The treatment weights.
outcome.fit: The regression objects returned for each stage outcome regression.
Y: The pseudo-outcomes.
regret: Estimates of the regret for each subject based on observed treatment and blip parameter estimates.
opt.treat: Optimal treatment decisions for each subject at each stage of treatment.
opt.Y: Predicted optimal outcome under recommended regimen.

call:

The original function call.

The functions coef(), predict() and confint() may be used with such model objects. The first two have specific help files for their implementation, while confint() is used in the same way as the standard confint() command, with the exception of the parm option, which is not available.

References

Simoneau, G., Moodie, E. E. M., Wallace, M.P., Platt, R. W. (2020) Optimal Dynamic Treatment Regimes with Survival Endpoints: Introducing DWSurv in the R package DTRreg. Journal of Statistical Computation and Simulation. 90, 2991-3008. (doi:10.1080/00949655.2020.1793341)

Simoneau, G., Moodie, E. E. M., Nijjar, J. S., Platt, R. W. (2019) Estimating Optimal Dynamic Treatment with Survival Outcomes. Journal of the American Statistical Association, 115, 1531-1539 (doi:10.1080/01621459.2019.1629939).

Wallace, M. P., Moodie, E. E. M., Stephens, D. A. (2017) Dynamic Treatment Regimen Estimation via Regression-Based Techniques: Introducing R Package DTRreg. Journal of Statistical Software 80(2), 1–20 (doi:10.18637/jss.v080.i02).

Simoneau, G., Moodie, E. E. M., Nijjar, J. S., and Platt, R. W. (2020) Finite Sample Variance Estimation for Optimal Dynamic Treatment Regimes of Survival Outcomes. Statistics in Medicine 39, 4466-4479.

Efron, B., and Tibshirani, R. (1986) Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy Source: Statistical Science 1 54-75.

Examples

#### example single run of a 2-stage DWSurv analysis
data(twoStageCens)
mod <- DWSurv(time = list(~ T1, ~ T2), 
              blip.mod = list(~ X11, ~ X21), 
              treat.mod = list(A1 ~ X11, A2 ~ 1), 
              tf.mod = list(~ X11 + X12, ~ X21 + X22 + X11), 
              cens.mod = list(delta ~ 1, delta ~ X11), 
              var.estim = "sandwich", 
              data = twoStageCens)
mod
  
#### example in the absence of censoring
data(twoStageSurv)
mod_nocensoring <- DWSurv(time = list(~ T1, ~ T2), 
                          blip.mod = list(~ X11, ~ X21), 
                          treat.mod = list(A1 ~ X11, A2 ~ 1), 
                          tf.mod = list(~ X11 + X12, ~ X21 + X22 + X11), 
                          cens.mod = list(~ delta, ~ delta), 
                          var.estim = "sandwich", 
                          data = twoStageSurv)
mod_nocensoring

[Package DTRreg version 2.2 Index]