DynTxRegime An R Package for Dynamic Treatment Regimes Shannon - - PowerPoint PPT Presentation

▶
dyntxregime
SMART_READER_LITE
LIVE PREVIEW

DynTxRegime An R Package for Dynamic Treatment Regimes Shannon - - PowerPoint PPT Presentation

Feb 24, 2023 121 likes •1.02k views

DynTxRegime An R Package for Dynamic Treatment Regimes Shannon Holloway Department of Statistics; NCSU September 10, 2020 Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 1 / 88 Background Background

slide-1
SLIDE 1

DynTxRegime

An R Package for Dynamic Treatment Regimes Shannon Holloway

Department of Statistics; NCSU

September 10, 2020

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 1 / 88

slide-2
SLIDE 2

Background

Background

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 2 / 88

slide-3
SLIDE 3

Background

Acknowledgements

Innovative Methods Program for Advancing Clinical Trials A joint venture of Duke, UNC-Chapel Hill, and NC State Supported by NCI Program Project P01 CA142538 (2010-2021) http://impact.unc.edu Statistical methods for precision cancer medicine

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 3 / 88

slide-4
SLIDE 4

Background

Goal

DynTxRegime was envisioned to be a toolkit for all things related to estimating DTRs The guiding principles to its development are User-friendly straightforward, consistent inputs and post-processing General minimize artificial limitation of code choices Expandable design task oriented design (S4 dominated structure)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 4 / 88

slide-5
SLIDE 5

Background

Goal

DynTxRegime was envisioned to be a toolkit for all things related to estimating DTRs The guiding principles to its development are User-friendly straightforward, consistent inputs and post-processing General minimize artificial limitation of code choices Expandable design task oriented design (S4 dominated structure) Efficient? generalization first priority

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 5 / 88

slide-6
SLIDE 6

Background

Key Challenge

How to allow for general specification of regression steps?

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 6 / 88

slide-7
SLIDE 7

Background

Key Challenge

How to allow for general specification of regression steps?

Our solution – the modelObj package. We start here Do you have the package installed?

library(package = modelObj)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 7 / 88

slide-8
SLIDE 8

modelObj

modelObj

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 8 / 88

slide-9
SLIDE 9

modelObj

What is it?

modelObj implements the design of a modeling object “modelObj” is a Class with state variables and behaviors The state variables define the model, method for obtaining parameter estimates, and method for obtaining predictions The only behavior is ’obtain parameter estimates’

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 9 / 88

slide-10
SLIDE 10

modelObj

How does it help?

For developers, this paradigm separates the details of a regression step from its use in a method For users, this scheme encapsulates each regression step into a compact object Users interface with package modelObj through function buildModelObj() buildModelObj() is used to specify each regression step of the methods implemented in DynTxRegime

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 10 / 88

slide-11
SLIDE 11

buildModelObj()

buildModelObj()

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 11 / 88

slide-12
SLIDE 12

buildModelObj()

str(object = buildModelObj) function (model, solver.method = NULL, solver.args = NULL, predict.method = NULL, predict.args = NULL)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 12 / 88

slide-13
SLIDE 13

buildModelObj()

str(object = buildModelObj) function (model, solver.method = NULL, solver.args = NULL, predict.method = NULL, predict.args = NULL) model: “formula”; the model Standard R “formula” object y ~ x Left-hand-side variable ignored; can (should) be omitted

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 13 / 88

slide-14
SLIDE 14

buildModelObj()

str(object = buildModelObj) function (model, solver.method = NULL, solver.args = NULL, predict.method = NULL, predict.args = NULL) solver.method: “character”; name of the R function to be used to obtain parameter estimates ‘lm’, ‘glm’, ‘nls’ solver.args: “list”; additional arguments for solver.method

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 14 / 88

slide-15
SLIDE 15

buildModelObj()

str(object = buildModelObj) function (model, solver.method = NULL, solver.args = NULL, predict.method = NULL, predict.args = NULL) solver.method: “character”; name of the R function to be used to obtain parameter estimates ‘lm’, ‘glm’, ‘nls’ solver.args: “list”; additional arguments for solver.method str(object = glm) function (formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart,

  • ffset, control = list(...), model = TRUE, method = "glm.fit",

x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, ...) solver.args = list(family = “binomial”)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 15 / 88

slide-16
SLIDE 16

buildModelObj()

str(object = buildModelObj) function (model, solver.method = NULL, solver.args = NULL, predict.method = NULL, predict.args = NULL) predict.method: “character”; name of the R function to be used to obtain predictions ‘predict.lm’, ‘predict.glm’, ‘predict’ predict.args: “list”; additional arguments for predict.method

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 16 / 88

