fit.subgroup | R Documentation |
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_xgboost", "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 A-learning 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 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")) } |
fit.custom.loss |
A function which minimizes a user-specified
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"
.
predict |
A function that returns predictions of the covariate-conditional treatment effects |
model |
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 |
coefficients |
If the underlying subgroup identification model is parametric, |
call |
The call that produced the returned object. If |
family |
The family corresponding to the outcome provided |
loss |
The loss function used |
method |
The method used (either weighting or A-learning) |
propensity.func |
The propensity score function used |
larger.outcome.better |
If larger outcomes are preferred for this model |
cutpoint |
Benefit score cutoff value used for determining subgroups |
var.names |
The names of all variables used |
n.trts |
The number of treatment levels |
comparison.trts |
All treatment levels other than the reference level |
reference.trt |
The reference level for the treatment. This should usually be the control group/level |
trts |
All treatment levels |
trt.received |
The vector of treatment assignments |
pi.x |
A vector of propensity scores |
y |
A vector of outcomes |
benefit.scores |
A vector of conditional treatment effects, i.e. benefit scores |
recommended.trts |
A vector of treatment recommendations (i.e. for each patient, which treatment results in the best expected potential outcomes) |
subgroup.trt.effects |
(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 |
individual.trt.effects |
estimates of the individual treatment effects as returned by
|
Huling. J.D. and Yu, M. (2021), Subgroup Identification Using the personalized Package. Journal of Statistical Software 98(5), 1-60. doi:10.18637/jss.v098.i05
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 doi: 10.1111/biom.12676
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 doi: 10.1111/biom.12322
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
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 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 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))) # 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 } #################### Continuous outcomes ################################ subgrp.model <- fit.subgroup(x = x, y = y, trt = trt01, propensity.func = prop.func, loss = "sq_loss_lasso", # option for cv.glmnet, # better to use 'nfolds=10' nfolds = 3) summary(subgrp.model) # estimates of the individual-specific # treatment effect estimates: subgrp.model$individual.trt.effects # 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 #################### Using an augmentation function ##################### ## augmentation funcions involve modeling the conditional mean E[Y|T, X] ## and returning predictions that are averaged over the treatment values ## return <- 1/2 * (hat{E}[Y|T=1, X] + hat{E}[Y|T=-1, X]) ########################################################################## augment.func <- function(x, y, trt) { data <- data.frame(x, y, trt) xm <- model.matrix(y~trt*x-1, data = data) lmod <- cv.glmnet(y = y, x = xm) ## get predictions when trt = 1 data$trt <- 1 xm <- model.matrix(y~trt*x-1, data = data) preds_1 <- predict(lmod, xm, s = "lambda.min") ## get predictions when trt = -1 data$trt <- -1 xm <- model.matrix(y~trt*x-1, data = data) preds_n1 <- predict(lmod, xm, s = "lambda.min") ## return predictions averaged over trt return(0.5 * (preds_1 + preds_n1)) } subgrp.model.aug <- fit.subgroup(x = x, y = y, trt = trt01, propensity.func = prop.func, augment.func = augment.func, loss = "sq_loss_lasso", # option for cv.glmnet, # better to use 'nfolds=10' nfolds = 3) # option for cv.glmnet summary(subgrp.model.aug) #################### 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 = 3) # option for cv.glmnet subgrp.model.bin #################### 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 = 3) # option for cv.glmnet subgrp.model.poisson #################### Time-to-event outcomes ############################# 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 = 3) # option for cv.glmnet subgrp.model.cox #################### 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) } subgrp.model.bin.cust <- fit.subgroup(x = x, y = y.binary, trt = trt01, propensity.func = prop.func, fit.custom.loss = fit.custom.loss.bin) subgrp.model.bin.cust ## 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) } # 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.