fit.subgroup.Rd
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", "poisson_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", "poisson_loss_lasso_gam", "logistic_loss_lasso_gam", "sq_loss_gam", "poisson_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", "poisson_loss_gbm", "logistic_loss_gbm", "cox_loss_gbm", "custom"), method = c("weighting", "a_learning"), match.id = NULL, augment.func = NULL, fit.custom.loss = NULL, cutpoint = 0, larger.outcome.better = TRUE, reference.trt = NULL, retcall = TRUE, ...)
x  The design matrix (not including intercept term) 

y  The response vector 
trt  treatment vector with each element equal to a 0 or a 1, with 1 indicating treatment status is active. 
propensity.func  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.:

loss  choice of both the M function from Chen, et al (2017) and potentially the penalty used for variable selection.
All

method  subgroup ID model type. Either the weighting or Alearning method of Chen et al, (2017) 
match.id  a (character, factor, or integer) vector with length equal to the number of observations in 
augment.func  function which inputs the response Example 1: Example 2: augment.func < function(x, y, trt) { data < data.frame(x, y, trt) lmod < lm(y ~ x * trt) ## get predictions when trt = 1 data$trt < 1 preds_1 < predict(lmod, data) ## get predictions when trt = 1 data$trt < 1 preds_n1 < predict(lmod, data) ## return predictions averaged over trt return(0.5 * (preds_1 + preds_n1)) } For binary and timetoevent 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")) } 
fit.custom.loss  A function which minimizes a userspecified
custom loss function M(y,v) to be used in model fitting.
If provided, The loss function \(M(y, v)\) to be minimized MUST meet the following two criteria:
An example of a valid loss function is \(M(y, v) = (y  v)^2\). In this case \(D_M(y, v) = 2(y  v)\). See Chen et al. (2017) for more details on the restrictions on the loss function \(M(y, v)\). The provided function MUST return a list with the following elements:
The provided function MUST be a function with the following arguments:
The provided function can also optionally take the following arguments:
Example 1: Here we minimize \(M(y, v) = (y  v)^2\) fit.custom.loss < function(x, y, weights, ...) { df < data.frame(y = y, x) # minimize squared error loss with NO lasso penalty lmf < lm(y ~ x  1, weights = weights, data = df, ...) # save coefficients cfs = coef(lmf) # create prediction function. Notice # how a column of 1's is appended # to ensure treatment main effects are included # in predictions prd = function(x, type = "response") { dfte < cbind(1, x) colnames(dfte) < names(cfs) predict(lmf, data.frame(dfte)) } # return lost of required components list(predict = prd, model = lmf, coefficients = cfs) } Example 2: \(M(y, v) = y\exp(v)\) fit.expo.loss < function(x, y, weights, ...) { ## define loss function to be minimized expo.loss < function(beta, x, y, weights) { sum(weights * y * exp(drop(tcrossprod(x, t(beta) ))) } # use optim() to minimize loss function opt < optim(rep(0, NCOL(x)), fn = expo.loss, x = x, y = y, weights = weights) coefs < opt$par pred < function(x, type = "response") { tcrossprod(cbind(1, x), t(coefs)) } # return list of required components list(predict = pred, model = opt, coefficients = coefs) } 
cutpoint  numeric value for patients with benefit scores above which
(or below which if 
larger.outcome.better  boolean value of whether a larger outcome is better/preferable. Set to 
reference.trt  which treatment should be treated as the reference treatment. Defaults to the first level of 
retcall  boolean value. if 
...  options to be passed to underlying fitting function. For all For all 
An object of class "subgroup_fitted"
.
A function that returns predictions of the covariateconditional treatment effects
An object returned by the underlying fitting function used. For example, if the lasso use used to fit
the underlying subgroup identification model, this will be an object returned by cv.glmnet
.
If the underlying subgroup identification model is parametric, coefficients
will contain
the estimated coefficients of the model.
The call that produced the returned object. If retcall = TRUE
, this will contain all objects
supplied to fit.subgroup()
The family corresponding to the outcome provided
The loss function used
The method used (either weighting or Alearning)
The propensity score function used
If larger outcomes are preferred for this model
Benefit score cutoff value used for determining subgroups
The names of all variables used
The number of treatment levels
All treatment levels other than the reference level
The reference level for the treatment. This should usually be the control group/level
All treatment levels
The vector of treatment assignments
A vector of propensity scores
A vector of outcomes
A vector of conditional treatment effects, i.e. benefit scores
A vector of treatment recommendations (i.e. for each patient, which treatment results in the best expected potential outcomes)
(Biased) estimates of the conditional treatment effects and conditional outcomes. These are essentially just empirical averages within different combinations of treatment assignments and treatment recommendations
estimates of the individual treatment effects as returned by
treat.effects
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), 645653. 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), 11061118. 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
.
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 nonrandomized 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 below drives treatment effect heterogeneity 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 ) # count outcomes y.count < round(abs(xbeta + rnorm(n.obs, sd = 2))) # timetoevent 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 } #################### Continuous outcomes ################################ 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 #> cutpoint: 0 #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 5.4606 (n = 81) 27.2778 (n = 120) #> Received 1 26.8929 (n = 124) 1.6271 (n = 175) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 21.4323 (n = 205) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 25.6507 (n = 295) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 32.002 4.996 2.423 9.251 29.137 #> #>  #> #> Summary of individual treatment effects: #> E[YT=1, X]  E[YT=0, X] #> #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 64.004 9.993 4.846 4.793 18.502 58.274 #> #>  #> #> 10 out of 15 interactions selected in total by the lasso (cross validation criterion). #> #> The first estimate is the treatment main effect, which is always selected. #> Any other variables selected represent treatmentcovariate interactions. #> #> Trt1 V2 V3 V5 V6 V7 V8 V10 V11 #> Estimate 2.6463 2.2645 2.1415 0.1068 0.43 0.0028 0.5195 0.5581 0.8164 #> V13 V15 #> Estimate 1.3847 0.5343# estimates of the individualspecific # treatment effect estimates: subgrp.model$individual.trt.effects#> Summary of individual treatment effects: #> E[YT=1, X]  E[YT=0, X] #> #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 64.004 9.993 4.846 4.793 18.502 58.274# fit lasso + gam model with REML option for gam # \donttest{ 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 #> cutpoint: 0 #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 12.6754 (n = 132) 25.3863 (n = 69) #> Received 1 21.0867 (n = 141) 2.7467 (n = 158) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 8.4113 (n = 273) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 22.6397 (n = 227) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 58.022 11.497 2.069 9.130 45.392 #> #>  #> #> Summary of individual treatment effects: #> E[YT=1, X]  E[YT=0, X] #> #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 116.043 22.994 4.138 2.356 18.260 90.785# } #################### Using an augmentation function ##################### ## augmentation funcions involve modeling the conditional mean E[YT, X] ## and returning predictions that are averaged over the treatment values ## return < 1/2 * (hat{E}[YT=1, X] + hat{E}[YT=1, X]) ########################################################################## augment.func < function(x, y, trt) { data < data.frame(x, y, trt) xm < model.matrix(y~trt*x1, data = data) lmod < cv.glmnet(y = y, x = xm) ## get predictions when trt = 1 data$trt < 1 xm < model.matrix(y~trt*x1, data = data) preds_1 < predict(lmod, xm, s = "lambda.min") ## get predictions when trt = 1 data$trt < 1 xm < model.matrix(y~trt*x1, data = data) preds_n1 < predict(lmod, xm, s = "lambda.min") ## return predictions averaged over trt return(0.5 * (preds_1 + preds_n1)) } # \donttest{ subgrp.model.aug < fit.subgroup(x = x, y = y, trt = trt01, propensity.func = prop.func, augment.func = augment.func, loss = "sq_loss_lasso", nfolds = 10) # option for cv.glmnet summary(subgrp.model.aug)#> family: gaussian #> loss: sq_loss_lasso #> method: weighting #> cutpoint: 0 #> augmentation #> function: augment.func #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 8.2593 (n = 83) 23.8274 (n = 118) #> Received 1 26.798 (n = 117) 2.4579 (n = 182) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 18.5387 (n = 200) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 21.3695 (n = 300) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 33.022 5.714 2.433 10.558 37.633 #> #>  #> #> Summary of individual treatment effects: #> E[YT=1, X]  E[YT=0, X] #> #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 66.045 11.427 4.865 5.851 21.116 75.265 #> #>  #> #> 10 out of 15 interactions selected in total by the lasso (cross validation criterion). #> #> The first estimate is the treatment main effect, which is always selected. #> Any other variables selected represent treatmentcovariate interactions. #> #> Trt1 V2 V3 V5 V6 V8 V9 V10 V11 #> Estimate 3.061 2.1429 2.4738 0.0522 0.6815 0.9514 0.3491 0.6354 1.0507 #> V13 V15 #> Estimate 1.2763 1.0572# } #################### Binary outcomes #################################### # use logistic loss for binary outcomes 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 #> cutpoint: 0 #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 0.5544 (n = 84) 0.1651 (n = 117) #> Received 1 0.1801 (n = 116) 0.526 (n = 183) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 0.3743 (n = 200) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 0.361 (n = 300) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 1.4737 0.1959 0.1478 0.5169 1.6003 #> #>  #> #> Summary of individual treatment effects: #> E[YT=1, X]  E[YT=0, X] #> #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 0.62723 0.09763 0.07376 0.07294 0.25283 0.66411#################### Count outcomes ##################################### # use poisson loss for count/poisson outcomes subgrp.model.poisson < fit.subgroup(x = x, y = y.count, trt = trt01, propensity.func = prop.func, loss = "poisson_loss_lasso", type.measure = "mse", # option for cv.glmnet nfolds = 5) # option for cv.glmnet subgrp.model.poisson#> family: poisson #> loss: poisson_loss_lasso #> method: weighting #> cutpoint: 0 #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 31.2336 (n = 113) 17.191 (n = 88) #> Received 1 17.2722 (n = 140) 26.9694 (n = 159) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 13.9615 (n = 253) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 9.7785 (n = 247) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 4.09981 1.18199 0.01039 0.90901 4.15009 #> #>  #> #> Summary of individual treatment effects: #> E[YT=1, X]  E[YT=0, X] #> #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 60.31246 2.95433 0.02079 0.81282 2.07903 63.42390#################### Timetoevent outcomes ############################# library(survival) # \donttest{ 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 #> cutpoint: 0 #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 175.2226 (n = 140) 8.4595 (n = 61) #> Received 1 180.6631 (n = 203) 134.9956 (n = 96) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 5.4405 (n = 343) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 126.5361 (n = 157) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 0.87626 0.34859 0.14888 0.06438 0.87688 #> #>  #> #> Summary of individual treatment effects: #> E[YT=1, X] / E[YT=0, X] #> #> Note: for survival outcomes, the above ratio is #> E[g(Y)T=1, X] / E[g(Y)T=0, X], #> where g() is a monotone increasing function of Y, #> the survival time #> #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 0.4161 0.9376 1.1605 1.2129 1.4171 2.4019# } #################### Using custom loss functions ######################## ## Use custom loss function for binary outcomes fit.custom.loss.bin < function(x, y, weights, offset, ...) { df < data.frame(y = y, x) # minimize logistic loss with NO lasso penalty # with allowance for efficiency augmentation glmf < glm(y ~ x  1, weights = weights, offset = offset, # offset term allows for efficiency augmentation family = binomial(), ...) # save coefficients cfs = coef(glmf) # create prediction function. prd = function(x, type = "response") { dfte < cbind(1, x) colnames(dfte) < names(cfs) ## predictions must be returned on the scale ## of the linear predictor predict(glmf, data.frame(dfte), type = "link") } # return lost of required components list(predict = prd, model = glmf, coefficients = cfs) } # \donttest{ subgrp.model.bin.cust < fit.subgroup(x = x, y = y.binary, trt = trt01, propensity.func = prop.func, fit.custom.loss = fit.custom.loss.bin)#> Warning: noninteger #successes in a binomial glm!subgrp.model.bin.cust#> family: gaussian #> loss: custom #> method: weighting #> cutpoint: 0 #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 0.1761 (n = 109) 0.4712 (n = 92) #> Received 1 0.1818 (n = 149) 0.5965 (n = 150) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 0.0057 (n = 258) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 0.1254 (n = 242) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 3.32689 0.83519 0.03229 0.72308 3.06675 #> #>  #># } ## try exponential loss for ## positive outcomes fit.expo.loss < function(x, y, weights, ...) { expo.loss < function(beta, x, y, weights) { sum(weights * y * exp(drop(x %*% beta))) } # use optim() to minimize loss function opt < optim(rep(0, NCOL(x)), fn = expo.loss, x = x, y = y, weights = weights) coefs < opt$par pred < function(x, type = "response") { tcrossprod(cbind(1, x), t(coefs)) } # return list of required components list(predict = pred, model = opt, coefficients = coefs) } # \donttest{ # use exponential loss for positive outcomes subgrp.model.expo < fit.subgroup(x = x, y = y.count, trt = trt01, propensity.func = prop.func, fit.custom.loss = fit.expo.loss) subgrp.model.expo#> family: gaussian #> loss: custom #> method: weighting #> cutpoint: 0 #> propensity #> function: propensity.func #> #> benefit score: f(x), #> Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint' #> #> Average Outcomes: #> Recommended 0 Recommended 1 #> Received 0 30.2422 (n = 143) 15.1636 (n = 58) #> Received 1 14.8942 (n = 102) 27.2527 (n = 197) #> #> Treatment effects conditional on subgroups: #> Est of E[YT=0,Recom=0]E[YT=/=0,Recom=0] #> 15.3481 (n = 245) #> Est of E[YT=1,Recom=1]E[YT=/=1,Recom=1] #> 12.0891 (n = 255) #> #> NOTE: The above average outcomes are biased estimates of #> the expected outcomes conditional on subgroups. #> Use 'validate.subgroup()' to obtain unbiased estimates. #> #>  #> #> Benefit score quantiles (f(X) for 1 vs 0): #> 0% 25% 50% 75% 100% #> 1.41061 0.30465 0.01095 0.36604 1.82677 #> #>  #># }