slide-17
SLIDE 17

buildModelObj()

str(object = buildModelObj) function (model, solver.method = NULL, solver.args = NULL, predict.method = NULL, predict.args = NULL) predict.method: “character”; name of the R function to be used to obtain predictions ‘predict.lm’, ‘predict.glm’, ‘predict’ predict.args: “list”; additional arguments for predict.method str(object = predict.glm) function (object, newdata = NULL, type = c("link", "response", "terms"), se.fit = FALSE, dispersion = NULL, terms = NULL, na.action = na.pass, ...) predict.args = list(type = “response”)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 17 / 88

slide-18
SLIDE 18

buildModelObj()

Example 1

Assume a dataset contains covariates (x1, x2, x3) and a continuous response variable y First - specify state variables Postulate model: y ∼ b0 + b1x1 + b2x2 + b3x3 Obtain parameter estimates using ordinary least square: lm() Obtain predictions: predict.lm() (dispatches predict()) Use default settings of fitting function and prediction method

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 18 / 88

slide-19
SLIDE 19

buildModelObj()

Example 1

Assume a dataset contains covariates (x1, x2, x3) and a continuous response variable y First - specify state variables Postulate model: y ∼ b0 + b1x1 + b2x2 + b3x3 Obtain parameter estimates using ordinary least square: lm() Obtain predictions: predict.lm() (dispatches predict()) Use default settings of fitting function and prediction method Second - define model object buildModelObj(mode = ~x1+x2+x3, solver.method = "lm", predict.method = "predict.lm")

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 19 / 88

slide-20
SLIDE 20

buildModelObj()

Example 2

Assume a dataset contains covariates (x1, x2, x3) and a binary response variable y First - specify state variables Postulate model: logit(y) ∼ b0 + b1x1 + b2x2 + b3x3 Obtain parameter estimates using ordinary least square: glm() Obtain predictions: predict.glm() (dispatches predict()) Use default settings of fitting function and prediction method

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 20 / 88

slide-21
SLIDE 21

buildModelObj()

Example 2

Assume a dataset contains covariates (x1, x2, x3) and a binary response variable y First - specify state variables Postulate model: logit(y) ∼ b0 + b1x1 + b2x2 + b3x3 Obtain parameter estimates using ordinary least square: glm() Obtain predictions: predict.glm() (dispatches predict()) Use default settings of fitting function and prediction method Second - define model object buildModelObj(mode = ~x1+x2+x3, solver.method = "glm", solver.args = list(family = "binomial"), predict.method = "predict.glm")

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 21 / 88

slide-22
SLIDE 22

buildModelObj()

There’s a lot more . . .

We’ve glossed over a few details. For the applications in this class, they are not important. For the record: solver.method must have a corresponding predict.method It does not have to have coef(), residuals(), plot(), etc. Function specified in solver.method must take

A “formula” object and a “data.frame” object; OR A design matrix (X) and a response vector (Y)

Function specified by predict.method must take the regression value object and a “data.frame” object Can define your own regression and predictions methods.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 22 / 88

slide-23
SLIDE 23

Requirements for modelObj in DynTxRegime

Requirements for modelObj in DynTxRegime

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 23 / 88

slide-24
SLIDE 24

Requirements for modelObj in DynTxRegime

Propensity Regression

“modelObj” must specify that predictions are returned on the scale of the probability, i.e., in the interval (0,1) For glm()/predict.glm(), default is scale of linear predictors str(object = predict.glm) function (object, newdata = NULL, type = c("link", "response", "terms"), se.fit = FALSE, dispersion = NULL, terms = NULL, na.action = na.pass, ...) buildModelObj(mode = ~x1+x2+x3, solver.method = "glm", solver.args = list(family = "binomial"), predict.method = "predict.glm", predict.args = list(type = "response"))

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 24 / 88

slide-25
SLIDE 25

Requirements for modelObj in DynTxRegime

Outcome Regression

“modelObj” must specify that predictions are returned on the scale of the response For lm()/predict.lm(), default is scale of the response str(object = predict.lm) function (object, newdata, se.fit = FALSE, scale = NULL, df = Inf, interval = c("none", "confidence", "prediction"), level = 0.95, type = c("response", "terms"), terms = NULL, na.action = na.pass, pred.var = res.var/weights, weights = 1, ...)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 25 / 88

slide-26
SLIDE 26

Requirements for modelObj in DynTxRegime

Outcome Regression

