Fits subgroup identification model class of Chen, et al (2017)

fit.subgroup(x, y, trt, propensity.func = NULL, loss = c("sq_loss_lasso",
"logistic_loss_lasso", "cox_loss_lasso", "owl_logistic_loss_lasso",
"owl_logistic_flip_loss_lasso", "owl_hinge_loss", "owl_hinge_flip_loss",
"sq_loss_lasso_gam", "logistic_loss_lasso_gam", "sq_loss_gam",
"logistic_loss_gam", "owl_logistic_loss_gam", "owl_logistic_flip_loss_gam",
"owl_logistic_loss_lasso_gam", "owl_logistic_flip_loss_lasso_gam",
"sq_loss_gbm", "abs_loss_gbm", "logistic_loss_gbm", "cox_loss_gbm"),
method = c("weighting", "a_learning"), match.id = NULL,
augment.func = NULL, cutpoint = 0, larger.outcome.better = TRUE,
reference.trt = NULL, retcall = TRUE, ...)

## Arguments

x The design matrix (not including intercept term) The response vector treatment vector with each element equal to a 0 or a 1, with 1 indicating treatment status is active. function that inputs the design matrix x and the treatment vector trt and outputs the propensity score, ie Pr(trt = 1 | X = x). Function should take two arguments 1) x and 2) trt. See example below. For a randomized controlled trial this can simply be a function that returns a constant equal to the proportion of patients assigned to the treatment group, i.e.: propensity.func = function(x, trt) 0.5. choice of both the M function from Chen, et al (2017) and potentially the penalty used for variable selection. All loss options starting with sq_loss use M(y, v) = (v - y) ^ 2, all options starting with logistic_loss use the logistic loss: M(y, v) = y * log(1 + exp{-v}), and all options starting with cox_loss use the negative partial likelihood loss for the Cox PH model. All options ending with lasso have a lasso penalty added to the loss for variable selection. sq_loss_lasso_gam and logistic_loss_lasso_gam first use the lasso to select variables and then fit a generalized additive model with nonparametric additive terms for each selected variable. sq_loss_gam involves a squared error loss with a generalized additive model and no variable selection. sq_loss_gbm involves a squared error loss with a gradient-boosted decision trees model for the benefit score; this allows for flexible estimation using machine learning and can be useful when the underlying treatment-covariate interaction is complex. Continuous Outcomes "sq_loss_lasso" - M(y, v) = (v - y) ^ 2 with linear model and lasso penalty "owl_logistic_loss_lasso" - M(y, v) = ylog(1 + exp{-v}) (method of Regularized Outcome Weighted Subgroup Identification) "owl_logistic_flip_loss_lasso" - M(y, v) = |y|log(1 + exp{-sign(y)v}) "owl_hinge_loss" - M(y, v) = ymax(0, 1 - v) (method of Estimating individualized treatment rules using outcome weighted learning) "owl_hinge_flip_loss" - M(y, v) = |y|max(0, 1 - sign(y)v) "sq_loss_lasso_gam" - M(y, v) = (v - y) ^ 2 with variables selected by lasso penalty and generalized additive model fit on the selected variables "sq_loss_gam" - M(y, v) = (v - y) ^ 2 with generalized additive model fit on all variables "owl_logistic_loss_gam" - M(y, v) = ylog(1 + exp{-v}) with generalized additive model fit on all variables "owl_logistic_flip_loss_gam" - M(y, v) = |y|log(1 + exp{-sign(y)v}) with generalized additive model fit on all variables "owl_logistic_loss_lasso_gam" - M(y, v) = ylog(1 + exp{-v}) with variables selected by lasso penalty and generalized additive model fit on the selected variables "owl_logistic_flip_loss_lasso_gam" - M(y, v) = |y|log(1 + exp{-sign(y)v}) with variables selected by lasso penalty and generalized additive model fit on the selected variables "sq_loss_gbm" - M(y, v) = (v - y) ^ 2 with gradient-boosted decision trees model "abs_loss_gbm" - M(y, v) = |v - y| with gradient-boosted decision trees model Binary Outcomes All losses for continuous outcomes can be used plus the following: "logistic_loss_lasso" - M(y, v) = -[yv - log(1 + exp{-v})] with with linear model and lasso penalty "logistic_loss_lasso_gam" - M(y, v) = y * log(1 + exp{-v}) with variables selected by lasso penalty and generalized additive model fit on the selected variables "logistic_loss_gam" - M(y, v) = y * log(1 + exp{-v}) with generalized additive model fit on all variables "logistic_loss_gbm" - M(y, v) = -[yv - log(1 + exp{-v})] with gradient-boosted decision trees model Time-to-Event Outcomes "cox_loss_lasso" - M corresponds to the negative partial likelihood of the cox model with linear model and additionally a lasso penalty "cox_loss_gbm" - M corresponds to the negative partial likelihood of the cox model with gradient-boosted decision trees model subgroup ID model type. Either the weighting or A-learning method of Chen et al, (2017) a (character, factor, or integer) vector with length equal to the number of observations in x indicating using integers or levels of a factor vector which patients are in which matched groups. Defaults to NULL and assumes the samples are not from a matched cohort. Matched case-control groups can be created using any method (propensity score matching, optimal matching, etc). If each case is matched with a control or multiple controls, this would indicate which case-control pairs or groups go together. If match.id is supplied, then it is unecessary to specify a function via the propensity.func argument. A quick usage example: if the first patient is a case and the second and third are controls matched to it, and the fouth patient is a case and the fifth through seventh patients are matched with it, then the user should specify match.id = c(1,1,1,2,2,2,2) or match.id = c(rep("Grp1", 3),rep("Grp2", 4))  function which inputs the response y, the covariates x, and trt and outputs predicted values (on the link scale) for the response using a model constructed with x. augment.func() can also be simply a function of x and y. This function is used for efficiency augmentation. When the form of the augmentation function is correct, it can provide efficient estimation of the subgroups Example 1: augment.func <- function(x, y) {lmod <- lm(y ~ x); return(fitted(lmod))} Example 2: augment.func <- function(x, y, trt) {lmod <- lm(y ~ x + trt); return(fitted(lmod))} For binary and time-to-event outcomes, make sure that predictions are returned on the scale of the predictors Example 3: augment.func <- function(x, y) { bmod <- glm(y ~ x, family = binomial()); return(predict(bmod, type = "link"))}  numeric value for patients with benefit scores above which (or below which if larger.outcome.better = FALSE) will be recommended to be in the treatment group. Can also set cutpoint = "median", which will use the median value of the benefit scores as the cutpoint or can set specific quantile values via "quantx" where "x" is a number between 0 and 100 representing the quantile value; e.g. cutpoint = "quant75" will use the 75th perent upper quantile of the benefit scores as the quantile. boolean value of whether a larger outcome is better/preferable. Set to TRUE if a larger outcome is better/preferable and set to FALSE if a smaller outcome is better/preferable. Defaults to TRUE. which treatment should be treated as the reference treatment. Defaults to the first level of trt if trt is a factor or the first alphabetical or numerically first treatment level. Not used for multiple treatment fitting with OWL-type losses. boolean value. if TRUE then the passed arguments will be saved. Do not set to FALSE if the validate.subgroup() function will later be used for your fitted subgroup model. Only set to FALSE if memory is limited as setting to TRUE saves the design matrix to the fitted object options to be passed to underlying fitting function. For all loss options with 'lasso', this will be passed to cv.glmnet. For all loss options with 'gam', this will be passed to gam from the mgcv package Note that for all loss options that use gam() from the mgcv package, the user cannot supply the gam argument method because it is also an argument of fit.subgroup, so instead, to change the gam method argument, supply method.gam, ie method.gam = "REML". For all loss options with 'hinge', this will be passed to both weighted.ksvm and ipop from the kernlab package

