Nothing
#' Calculation of covariate-conditional treatment effects
#'
#' @description Calculates covariate conditional treatment effects using estimated benefit scores
#' @param ... not used
#' @return A List with elements \code{delta} (if the treatment effects are a difference/contrast,
#' i.e. \eqn{E[Y|T=1, X] - E[Y|T=-1, X]}) and \code{gamma} (if the treatment effects are a ratio,
#' i.e. \eqn{E[Y|T=1, X] / E[Y|T=-1, X]})
#' @export
treatment.effects <- function(x, ...) UseMethod("treatment.effects")
all_losses <- 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")
calc_treatment_effects <- function(benefit.scores,
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"),
pi.x = NULL)
{
loss <- match.arg(loss)
method <- match.arg(method)
benefit.scores <- drop(benefit.scores)
if (!is.null(pi.x))
{
pi.x <- drop(pi.x)
}
if (method == "a_learning" & is.null(pi.x) &
!(grepl("sq_loss", loss) | grepl("cox_loss", loss)) )
{
stop("propensity score vector 'pi.x' must be supplied for method = 'a_learning'")
}
trt_eff_delta <- trt_eff_gamma <- NA
if (grepl("sq_loss", loss))
{
if (method == "weighting")
{
trt_eff_delta <- 2 * benefit.scores
} else
{
trt_eff_delta <- benefit.scores
}
} else if (grepl("owl_logistic_loss", loss) )
{
if (method == "weighting")
{
trt_eff_gamma <- exp(benefit.scores)
} else
{
trt_eff_gamma <- (1 + exp((1 - pi.x) * benefit.scores)) / (1 + exp(-pi.x * benefit.scores))
}
} else if (grepl("cox_loss", loss) )
{
trt_eff_gamma <- exp(-benefit.scores)
} else if (grepl("poisson_loss", loss) )
{
if (method == "weighting")
{
trt_eff_delta <- exp(benefit.scores) - exp(-benefit.scores)
} else
{
trt_eff_delta <- exp((1 - pi.x) * benefit.scores) - exp(-pi.x * benefit.scores)
}
} else if (grepl("logistic_loss", loss) & !grepl("owl_", loss))
{
if (method == "weighting")
{
exp_bs <- exp(benefit.scores)
trt_eff_delta <- (exp_bs - 1) / (exp_bs + 1)
} else
{
exp_bs <- exp(benefit.scores)
exp_pi_bs <- exp(pi.x * benefit.scores)
trt_eff_delta <- ((exp_bs - 1) / (exp_pi_bs + 1)) * (exp_pi_bs / (exp_pi_bs + exp_bs))
}
} else
{
warning(paste("treatment effects not available for loss:", loss) )
}
effects <- list(delta = trt_eff_delta,
gamma = trt_eff_gamma)
class(effects) <- c("individual_treatment_effects", class(effects) )
effects
}
#' @rdname treatment.effects
#' @export
treatment.effects.default <- function(x, ...)
{
treatment.effects.subgroup_fitted(x, ...)
}
#' @param benefit.scores vector of estimated benefit scores
#' @param loss loss choice USED TO CALCULATE \code{benefit.scores} of both the M function from Chen, et al (2017) and
#' potentially the penalty used for variable selection. See \code{\link[personalized]{fit.subgroup}} for more details.
#' @param method method choice USED TO CALCULATE \code{benefit.scores}. Either the \code{"weighting"} method or
#' \code{"a_learning"} method. See \code{\link[personalized]{fit.subgroup}} for more details
#' @param pi.x The propensity score for each observation
#' @rdname treatment.effects
#' @export
treat.effects <- function(benefit.scores,
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"),
pi.x = NULL, ...)
{
loss <- match.arg(loss)
method <- match.arg(method)
calc_treatment_effects(benefit.scores, loss, method, pi.x)
}
#' @seealso \code{\link[personalized]{fit.subgroup}} for function which fits subgroup identification models.
#' @param x a fitted object from \code{fit.subgroup()}
#' @seealso \code{\link[personalized]{print.individual_treatment_effects}} for printing of objects returned by
#' \code{treat.effects} or \code{treatment.effects}
#' @rdname treatment.effects
#' @export
#' @examples
#' library(personalized)
#'
#' set.seed(123)
#' n.obs <- 500
#' n.vars <- 25
#' x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
#'
#'
#' # simulate non-randomized treatment
#' xbetat <- 0.5 + 0.5 * x[,21] - 0.5 * x[,11]
#' 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]
#' xbeta <- xbeta + delta * trt
#'
#' # continuous outcomes
#' y <- drop(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
#' }
#'
#' subgrp.model <- fit.subgroup(x = x, y = y,
#' trt = trt01,
#' propensity.func = prop.func,
#' loss = "sq_loss_lasso",
#' nfolds = 3) # option for cv.glmnet
#'
#' trt_eff <- treatment.effects(subgrp.model)
#' str(trt_eff)
#'
#' trt_eff
#'
#'
#' 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
#'
#' trt_eff_c <- treatment.effects(subgrp.model.cox)
#' str(trt_eff_c)
#'
#' trt_eff_c
#'
treatment.effects.subgroup_fitted <- function(x, ...)
{
treat.effects(x$benefit.scores,
x$loss,
x$method,
x$pi.x)
}
#' Printing individualized treatment effects
#'
#' @description Prints results for estimated subgroup treatment effects
#'
#' @param x a fitted object from either \code{\link[personalized]{treat.effects}} or \code{\link[personalized]{treatment.effects}}
#' @param digits minimal number of significant digits to print.
#' @param ... further arguments passed to or from \code{\link[base]{print.default}}.
#' @export
print.individual_treatment_effects <- function(x, digits = max(getOption('digits')-3, 3), ...)
{
if (!is.na(x$delta[1]))
{
comp <- attr(x$delta, "comparison.trts")
ref <- attr(x$delta, "reference.trt")
trts <- attr(x$delta, "trts")
if (NCOL(x$delta) == 1)
{
cat(paste0("Summary of individual treatment effects: \nE[Y|T=", comp, ", X] - E[Y|T=", ref, ", X]\n\n"))
print(summary(x$delta), digits = digits, ...)
} else if (NCOL(x$delta) == length(trts)) ##multinomial model case
{
cat(paste0("Summary of individual treatment effects: \nE[Y|T=", "trt", ", X]\nwhere 'trt' is ", paste(trts, collapse = " and "), "\n\n"))
print(summary(x$delta), digits = digits, ...)
} else
{
cat(paste0("Summary of individual treatment effects: \nE[Y|T=", "trt", ", X] - E[Y|T=", ref, ", X]\nwhere 'trt' is ", paste(comp, collapse = " and "), "\n\n"))
print(summary(x$delta), digits = digits, ...)
}
}
if (!is.na(x$gamma[1]))
{
comp <- attr(x$gamma, "comparison.trts")
ref <- attr(x$gamma, "reference.trt")
trts <- attr(x$gamma, "trts")
if (NCOL(x$gamma) == 1)
{
cat(paste0("Summary of individual treatment effects: \nE[Y|T=", comp, ", X] / E[Y|T=", ref, ", X]\n\n"))
cat(paste0("Note: for survival outcomes, the above ratio is \nE[g(Y)|T=", comp, ", X] / E[g(Y)|T=", ref, ", X], \nwhere g() is a monotone increasing function of Y, \nthe survival time\n\n"))
print(summary(x$gamma), digits = digits, ...)
} else if (NCOL(x$gamma) == length(trts)) ##multinomial model case
{
cat(paste0("Summary of individual treatment effects: \nE[Y|T=", "trt", ", X]\nwhere 'trt' is ", paste(trts, collapse = " and "), "\n\n"))
print(summary(x$gamma), digits = digits, ...)
} else
{
cat(paste0("Summary of individual treatment effects: \nE[Y|T=", "trt", ", X] / E[Y|T=", ref, ", X]\nwhere 'trt' is ", paste(comp, collapse = " and "), "\n\n"))
cat("Note: for survival outcomes, the above ratio is \nE[g(Y)|T=1, X] / E[g(Y)|T=-1, X], \nwhere g() is a monotone increasing function of Y, \nthe survival time\n\n")
print(summary(x$gamma), digits = digits, ...)
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.