Specified using two “modelObj” objects y ∼ µ + A C Specify µ as a “modelObj” object and/or C as a “modelObj” object Can specify different model structures for each component, e.g., a linear model for µ and a non-linear model for C µ and C are combined into a single model if same model structure or fit separately if different model structures (iterative algorithm)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 26 / 88

slide-27
SLIDE 27

Requirements for modelObj in DynTxRegime

Responsible Model Selection

WARNING!!

DynTxRegime and modelObj do NOT check for appropriateness of models Responsible model selection is the responsibility of the user

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 27 / 88

slide-28
SLIDE 28

Toy Dataset

Toy Dataset

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 28 / 88

slide-29
SLIDE 29

Toy Dataset

https::www4.stat.ncsu.edu/~sthollow/ST790/ Download st790Data.txt Load Data df <- read.csv(file = "path/st790Data.txt", header = TRUE) Summary Statistics

summary(object = df) x1 x2 x3 Min. :-3.40000 Min. :0.000 Min. : 0.640 1st Qu.:-0.67000 1st Qu.:0.000 1st Qu.: 7.997 Median :-0.04000 Median :0.000 Median :10.120 Mean :-0.02656 Mean :0.306 Mean :10.072 3rd Qu.: 0.62000 3rd Qu.:1.000 3rd Qu.:12.105 Max. : 3.20000 Max. :1.000 Max. :19.500 A y Min. :0.000 Min. :-5.0600 1st Qu.:0.000 1st Qu.:-1.3100 Median :0.000 Median :-0.3000 Mean :0.499 Mean :-0.3138 3rd Qu.:1.000 3rd Qu.: 0.7025 Max. :1.000 Max. : 3.9900

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 29 / 88

slide-30
SLIDE 30

DynTxRegime

Section 1 DynTxRegime

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 30 / 88

slide-31
SLIDE 31

DynTxRegime

Ready?

Do you have the package installed? library(package = DynTxRegime)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 31 / 88

slide-32
SLIDE 32

DynTxRegime Outcome Regression (Q-Learning)

Outcome Regression (Q-Learning)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 32 / 88

slide-33
SLIDE 33

DynTxRegime Outcome Regression (Q-Learning)

str(object = qLearn) function (..., moMain, moCont, data, response, txName, fSet = NULL, iter = 0L, verbose = TRUE) argument class description ...

  • Ignored. Included to require named input.

moMain "modelObj" mu of outcome regression. moCont "modelObj" C of outcome regression. data "data.frame" covariates and treatment histories response "vector"

  • utcome of interest

txName "character" treatment variable name fSet not used in single stage iter "numeric" if >0, iterative methods verbose "logical" if FALSE, screen prints suppressed.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 33 / 88

slide-34
SLIDE 34

DynTxRegime Outcome Regression (Q-Learning)

str(object = qLearn) function (..., moMain, moCont, data, response, txName, fSet = NULL, iter = 0L, verbose = TRUE) Postulate a model for Q1(h1, a1). Recall, ∼ µ + AC Y ∼

µ

  • β0 + β1x1 + β2x2 +A β3 + β4x2 + β5x3
  • C

Use least squares with a single model Try to create the “modelObj”. What should the value of iter be?

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 34 / 88

slide-35
SLIDE 35

DynTxRegime Outcome Regression (Q-Learning)

str(object = qLearn) function (..., moMain, moCont, data, response, txName, fSet = NULL, iter = 0L, verbose = TRUE) Postulate a model for Q1(h1, a1). Recall, ∼ µ + AC Y ∼

µ

  • β0 + β1x1 + β2x2 +A β3 + β4x2 + β5x3
  • C

Use least squares with a single model moMain <- buildModelObj(model = ~x1+x2, solver.method = 'lm', predict.method = 'predict.lm') moCont <- buildModelObj(mode = ~x2+x3, solver.method = 'lm', predict.method = 'predict.lm') iter <- 0L

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 35 / 88

slide-36
SLIDE 36

DynTxRegime Outcome Regression (Q-Learning)

str(object = qLearn) function (..., moMain, moCont, data, response, txName, fSet = NULL, iter = 0L, verbose = TRUE) data is a complete dataset, (df) response is the “vector” outcome of interest (df$y) txName is the treatment variable name, (“A”)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 36 / 88

slide-37
SLIDE 37

DynTxRegime Outcome Regression (Q-Learning)

moMain <- buildModelObj(model = ~x1+x2, solver.method = 'lm', predict.method = 'predict.lm') moCont <- buildModelObj(mode = ~x2+x3, solver.method = 'lm', predict.method = 'predict.lm') qObj <- qLearn(moMain = moMain, moCont = moCont, iter = 0L, data = df, response = df$y, txName = 'A', verbose = TRUE)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 37 / 88