## References

Chen, S., Tian, L., Cai, T. and Yu, M. (2017), A general statistical framework for subgroup identification and comparative treatment scoring. Biometrics. doi:10.1111/biom.12676 http://onlinelibrary.wiley.com/doi/10.1111/biom.12676/abstract

Xu, Y., Yu, M., Zhao, Y. Q., Li, Q., Wang, S., & Shao, J. (2015), Regularized outcome weighted subgroup identification for differential treatment effects. Biometrics, 71(3), 645-653. doi: 10.1111/biom.12322 http://onlinelibrary.wiley.com/doi/10.1111/biom.12322/full

Zhao, Y., Zeng, D., Rush, A. J., & Kosorok, M. R. (2012), Estimating individualized treatment rules using outcome weighted learning. Journal of the American Statistical Association, 107(499), 1106-1118. doi: 10.1080/01621459.2012.695674 http://dx.doi.org/10.1080/01621459.2012.695674

validate.subgroup for function which creates validation results for subgroup identification models, predict.subgroup_fitted for a prediction function for fitted models from fit.subgroup, plot.subgroup_fitted for a function which plots results from fitted models, and print.subgroup_fitted for arguments for printing options for fit.subgroup(). from fit.subgroup.

## Examples

library(personalized)

set.seed(123)
n.obs  <- 500
n.vars <- 15
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)

# simulate non-randomized treatment
xbetat   <- 0.5 + 0.5 * x[,7] - 0.5 * x[,9]
trt.prob <- exp(xbetat) / (1 + exp(xbetat))
trt01    <- rbinom(n.obs, 1, prob = trt.prob)

trt      <- 2 * trt01 - 1