slide-38
SLIDE 38

DynTxRegime Outcome Regression (Q-Learning)

qObj <- qLearn(moMain = moMain, moCont = moCont, iter = 0L, data = df, response = df$y, txName = 'A', verbose = TRUE) First step of the Q-Learning Algorithm. Outcome regression. Combined outcome regression model: ~ x1+x2 + A + A:(x2+x3) . Regression analysis for Combined: Call: lm(formula = YinternalY ~ x1 + x2 + A + x2:A + A:x3, data = data) Coefficients: (Intercept) x1 x2 A 0.08323 0.74739 0.21868 0.99709 x2:A A:x3 0.31114

  • 0.19662

Recommended Treatments: 1 911 89 Estimated value: 0.1537753

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 38 / 88

slide-39
SLIDE 39

DynTxRegime Outcome Regression (Q-Learning)

Methods available: method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by regression analysis plot(x, suppress) plot fit results summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 39 / 88

slide-40
SLIDE 40

DynTxRegime Outcome Regression (Q-Learning)

Methods available: Model Diagnostics method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by regression analysis plot(x, suppress) plot fit results summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 40 / 88

slide-41
SLIDE 41

DynTxRegime Outcome Regression (Q-Learning)

coef(object = qObj) $outcome $outcome$Combined (Intercept) x1 x2 A x2:A 0.08322797 0.74739403 0.21868299 0.99708958 0.31113694 A:x3

  • 0.19662235

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 41 / 88

slide-42
SLIDE 42

DynTxRegime Outcome Regression (Q-Learning) fitObj <- fitObject(object = qObj) fitObj $outcome $outcome$Combined Call: lm(formula = YinternalY ~ x1 + x2 + A + x2:A + A:x3, data = data) Coefficients: (Intercept) x1 x2 A 0.08323 0.74739 0.21868 0.99709 x2:A A:x3 0.31114

  • 0.19662

is(object = fitObj$outcome$Combined) [1] "lm" "oldClass" utils::methods(class = is(object = fitObj$outcome$Combined)[1L]) [1] add1 alias anova [4] case.names coerce confint [7] cooks.distance deviance dfbeta [10] dfbetas drop1 dummy.coef [13] effects extractAIC family [16] formula hatvalues influence [19] initialize kappa labels [22] logLik model.frame model.matrix [25] nobs plot predict [28] print proj qr [31] residuals rstandard rstudent [34] show simulate slotsFromS3 [37] summary variable.names vcov see '?methods' for accessing help and source code Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 42 / 88

slide-43
SLIDE 43

DynTxRegime Outcome Regression (Q-Learning)

Methods available: Training Diagnostics method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by regression analysis plot(x, suppress) plot fit results summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 43 / 88

slide-44
SLIDE 44

DynTxRegime Outcome Regression (Q-Learning)

  • t <- optTx(x = qObj)

names(ot) [1] "optimalTx" "decisionFunc" table(ot$optimalTx) 1 911 89 head(ot$decisionFunc) 1 [1,] -0.8211188 -2.8480810 [2,] 0.2924983 0.1275498 [3,] 0.8904135 0.3381190 [4,] -1.6731480 -2.3139226 [5,] 0.4046074 -0.3816677 [6,] 0.4643989 -0.2825517 estimator(x = qObj) [1] 0.1537753

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 44 / 88

slide-45
SLIDE 45

DynTxRegime Outcome Regression (Q-Learning)

Methods available: Predictions method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by regression analysis plot(x, suppress) plot fit results summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 45 / 88

slide-46
SLIDE 46

DynTxRegime Outcome Regression (Q-Learning)

newPatient <- data.frame(x1 = 2, x2 = 0, x3 = 8)

  • ptTx(x = qObj, newdata = newPatient)

$optimalTx [1] 0 $decisionFunc 1 [1,] 1.578016 1.002127

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 46 / 88

slide-47
SLIDE 47

DynTxRegime Value Search

Value Search

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 47 / 88

slide-48
SLIDE 48

DynTxRegime Value Search

/

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 48 / 88

slide-49
SLIDE 49

DynTxRegime Value Search

str(object = optimalSeq) function (..., moPropen, moMain, moCont, data, response, txName, regimes, fSet = NULL, refit = FALSE, iter = 0L, verbose = TRUE) argument class description ... Additional inputs for rgenoud() moPropen "modelObj" propensity regression moMain "modelObj" mu of outcome regression. moCont "modelObj" C of outcome regression. data "data.frame" covariates and treatment histories response "vector"

  • utcome of interest