# simulate response
delta <- 2 * (0.5 + x[,2] - x[,3] - x[,11] + x[,1] * x[,12] )
xbeta <- x[,1] + x[,11] - 2 * x[,12]^2 + x[,13] + 0.5 * x[,15] ^ 2
xbeta <- xbeta + delta * trt

# continuous outcomes
y <- drop(xbeta) + rnorm(n.obs, sd = 2)

# binary outcomes
y.binary <- 1 * (xbeta + rnorm(n.obs, sd = 2) > 0 )

# time-to-event outcomes
surv.time <- exp(-20 - xbeta + rnorm(n.obs, sd = 1))
cens.time <- exp(rnorm(n.obs, sd = 3))
y.time.to.event  <- pmin(surv.time, cens.time)
status           <- 1 * (surv.time <= cens.time)

# create function for fitting propensity score model
prop.func <- function(x, trt)
{
# fit propensity score model
propens.model <- cv.glmnet(y = trt,
x = x, family = "binomial")
pi.x <- predict(propens.model, s = "lambda.min",
newx = x, type = "response")[,1]
pi.x
}

subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss   = "sq_loss_lasso",
nfolds = 10)              # option for cv.glmnet

summary(subgrp.model)#> family:  gaussian
#> loss:    sq_loss_lasso
#> method:  weighting
#> propensity
#> function: propensity.func
#>
#> Average Outcomes:
#>                 Recommended 0      Recommended 1
#> Received 0   -6.2631 (n = 82) -26.8532 (n = 119)
#> Received 1 -26.4854 (n = 128)  -1.6497 (n = 171)
#>
#> 0 effect among recommended 0 1 effect among recommended 1
#>            20.2223 (n = 210)            25.2035 (n = 290)
#>
#> Benefit score quantiles:
#>      0%     25%     50%     75%    100%
#> -36.451  -6.808   2.337  10.538  37.786
#>
#> 12 out of 15 variables selected in total by the lasso (cross validation criterion).
#>
#>      Estimate
#> Trt1   2.4460
#> V1     0.1780
#> V2     2.5630
#> V3    -2.5829
#> V4    -0.1536
#> V5    -0.2795
#> V6     0.7616
#> V7    -0.3729
#> V8     0.7436
#> V10    0.8779
#> V11   -1.1277
#> V13    1.6186
#> V15    0.7653
# fit lasso + gam model with REML option for gam

subgrp.modelg <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss   = "sq_loss_lasso_gam",
method.gam = "REML",     # option for gam
nfolds = 5)              # option for cv.glmnet

subgrp.modelg#> family:  gaussian
#> loss:    sq_loss_lasso_gam
#> method:  weighting
#> propensity
#> function: propensity.func
#>
#> Average Outcomes:
#>                 Recommended 0     Recommended 1
#> Received 0 -12.6764 (n = 131) -25.1851 (n = 70)
#> Received 1 -21.0034 (n = 141)  -2.777 (n = 158)
#>
#> 0 effect among recommended 0 1 effect among recommended 1
#>              8.327 (n = 272)            22.4081 (n = 228)
#>
#> Benefit score quantiles:
#>      0%     25%     50%     75%    100%
#> -57.584 -11.433  -1.937   8.934  44.564
subgrp.model.bin <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
loss   = "logistic_loss_lasso",
type.measure = "auc",    # option for cv.glmnet
nfolds = 5)              # option for cv.glmnet

subgrp.model.bin#> family:  binomial
#> loss:    logistic_loss_lasso
#> method:  weighting
#> propensity
#> function: propensity.func
#>
#> Average Outcomes:
#>               Recommended 0    Recommended 1
#> Received 0  0.5161 (n = 85) 0.1788 (n = 116)
#> Received 1 0.1778 (n = 130) 0.5541 (n = 169)
#>
#> 0 effect among recommended 0 1 effect among recommended 1
#>             0.3383 (n = 215)             0.3753 (n = 285)
#>
#> Benefit score quantiles:
#>      0%     25%     50%     75%    100%
#> -2.0266 -0.3695  0.1284  0.6247  2.2427
library(survival)
subgrp.model.cox <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = trt01,
propensity.func = prop.func,
loss   = "cox_loss_lasso",
nfolds = 5)              # option for cv.glmnet

subgrp.model.cox#> family:  cox
#> loss:    cox_loss_lasso
#> method:  weighting
#> propensity
#> function: propensity.func
#>
#> Average Outcomes:
#>                 Recommended 0     Recommended 1
#> Received 0 115.7397 (n = 137)   2.2659 (n = 64)
#> Received 1  78.6273 (n = 196) 40.6071 (n = 103)
#>
#> 0 effect among recommended 0 1 effect among recommended 1
#>            37.1124 (n = 333)            38.3413 (n = 167)
#>
#> Benefit score quantiles:
#>      0%     25%     50%     75%    100%
#> -1.0584 -0.4067 -0.1596  0.1056  1.0928