txName "character" treatment variable name regimes "function" definition of restricted class of regimes fSet not used in single stage refit deprecated iter "numeric" if >0, iterative methods verbose "logical" if FALSE, screen prints suppressed.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 49 / 88

slide-50
SLIDE 50

DynTxRegime Value Search

str(object = optimalSeq) function (..., moPropen, moMain, moCont, data, response, txName, regimes, fSet = NULL, refit = FALSE, iter = 0L, verbose = TRUE) Postulate a model for π1(h1, a1). logit(A) ∼ γ1 Use maximum likelihood. Try to create the “modelObj”.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 50 / 88

slide-51
SLIDE 51

DynTxRegime Value Search

str(object = optimalSeq) function (..., moPropen, moMain, moCont, data, response, txName, regimes, fSet = NULL, refit = FALSE, iter = 0L, verbose = TRUE) Postulate a model for π1(h1, a1). logit(A) ∼ γ1 Use maximum likelihood. moPropen <- buildModelObj(model = ~1, solver.method = 'glm', solver.args = list(family = 'binomial'), predict.method = 'predict.glm', predict.args = list(type = "response"))

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 51 / 88

slide-52
SLIDE 52

DynTxRegime Value Search

str(object = optimalSeq) function (..., moPropen, moMain, moCont, data, response, txName, regimes, fSet = NULL, refit = FALSE, iter = 0L, verbose = TRUE) Postulate a model for Q1(h1, a1). Recall, ∼ µ + AC Y ∼

µ

  • β0 + β1x1 + β2x2 +A β3 + β4x2 + β5x3
  • C

Use least squares with a single model We’ll use our previously defined model objects.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 52 / 88

slide-53
SLIDE 53

DynTxRegime Value Search

str(object = optimalSeq) function (..., moPropen, moMain, moCont, data, response, txName, regimes, fSet = NULL, refit = FALSE, iter = 0L, verbose = TRUE) data is a complete dataset, (df) response is the “vector” outcome of interest (df$y) txName is the treatment variable name, (“A”)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 53 / 88

slide-54
SLIDE 54

DynTxRegime Value Search

str(object = optimalSeq) function (..., moPropen, moMain, moCont, data, response, txName, regimes, fSet = NULL, refit = FALSE, iter = 0L, verbose = TRUE) Define the elements of the restricted class of regimes as d1(h1, η) = I(x1 < η1 & x3 < η2) regimes <- function(eta1, eta2, data){ d1 <- {data$x1 < eta1} & {data$x3 < eta2} return(as.integer(x = d1)) }

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 54 / 88

slide-55
SLIDE 55

DynTxRegime Value Search

str(object = optimalSeq) function (..., moPropen, moMain, moCont, data, response, txName, regimes, fSet = NULL, refit = FALSE, iter = 0L, verbose = TRUE) Uses genetic algorithm implemented by genoud() Additional inputs forgenoud() passed through ellipsis. Must include starting.values initial estimates for regime parameters Domains matrix defining search space pop.size population size starting.values <- c(0,0) Domains <- matrix(data = c(-10,-10,10,10), ncol = 2L) pop.size <- 500 #TOO SMALL

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 55 / 88

slide-56
SLIDE 56

DynTxRegime Value Search

moPropen <- buildModelObj(model = ~ 1, solver.method = 'glm', solver.args = list(family = 'binomial'), predict.method = 'predict.glm', predict.args = list(type = 'response')) moMain <- buildModelObj(model = ~ x1 + x2, solver.method = 'lm') moCont <- buildModelObj(model = ~ x2 + x3, solver.method = 'lm') regimes <- function(eta1, eta2, data){ d1 <- {data$x1 < eta1} & {data$x3 < eta2} return( as.integer(x = d1) ) } vsObj <- optimalSeq(moPropen = moPropen, moMain = moMain, moCont = moCont, iter = 0L, data = df, response = df$y, txName = 'A', regimes = regimes, Domains = matrix(data = c(-10,-10,10,10), ncol = 2L), starting.values = c(0,0), pop.size = 500, verbose = TRUE)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 56 / 88

slide-57
SLIDE 57

DynTxRegime Value Search vsObj <- optimalSeq(moPropen = moPropen, moMain = moMain, moCont = moCont, iter = 0L, data = df, response = df$y, txName = 'A', regimes = regimes, Domains = matrix(data = c(-10,-10,10,10), ncol = 2L), starting.values = c(0,0), pop.size = 500, verbose = TRUE) Value Search - Missing Data Perspective. Propensity for treatment regression. Regression analysis for moPropen: Call: glm(formula = YinternalY ~ 1, family = "binomial", data = data) Coefficients: (Intercept)

  • 0.004

Degrees of Freedom: 999 Total (i.e. Null); 999 Residual Null Deviance: 1386 Residual Deviance: 1386 AIC: 1388 Outcome regression. Combined outcome regression model: ~ x1+x2 + A + A:(x2+x3) . Regression analysis for Combined: Call: lm(formula = YinternalY ~ x1 + x2 + A + x2:A + A:x3, data = data) Coefficients: (Intercept) x1 x2 A 0.08323 0.74739 0.21868 0.99709 x2:A A:x3 0.31114

  • 0.19662

Wed Sep 9 17:34:34 2020 Domains:

  • 1.000000e+01

<= X1 <= 1.000000e+01

  • 1.000000e+01

<= X2 <= 1.000000e+01 Data Type: Floating Point Operators (code number, name, population) Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 57 / 88

slide-58
SLIDE 58

DynTxRegime Value Search

Methods available: method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s genetic(object) retrieve value returned by genoud()

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis regimeCoef(object) retrieve estimated regime parameters summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 58 / 88

slide-59
SLIDE 59

DynTxRegime Value Search

Methods available: Model Diagnostics method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s genetic(object) retrieve value returned by genoud()

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis regimeCoef(object) retrieve estimated regime parameters summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 59 / 88

slide-60
SLIDE 60

DynTxRegime Value Search

coef(object = vsObj) $propensity (Intercept)

  • 0.004000005

$outcome $outcome$Combined (Intercept) x1 x2 A x2:A 0.08322797 0.74739403 0.21868299 0.99708958 0.31113694 A:x3

  • 0.19662235

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 60 / 88

slide-61
SLIDE 61

DynTxRegime Value Search fitObj <- fitObject(object = vsObj) print(x = fitObj) $propensity Call: glm(formula = YinternalY ~ 1, family = "binomial", data = data) Coefficients: (Intercept)

  • 0.004

Degrees of Freedom: 999 Total (i.e. Null); 999 Residual Null Deviance: 1386 Residual Deviance: 1386 AIC: 1388 $outcome $outcome$Combined Call: lm(formula = YinternalY ~ x1 + x2 + A + x2:A + A:x3, data = data) Coefficients: (Intercept) x1 x2 A 0.08323 0.74739 0.21868 0.99709 x2:A A:x3 0.31114

  • 0.19662

is(object = fitObj$propensity) [1] "glm" "lm" "oldClass" is(object = fitObj$outcome$Combined) [1] "lm" "oldClass" Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 61 / 88

slide-62
SLIDE 62

DynTxRegime Value Search genetic(object = vsObj) $value [1] 0.1534204 $par [1] 0.4259397 4.6884133 $gradients [1] NA NA $generations [1] 12 $peakgeneration [1] 1 $popsize [1] 500 $operators [1] 65 62 62 62 62 62 62 62 Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 62 / 88

slide-63
SLIDE 63

DynTxRegime Value Search

Methods available: Training Diagnostics method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s genetic(object) retrieve value returned by genoud()

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis regimeCoef(object) retrieve estimated regime parameters summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 63 / 88

slide-64
SLIDE 64

DynTxRegime Value Search

  • t <- optTx(x = vsObj)

names(ot) [1] "optimalTx" "decisionFunc" table(ot$optimalTx) 1 970 30

  • t$decisionFunc

[1] NA estimator(x = vsObj) [1] 0.1534204 regimeCoef(object = vsObj) eta1 eta2 0.4259397 4.6884133

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 64 / 88

slide-65
SLIDE 65

DynTxRegime Value Search

Methods available: Predictions method description Call(name) retrieve the unevaluated call coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s genetic(object) retrieve value returned by genoud()

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis regimeCoef(object) retrieve estimated regime parameters summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 65 / 88

slide-66
SLIDE 66

DynTxRegime Value Search

newPatient <- data.frame(x1 = 2, x2 = 0, x3 = 8)

  • ptTx(x = vsObj, newdata = newPatient)

$optimalTx [1] 0 $decisionFunc [1] NA

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 66 / 88

slide-67
SLIDE 67

DynTxRegime Classification

Classification

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 67 / 88

slide-68
SLIDE 68

DynTxRegime Classification

/

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 68 / 88

slide-69
SLIDE 69

DynTxRegime Classification

str(object = optimalClass) function (..., moPropen, moMain, moCont, moClass, data, response, txName, iter = 0L, fSet = NULL, verbose = TRUE) argument class description ...

  • Ignored. Included to require named input.

moPropen "modelObj" propensity regression moMain "modelObj" m of outcome regression. moCont "modelObj" C of outcome regression. moClass "modelObj" classification regression data "data.frame" covariates and treatment histories response "vector"

  • utcome of interest

txName "character" treatment variable name iter "numeric" if >0, iterative methods fSet not used in single stage verbose "logical" if FALSE, screen prints suppressed.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 69 / 88

slide-70
SLIDE 70

DynTxRegime Classification

str(object = optimalClass) function (..., moPropen, moMain, moCont, moClass, data, response, txName, iter = 0L, fSet = NULL, verbose = TRUE) Postulate a model for π1(h1, a1). logit(A) ∼ γ1 Use maximum likelihood. We’ll use the same "modelObj" previously defined.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 70 / 88

slide-71
SLIDE 71

DynTxRegime Classification

str(object = optimalClass) function (..., moPropen, moMain, moCont, moClass, data, response, txName, iter = 0L, fSet = NULL, verbose = TRUE) Postulate a model for Q1(h1, a1). Recall, ∼ µ + AC Y ∼

µ

  • β0 + β1x1 + β2x2 +A β3 + β4x2 + β5x3
  • C

Use least squares with a single model We’ll use the same "modelObj" previously defined.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 71 / 88

slide-72
SLIDE 72

DynTxRegime Classification

str(object = optimalClass) function (..., moPropen, moMain, moCont, moClass, data, response, txName, iter = 0L, fSet = NULL, verbose = TRUE) Define the classification method. Method must use input weights and return predictions as the class. library(rpart) moClass <- buildModelObj(model = ~x1+x2+x3, solver.method = 'rpart', predict.method = 'predict', predict.args = list(type = "class"))

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 72 / 88

slide-73
SLIDE 73

DynTxRegime Classification

str(object = optimalClass) function (..., moPropen, moMain, moCont, moClass, data, response, txName, iter = 0L, fSet = NULL, verbose = TRUE) data is a complete dataset, (df) response is the “vector” outcome of interest (df$y) txName is the treatment variable name, (“A”)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 73 / 88

slide-74
SLIDE 74

DynTxRegime Classification

moPropen <- buildModelObj(model = ~ 1, solver.method = 'glm', solver.args = list(family = 'binomial'), predict.method = 'predict.glm', predict.args = list(type = 'response')) moMain <- buildModelObj(model = ~ x1 + x2, solver.method = 'lm') moCont <- buildModelObj(model = ~ x2 + x3, solver.method = 'lm') moClass <- buildModelObj(model = ~x1+x2+x3, solver.method = 'rpart', predict.method = 'predict', predict.args = list(type = "class")) clObj <- optimalClass(moPropen = moPropen, moMain = moMain, moCont = moCont, iter = 0L, moClass = moClass, data = df, response = df$y, txName = 'A', verbose = TRUE)

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 74 / 88

slide-75
SLIDE 75

DynTxRegime Classification clObj <- optimalClass(moPropen = moPropen, moMain = moMain, moCont = moCont, iter = 0L, moClass = moClass, data = df, response = df$y, txName = 'A', verbose = TRUE) AIPW value estimator First step of the Classification Algorithm. Classification Perspective. Propensity for treatment regression. Regression analysis for moPropen: Call: glm(formula = YinternalY ~ 1, family = "binomial", data = data) Coefficients: (Intercept)

  • 0.004

Degrees of Freedom: 999 Total (i.e. Null); 999 Residual Null Deviance: 1386 Residual Deviance: 1386 AIC: 1388 Outcome regression. Combined outcome regression model: ~ x1+x2 + A + A:(x2+x3) . Regression analysis for Combined: Call: lm(formula = YinternalY ~ x1 + x2 + A + x2:A + A:x3, data = data) Coefficients: (Intercept) x1 x2 A 0.08323 0.74739 0.21868 0.99709 x2:A A:x3 0.31114

  • 0.19662

Classification Analysis Regression analysis for moClass: n= 1000 Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 75 / 88

slide-76
SLIDE 76

DynTxRegime Classification

Methods available: method description Call(name) retrieve the unevaluated call classif(object) retrieve value returned by classification method coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 76 / 88

slide-77
SLIDE 77

DynTxRegime Classification

Methods available: Model Diagnostics method description Call(name) retrieve the unevaluated call classif(object) retrieve value returned by classification method coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 77 / 88

slide-78
SLIDE 78

DynTxRegime Classification

coef(object = clObj) $propensity (Intercept)

  • 0.004000005

$outcome $outcome$Combined (Intercept) x1 x2 A x2:A 0.08322797 0.74739403 0.21868299 0.99708958 0.31113694 A:x3

  • 0.19662235

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 78 / 88

slide-79
SLIDE 79

DynTxRegime Classification

classObj <- classif(object = clObj) classObj n= 1000 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 1000 0.254446300 0 (0.89567578 0.10432422) 2) x3>=8.48 701 0.134399100 0 (0.93011385 0.06988615) * 3) x3< 8.48 299 0.120047100 0 (0.76729621 0.23270379) 6) x3< 8.32 278 0.107693200 0 (0.78016306 0.21983694) 12) x3>=3.985 254 0.094892290 0 (0.79479473 0.20520527) * 13) x3< 3.985 24 0.005000009 1 (0.53368911 0.46631089) * 7) x3>=8.32 21 0.004657874 1 (0.52488472 0.47511528) 14) x1>=0.135 7 0.001631934 0 (0.88818501 0.11181499) * 15) x1< 0.135 14 0.000233792 1 (0.06005364 0.93994636) * is(object = classObj) [1] "rpart" utils::methods(class = is(object = classObj)[1L]) [1] labels meanvar model.frame plot [5] post predict print prune [9] residuals summary text see '?methods' for accessing help and source code

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 79 / 88

slide-80
SLIDE 80

DynTxRegime Classification

Methods available: Training Diagnostics method description Call(name) retrieve the unevaluated call classif(object) retrieve value returned by classification method coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 80 / 88

slide-81
SLIDE 81

DynTxRegime Classification

  • t <- optTx(x = clObj)

names(ot) [1] "optimalTx" "decisionFunc" table(ot$optimalTx) 1 962 38

  • t$decisionFunc

[1] NA estimator(x = clObj) [1] 0.1633692

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 81 / 88

slide-82
SLIDE 82

DynTxRegime Classification

Methods available: Predictions method description Call(name) retrieve the unevaluated call classif(object) retrieve value returned by classification method coef(object) retrieve parameter estimates DTRstep(object) print description of method estimator(x) retrieve estimated value fitObject(object) retrieve value returned by ’modelObj’s

  • ptTx(object)

retrieve recommended optimal treatments

  • ptTx(object, newdata)

estimate optimal treatments based on a prior analysis

  • utcome(object)

retrieve value returned by outcome regression analysis plot(x, suppress) plot fit results propen(object) retrieve value returned by propensity regression analysis summary(object) retrieve summary information

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 82 / 88

slide-83
SLIDE 83

DynTxRegime Classification

newPatient <- data.frame(x1 = 2, x2 = 0, x3 = 8)

  • ptTx(x = clObj, newdata = newPatient)

$optimalTx [1] 0 $decisionFunc [1] NA

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 83 / 88

slide-84
SLIDE 84

DynTxRegime Outcome Weighted Learning

Outcome Weighted Learning

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 84 / 88

slide-85
SLIDE 85

DynTxRegime Outcome Weighted Learning

/

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 85 / 88

slide-86
SLIDE 86

DynTxRegime Outcome Weighted Learning

str(object = owl) function (..., moPropen, data, reward, txName, regime, response, lambdas = 2, cvFolds = 0L, kernel = "linear", kparam = NULL, surrogate = "hinge", verbose = 2L) argument class description ...

  • Ignored. Included to require named input.

moPropen "modelObj" propensity regression data "data.frame" covariates and treatment histories reward "vector"

  • utcome of interest

txName "character" treatment variable name regime "formula" covariates of the decision function response "numeric"

  • utcome of interest

lambdas "numeric"

  • ne or more tuning parameters

cvFolds "integer" number of cross-validation steps kernel "character"

  • ne of {linear, poly, radial}

kparam "numeric" kernel parameter surrogate if >0, iterative methods verbose "logical" if FALSE, screen prints suppressed.

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 86 / 88

slide-87
SLIDE 87

DynTxRegime Conclusion

Conclusion

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 87 / 88

slide-88
SLIDE 88

DynTxRegime Conclusion

Packages are always under development. For a limited time, these slides are available from https://www4.stat.ncsu.edu/~sthollow/ST790/ If you have ANY questions or encounter any problems, please do not hesitate to contact me Shannon Holloway sthollow@ncsu.edu

Shannon Holloway (Department of Statistics; NCSU) DynTxRegime September 10, 2020 88 / 88