tests/testthat/test-fitsubgroup.R
In jaredhuling/personalized: Estimation and Validation Methods for Subgroup Identification and Personalized Medicine
context("fit.subgroup")
test_that("test fit.subgroup for continuous outcomes and various losses", {
library(kernlab)
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
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
}
prop.func2 <- 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")
pi.x
}
prop.func3 <- 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]
dim(pi.x) <- NROW(pi.x)
pi.x
}
prop.func.bad <- function(x, trt, something)
{
return(numeric(NROW(trt)))
}
prop.func.bad2 <- function(x, trt, something, somethingelse)
{
return(numeric(NROW(trt)))
}
prop.func.bad3 <- function(x, something)
{
return(numeric(NROW(x)))
}
prop.func.bad4 <- function(x, trt)
{
ret <- numeric(NROW(x))
ret[5] <- 100
ret
}
if (Sys.info()[[1]] != "windows")
{
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad2,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad3,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad4,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = y, y = y,
trt = trt01,
propensity.func = prop.func))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
print(subgrp.model)
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
larger.outcome.better = FALSE,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
print(subgrp.model)
#### TEST CUSTOM LOSS FUNCTIONS
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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss2 <- function(x, y, weights, family, match.id,
n.trts, offset, trt, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad <- function(x, y, weights, wrongarg, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad2 <- function(x, y, weights, wrongarg, wrongarg2, wrongarg3, wrongarg4, wrongarg5,
just, toooo, many, arguments, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad3 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad4 <- function(x, y, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad5 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, badelement = 1000)
}
fit.custom.loss.bad6 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, badelement = 1000)
}
## predict element returned not a function
fit.custom.loss.bad7 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = cfs, model = lmf, coefficients = cfs)
}
## predict function returns something that's not the right length
fit.custom.loss.bad8 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x)[1:10,], t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bin <- function(x, y, weights, offset, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
glmf <- glm(y ~ x - 1, weights = weights,
offset = offset, # offset term allows for efficiency augmentation
family = binomial(), ...)
# save coefficients
cfs = unname(coef(glmf))
# create prediction function.
prd = function(x, type = "response")
{
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = glmf, coefficients = cfs)
}
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss,
loss = "custom")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
trt01_bad <- trt01
trt01_bad[4] <- "THE TRTMENT"
expect_error(subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01_bad,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss,
loss = "custom"))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "custom",
fit.custom.loss = fit.custom.loss)
# use loss with all possible args provided
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "custom",
fit.custom.loss = fit.custom.loss2)
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "custom"))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad2))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad4))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad5))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad6))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad7))
expect_warning(expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad8)))
}
###############################
if (Sys.info()[[1]] != "windows")
{
# test for factor trt
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
expect_warning(fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
reference.trt = 2,
loss = "owl_logistic_loss_lasso",
nfolds = 3))
## only 1 trt level
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(rep(1, NROW(y))),
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3))
## too many trt levels
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(1:NROW(y)),
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "owl_logistic_flip_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
}
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
trt01_bad <- trt01
trt01_bad[4] <- "third_trt_level"
expect_error(subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01_bad,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
nfolds = 3))
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_loss",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_loss",
nfolds = -5) )
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "vanilladot")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "polydot",
kpar = list(degree = 3),
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "laplacedot",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "besseldot",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
suppressWarnings({
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "anovadot",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
suppressWarnings({
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "splinedot",
margin = 0.2,
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
larger.outcome.better = FALSE,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
if (Sys.info()[[1]] != "windows")
{
# test if pi.x is a matrix with 1 column
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func2,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
# test if pi.x is a matrix with 1 column
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
# no prop func
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
}
if (Sys.info()[[1]] != "windows")
{
## tests for gam argument specification
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
method.gam = "REML",
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
method.gam = "REML",
loss = "sq_loss_lasso_gam")
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
control = gam.control(epsilon = 1e-10),
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
optimizer = "bfgs",
loss = "sq_loss_lasso_gam")
expect_is(subgrp.model, "subgroup_fitted")
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "sq_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
}
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "abs_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
})
test_that("test fit.subgroup for time-to-event outcomes and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
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
}
prop.func2 <- 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")
pi.x
}
subgrp.model <- 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
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
if (Sys.info()[[1]] != "windows")
{
# test for factor trt
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = trt01,
larger.outcome.better = FALSE,
propensity.func = prop.func,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
# test if pi.x is a matrix with 1 column
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = trt01,
propensity.func = prop.func2,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
}
})
test_that("test fit.subgroup with augment.func for continuous outcomes and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
xbeta <- xbeta + delta * trt
# continuous outcomes
y <- drop(xbeta) + rnorm(n.obs, sd = 2)
# count outcomes
y.count <- round(abs(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
}
augment.func <- function(x, y) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func2 <- function(x, y, trt) {lmod <- lm(y ~ x + trt); return(fitted(lmod))}
augment.func.bad <- function(x, y, something) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func.bad2 <- function(x, y, something, somethingelse) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func.bad3 <- function(x, something) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func.bad4 <- function(x, y) {lmod <- lm(y ~ x); return(fitted(lmod)[1:10])}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
if (Sys.info()[[1]] != "windows")
{
## test for method which uses adjustment directly
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func,
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func2,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad2,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad3,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad4,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "sq_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "sq_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "abs_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
expect_warning(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_xgboost",
nrounds = 5,
params = list(max_depth = 1, eta = 0.1,
booster = "gbtree"),
nfold = 3))
# expect_warning(fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "abs_loss_gbm",
# n.trees = 5,
# cv.folds = 1,
# n.cores = 1))
# subgrp.model <- fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
# expect_warning(fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gbm",
# n.trees = 5,
# cv.folds = 1,
# n.cores = 1))
}
if (Sys.info()[[1]] != "windows")
{
# subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
# trt = trt01,
# propensity.func = prop.func,
# loss = "cox_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
# expect_warning(fit.subgroup(x = x, y = Surv(y.time.to.event, status),
# trt = trt01,
# propensity.func = prop.func,
# loss = "cox_loss_gbm",
# n.trees = 5,
# cv.folds = 1,
# n.cores = 1))
}
})
test_that("test fit.subgroup for continuous outcomes with no propensity function specified", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
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)
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
}
})
test_that("test fit.subgroup for binary outcomes and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1] + x[,5]
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(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.binary,
trt = trt01,
propensity.func = prop.func,
loss = "logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y.count,
trt = trt01,
propensity.func = prop.func,
loss = "poisson_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y.count,
trt = trt01,
propensity.func = prop.func,
loss = "poisson_loss_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y.count,
trt = trt01,
propensity.func = prop.func,
loss = "poisson_loss_lasso_gam")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
# subgrp.model <- fit.subgroup(x = x, y = y.count,
# trt = trt01,
# propensity.func = prop.func,
# loss = "poisson_loss_gbm")
#
# expect_is(subgrp.model, "subgroup_fitted")
#
# invisible(capture.output(print(subgrp.model, digits = 2)))
#
# invisible(capture.output(summary(subgrp.model)))
augment.func <- function(x, y) {
lmod <- glm(y ~ x, family = binomial());
return(predict(lmod, type = "link"))
}
subgrp.modela <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.modela, "subgroup_fitted")
invisible(capture.output(print(subgrp.modela, digits = 2)))
invisible(capture.output(summary(subgrp.modela)))
suppressWarnings({
subgrp.modela <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "logistic_loss_gam",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.modela, "subgroup_fitted")
invisible(capture.output(print(subgrp.modela, digits = 2)))
invisible(capture.output(summary(subgrp.modela)))
suppressWarnings({
subgrp.modelg <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
loss = "logistic_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.modelg, "subgroup_fitted")
invisible(capture.output(print(subgrp.modelg, digits = 2)))
invisible(capture.output(summary(subgrp.modelg)))
suppressWarnings({
subgrp.modelga <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "logistic_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.modelga, "subgroup_fitted")
invisible(capture.output(print(subgrp.modelga, digits = 2)))
invisible(capture.output(summary(subgrp.modelga)))
# subgrp.model <- fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gam")
#
# expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
}
# subgrp.model <- fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gbm", n.trees = 5,
# n.cores = 1)
#
# expect_is(subgrp.model, "subgroup_fitted")
#
# invisible(capture.output(print(subgrp.model, digits = 2)))
#
# invisible(capture.output(summary(subgrp.model)))
})
test_that("test fit.subgroup for continuous outcomes and match.id provided", {
set.seed(123)
n.obs <- 1000
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- -3 + 0.5 * x[,1] - 0.5 * x[,5]
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 + 5 * sin(x[,2] ^ 2) - x[,3] )
xbeta <- x[,1]
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",
type.measure = "auc")
pi.x <- predict(propens.model, s = "lambda.min",
newx = x, type = "response")[,1]
pi.x
}
n.matches <- 2
pscore <- prop.func(x, trt)
match.mat <- array(NA, dim = c(sum(trt01), n.matches) )
trt.idx <- which(trt01 == 1)
ctrl.idx <- which(trt01 == 0)
for (m in 1:n.matches)
{
for (i in 1:length(trt.idx))
{
best.match <- which.min(abs(pscore[trt.idx[i]] - pscore[ctrl.idx]))
# remove match from pool
ctrl.idx <- ctrl.idx[ctrl.idx != ctrl.idx[best.match]]
match.mat[i,m] <- ctrl.idx[best.match]
}
}
n.total <- nrow(match.mat) * (ncol(match.mat) + 1)
match.id.mat <- matrix(nrow = nrow(match.mat), ncol = ncol(match.mat))
for (cc in 1:ncol(match.mat))
{
match.id.mat[,cc] <- 1:nrow(match.mat)
}
match.id <- c(1:nrow(match.mat), as.vector(match.id.mat))
match.idx <- c(trt.idx, as.vector(match.mat))
x.m <- x[match.idx,]
y.m <- y[match.idx]
trt.m <- trt01[match.idx]
y.time.to.event.m <- y.time.to.event[match.idx]
status.m <- status[match.idx]
y.binary.m <- y.binary[match.idx]
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
invisible(capture.output(print(subgrp.model.m, digits = 2)))
invisible(capture.output(summary(subgrp.model.m)))
if (Sys.info()[[1]] != "windows")
{
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
suppressWarnings({
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "owl_hinge_loss",
kernel = "vanilladot",
C = c(0.1, 0.5, 1, 2, 10)[4:5],
eps = 1e-4,
margin = 5)
})
expect_is(subgrp.model.m, "subgroup_fitted")
suppressWarnings({
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "owl_hinge_flip_loss",
kernel = "vanilladot",
C = c(0.1, 0.5, 1, 2, 10)[4:5],
eps = 1e-4,
margin = 5)
})
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
penlty.factor = rep(1, ncol(x.m) + 1),
loss = "sq_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
expect_error(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "sq_loss_lasso_gam",
nfolds = 2))
expect_warning(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "sq_loss_lasso_gam",
foldid = sample(1:5, nrow(x.m), replace = TRUE)))
# subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
# trt = trt.m,
# match.id = match.id,
# loss = "sq_loss_gbm", n.trees = 5, n.cores = 1)
#
# expect_is(subgrp.model.m, "subgroup_fitted")
# subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
# trt = trt.m,
# match.id = match.id,
# loss = "abs_loss_gbm", n.trees = 5, n.cores = 1)
#
# expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso_gam")
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
propensity.func = prop.func,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.val <- validate.subgroup(subgrp.model.m, B = 2,
method = "train")
expect_is(subgrp.val, "subgroup_validated")
invisible(capture.output(print(subgrp.val, digits = 2)))
subgrp.val <- validate.subgroup(subgrp.model.m, B = 2,
method = "boot")
expect_is(subgrp.val, "subgroup_validated")
invisible(capture.output(print(subgrp.val, digits = 2)))
expect_warning(subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
foldid = sample(1:5, nrow(x.m), replace = TRUE)))
expect_is(subgrp.model.m, "subgroup_fitted")
expect_warning(subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
match.id = as.factor(match.id),
loss = "cox_loss_lasso",
foldid = sample(1:5, nrow(x.m), replace = TRUE)))
expect_is(subgrp.model.m, "subgroup_fitted")
expect_error(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
nfolds = 2) )
expect_error(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
propensity.func = prop.func,
match.id = rep(1, length(y.m)), # provide bad match.id (only 1 level)
loss = "sq_loss_lasso",
nfolds = 3) )
expect_error(fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
match.id = as.factor(match.id),
loss = "cox_loss_lasso",
nfolds = 2) )
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
penalty.factor = rep(1, ncol(x.m) + 1),
loss = "sq_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
match.id = as.factor(match.id),
loss = "cox_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
# subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
# trt = trt.m,
# match.id = as.factor(match.id),
# loss = "cox_loss_gbm", n.trees = 5, n.cores = 1)
#
# expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
augment.func = function(x,y,trt) {return(rep(1, NROW(x)))},
match.id = as.factor(match.id),
penalty.factor = rep(1, ncol(x.m) + 1),
loss = "cox_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
## parallel
subgrp.val <- validate.subgroup(subgrp.model.m, B = 10,
parallel = TRUE,
method = "training")
subgrp.val <- validate.subgroup(subgrp.model.m, B = 10,
parallel = TRUE,
method = "boot")
}
})
test_that("test fit.subgroup for continuous outcomes and multiple trts and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulated non-randomized treatment with multiple levels
xbetat_1 <- 0.15 + 0.5 * x[,1] - 0.25 * x[,5]
xbetat_2 <- 0.15 - 0.5 * x[,2] + 0.25 * x[,3]
trt.1.prob <- exp(xbetat_1) / (1 + exp(xbetat_1) + exp(xbetat_2))
trt.2.prob <- exp(xbetat_2) / (1 + exp(xbetat_1) + exp(xbetat_2))
trt.3.prob <- 1 - (trt.1.prob + trt.2.prob)
prob.mat <- cbind(trt.1.prob, trt.2.prob, trt.3.prob)
trt <- apply(prob.mat, 1, function(rr) rmultinom(1, 1, prob = rr))
trt <- apply(trt, 2, function(rr) which(rr == 1))
# simulate response
delta1 <- 8 * (0.5 + x[,2] - x[,3] )
delta2 <- 4 * (-0.5 + x[,1] - x[,5] + x[,4] )
xbeta <- x[,1] - 2 * x[,2] - 3 * x[,5]
xbeta <- xbeta + delta1 * ((trt == 1) - (trt == 3) ) + delta2 * ((trt == 2) - (trt == 3) )
# 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)
# use multinomial logistic regression model with lasso penalty for propensity
propensity.multinom.lasso <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- propens[,match(levels(trt), colnames(propens))]
probs
}
# use multinomial logistic regression model with lasso penalty for propensity
propensity.multinom.lasso.bad <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- propens[,match(levels(trt), colnames(propens))][,1:2]
probs
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "owl_logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
expect_warning(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "owl_logistic_loss_lasso",
method = "a_learning",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "owl_logistic_flip_loss_lasso",
nfolds = 3))
## test matching without prop score
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
match.id = rep(1:10, each = 10),
loss = "sq_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
summ <- summarize.subgroups(subgrp.model)
expect_is(summ, "data.frame")
invisible(capture.output(print(summ, digits = 2, p.value = 0.25)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso.bad,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) )
# provide incorrect reference trt
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
reference.trt = "something_wrong",
loss = "sq_loss_lasso",
nfolds = 3) )
# no prop func
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
larger.outcome.better = FALSE,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
}
# use multinomial logistic regression model with lasso penalty for propensity
propensity.multinom.lasso <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- propens[cbind(1:nrow(propens), match(levels(trt)[trt], colnames(propens)))]
probs
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
summ <- summarize.subgroups(subgrp.model)
expect_is(summ, "data.frame")
invisible(capture.output(print(summ, digits = 2, p.value = 0.25)))
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
larger.outcome.better = FALSE,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = prop.func,
loss = "sq_loss_lasso_gam",
nfolds = 3))
# different outcomes
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y.binary,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) )
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "cox_loss_lasso",
nfolds = 3) )
propensity.multinom.lasso.array <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- array(dim = c(dim(propens), 2, 4))
probs[is.na(probs)] <- 0.5
probs
}
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso.array,
loss = "cox_loss_lasso",
nfolds = 3) )
}
})
jaredhuling/personalized documentation built on Sept. 10, 2022, 11:35 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
context("fit.subgroup")
test_that("test fit.subgroup for continuous outcomes and various losses", {
library(kernlab)
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
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
}
prop.func2 <- 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")
pi.x
}
prop.func3 <- 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]
dim(pi.x) <- NROW(pi.x)
pi.x
}
prop.func.bad <- function(x, trt, something)
{
return(numeric(NROW(trt)))
}
prop.func.bad2 <- function(x, trt, something, somethingelse)
{
return(numeric(NROW(trt)))
}
prop.func.bad3 <- function(x, something)
{
return(numeric(NROW(x)))
}
prop.func.bad4 <- function(x, trt)
{
ret <- numeric(NROW(x))
ret[5] <- 100
ret
}
if (Sys.info()[[1]] != "windows")
{
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad2,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad3,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func.bad4,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = y, y = y,
trt = trt01,
propensity.func = prop.func))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
print(subgrp.model)
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
larger.outcome.better = FALSE,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
print(subgrp.model)
#### TEST CUSTOM LOSS FUNCTIONS
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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss2 <- function(x, y, weights, family, match.id,
n.trts, offset, trt, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad <- function(x, y, weights, wrongarg, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad2 <- function(x, y, weights, wrongarg, wrongarg2, wrongarg3, wrongarg4, wrongarg5,
just, toooo, many, arguments, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad3 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad4 <- function(x, y, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
lmf <- lm(y ~ x - 1, weights = weights, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bad5 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, badelement = 1000)
}
fit.custom.loss.bad6 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, badelement = 1000)
}
## predict element returned not a function
fit.custom.loss.bad7 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = cfs, model = lmf, coefficients = cfs)
}
## predict function returns something that's not the right length
fit.custom.loss.bad8 <- 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, ...)
# save coefficients
cfs = unname(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")
{
# roxygen2 removes the below, so
# using tcrossprod instead
#cbind(1, x)
tcrossprod(cbind(1, x)[1:10,], t(cfs))
}
# return lost of required components
list(predict = prd, model = lmf, coefficients = cfs)
}
fit.custom.loss.bin <- function(x, y, weights, offset, ...) {
df <- data.frame(y = y, x)
# minimize squared error loss with NO lasso penalty
glmf <- glm(y ~ x - 1, weights = weights,
offset = offset, # offset term allows for efficiency augmentation
family = binomial(), ...)
# save coefficients
cfs = unname(coef(glmf))
# create prediction function.
prd = function(x, type = "response")
{
#cbind(1, x)
tcrossprod(cbind(1, x), t(cfs))
}
# return lost of required components
list(predict = prd, model = glmf, coefficients = cfs)
}
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss,
loss = "custom")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
trt01_bad <- trt01
trt01_bad[4] <- "THE TRTMENT"
expect_error(subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01_bad,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss,
loss = "custom"))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "custom",
fit.custom.loss = fit.custom.loss)
# use loss with all possible args provided
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "custom",
fit.custom.loss = fit.custom.loss2)
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "custom"))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad2))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad4))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad5))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad6))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad7))
expect_warning(expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
fit.custom.loss = fit.custom.loss.bad8)))
}
###############################
if (Sys.info()[[1]] != "windows")
{
# test for factor trt
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
expect_warning(fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
reference.trt = 2,
loss = "owl_logistic_loss_lasso",
nfolds = 3))
## only 1 trt level
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(rep(1, NROW(y))),
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3))
## too many trt levels
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(1:NROW(y)),
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "owl_logistic_flip_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
}
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
trt01_bad <- trt01
trt01_bad[4] <- "third_trt_level"
expect_error(subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01_bad,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
nfolds = 3))
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_loss",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_loss",
nfolds = -5) )
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "vanilladot")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "polydot",
kpar = list(degree = 3),
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "laplacedot",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "besseldot",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
suppressWarnings({
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "anovadot",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
suppressWarnings({
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "owl_hinge_flip_loss",
kernel = "splinedot",
margin = 0.2,
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
larger.outcome.better = FALSE,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
if (Sys.info()[[1]] != "windows")
{
# test if pi.x is a matrix with 1 column
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func2,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
# test if pi.x is a matrix with 1 column
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
# no prop func
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
}
if (Sys.info()[[1]] != "windows")
{
## tests for gam argument specification
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
method.gam = "REML",
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
method.gam = "REML",
loss = "sq_loss_lasso_gam")
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
control = gam.control(epsilon = 1e-10),
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
optimizer = "bfgs",
loss = "sq_loss_lasso_gam")
expect_is(subgrp.model, "subgroup_fitted")
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "sq_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
}
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "abs_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
})
test_that("test fit.subgroup for time-to-event outcomes and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
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
}
prop.func2 <- 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")
pi.x
}
subgrp.model <- 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
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
if (Sys.info()[[1]] != "windows")
{
# test for factor trt
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = as.factor(trt01),
propensity.func = prop.func,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = trt01,
larger.outcome.better = FALSE,
propensity.func = prop.func,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
# test if pi.x is a matrix with 1 column
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = trt01,
propensity.func = prop.func2,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
}
})
test_that("test fit.subgroup with augment.func for continuous outcomes and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
xbeta <- xbeta + delta * trt
# continuous outcomes
y <- drop(xbeta) + rnorm(n.obs, sd = 2)
# count outcomes
y.count <- round(abs(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
}
augment.func <- function(x, y) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func2 <- function(x, y, trt) {lmod <- lm(y ~ x + trt); return(fitted(lmod))}
augment.func.bad <- function(x, y, something) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func.bad2 <- function(x, y, something, somethingelse) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func.bad3 <- function(x, something) {lmod <- lm(y ~ x); return(fitted(lmod))}
augment.func.bad4 <- function(x, y) {lmod <- lm(y ~ x); return(fitted(lmod)[1:10])}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
if (Sys.info()[[1]] != "windows")
{
## test for method which uses adjustment directly
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func,
propensity.func = prop.func,
loss = "owl_logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func2,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad2,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad3,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt01,
augment.func = augment.func.bad4,
propensity.func = prop.func,
loss = "sq_loss_lasso",
nfolds = 3))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "sq_loss_gam")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "sq_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model)))
invisible(capture.output(summary(subgrp.model)))
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "sq_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
# subgrp.model <- fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "abs_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
expect_warning(fit.subgroup(x = x, y = y,
trt = trt01,
propensity.func = prop.func,
loss = "sq_loss_xgboost",
nrounds = 5,
params = list(max_depth = 1, eta = 0.1,
booster = "gbtree"),
nfold = 3))
# expect_warning(fit.subgroup(x = x, y = y,
# trt = trt01,
# propensity.func = prop.func,
# loss = "abs_loss_gbm",
# n.trees = 5,
# cv.folds = 1,
# n.cores = 1))
# subgrp.model <- fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
# expect_warning(fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gbm",
# n.trees = 5,
# cv.folds = 1,
# n.cores = 1))
}
if (Sys.info()[[1]] != "windows")
{
# subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
# trt = trt01,
# propensity.func = prop.func,
# loss = "cox_loss_gbm",
# n.trees = 5,
# n.cores = 1)
#
# invisible(capture.output(print(subgrp.model)))
# invisible(capture.output(summary(subgrp.model)))
# expect_is(subgrp.model, "subgroup_fitted")
# expect_warning(fit.subgroup(x = x, y = Surv(y.time.to.event, status),
# trt = trt01,
# propensity.func = prop.func,
# loss = "cox_loss_gbm",
# n.trees = 5,
# cv.folds = 1,
# n.cores = 1))
}
})
test_that("test fit.subgroup for continuous outcomes with no propensity function specified", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1]
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)
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt01,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
}
})
test_that("test fit.subgroup for binary outcomes and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- 0.5 + 0.5 * x[,1] - 0.5 * x[,5]
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] )
xbeta <- x[,1] + x[,5]
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(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.binary,
trt = trt01,
propensity.func = prop.func,
loss = "logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y.count,
trt = trt01,
propensity.func = prop.func,
loss = "poisson_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y.count,
trt = trt01,
propensity.func = prop.func,
loss = "poisson_loss_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
subgrp.model <- fit.subgroup(x = x, y = y.count,
trt = trt01,
propensity.func = prop.func,
loss = "poisson_loss_lasso_gam")
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
# subgrp.model <- fit.subgroup(x = x, y = y.count,
# trt = trt01,
# propensity.func = prop.func,
# loss = "poisson_loss_gbm")
#
# expect_is(subgrp.model, "subgroup_fitted")
#
# invisible(capture.output(print(subgrp.model, digits = 2)))
#
# invisible(capture.output(summary(subgrp.model)))
augment.func <- function(x, y) {
lmod <- glm(y ~ x, family = binomial());
return(predict(lmod, type = "link"))
}
subgrp.modela <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.modela, "subgroup_fitted")
invisible(capture.output(print(subgrp.modela, digits = 2)))
invisible(capture.output(summary(subgrp.modela)))
suppressWarnings({
subgrp.modela <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "logistic_loss_gam",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.modela, "subgroup_fitted")
invisible(capture.output(print(subgrp.modela, digits = 2)))
invisible(capture.output(summary(subgrp.modela)))
suppressWarnings({
subgrp.modelg <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
loss = "logistic_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.modelg, "subgroup_fitted")
invisible(capture.output(print(subgrp.modelg, digits = 2)))
invisible(capture.output(summary(subgrp.modelg)))
suppressWarnings({
subgrp.modelga <- fit.subgroup(x = x, y = y.binary,
trt = trt01,
propensity.func = prop.func,
augment.func = augment.func,
loss = "logistic_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
})
expect_is(subgrp.modelga, "subgroup_fitted")
invisible(capture.output(print(subgrp.modelga, digits = 2)))
invisible(capture.output(summary(subgrp.modelga)))
# subgrp.model <- fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gam")
#
# expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
}
# subgrp.model <- fit.subgroup(x = x, y = y.binary,
# trt = trt01,
# propensity.func = prop.func,
# loss = "logistic_loss_gbm", n.trees = 5,
# n.cores = 1)
#
# expect_is(subgrp.model, "subgroup_fitted")
#
# invisible(capture.output(print(subgrp.model, digits = 2)))
#
# invisible(capture.output(summary(subgrp.model)))
})
test_that("test fit.subgroup for continuous outcomes and match.id provided", {
set.seed(123)
n.obs <- 1000
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulate non-randomized treatment
xbetat <- -3 + 0.5 * x[,1] - 0.5 * x[,5]
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 + 5 * sin(x[,2] ^ 2) - x[,3] )
xbeta <- x[,1]
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",
type.measure = "auc")
pi.x <- predict(propens.model, s = "lambda.min",
newx = x, type = "response")[,1]
pi.x
}
n.matches <- 2
pscore <- prop.func(x, trt)
match.mat <- array(NA, dim = c(sum(trt01), n.matches) )
trt.idx <- which(trt01 == 1)
ctrl.idx <- which(trt01 == 0)
for (m in 1:n.matches)
{
for (i in 1:length(trt.idx))
{
best.match <- which.min(abs(pscore[trt.idx[i]] - pscore[ctrl.idx]))
# remove match from pool
ctrl.idx <- ctrl.idx[ctrl.idx != ctrl.idx[best.match]]
match.mat[i,m] <- ctrl.idx[best.match]
}
}
n.total <- nrow(match.mat) * (ncol(match.mat) + 1)
match.id.mat <- matrix(nrow = nrow(match.mat), ncol = ncol(match.mat))
for (cc in 1:ncol(match.mat))
{
match.id.mat[,cc] <- 1:nrow(match.mat)
}
match.id <- c(1:nrow(match.mat), as.vector(match.id.mat))
match.idx <- c(trt.idx, as.vector(match.mat))
x.m <- x[match.idx,]
y.m <- y[match.idx]
trt.m <- trt01[match.idx]
y.time.to.event.m <- y.time.to.event[match.idx]
status.m <- status[match.idx]
y.binary.m <- y.binary[match.idx]
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
invisible(capture.output(print(subgrp.model.m, digits = 2)))
invisible(capture.output(summary(subgrp.model.m)))
if (Sys.info()[[1]] != "windows")
{
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
suppressWarnings({
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "owl_hinge_loss",
kernel = "vanilladot",
C = c(0.1, 0.5, 1, 2, 10)[4:5],
eps = 1e-4,
margin = 5)
})
expect_is(subgrp.model.m, "subgroup_fitted")
suppressWarnings({
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "owl_hinge_flip_loss",
kernel = "vanilladot",
C = c(0.1, 0.5, 1, 2, 10)[4:5],
eps = 1e-4,
margin = 5)
})
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
penlty.factor = rep(1, ncol(x.m) + 1),
loss = "sq_loss_lasso_gam",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
expect_error(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "sq_loss_lasso_gam",
nfolds = 2))
expect_warning(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = match.id,
loss = "sq_loss_lasso_gam",
foldid = sample(1:5, nrow(x.m), replace = TRUE)))
# subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
# trt = trt.m,
# match.id = match.id,
# loss = "sq_loss_gbm", n.trees = 5, n.cores = 1)
#
# expect_is(subgrp.model.m, "subgroup_fitted")
# subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
# trt = trt.m,
# match.id = match.id,
# loss = "abs_loss_gbm", n.trees = 5, n.cores = 1)
#
# expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso_gam")
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
propensity.func = prop.func,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.val <- validate.subgroup(subgrp.model.m, B = 2,
method = "train")
expect_is(subgrp.val, "subgroup_validated")
invisible(capture.output(print(subgrp.val, digits = 2)))
subgrp.val <- validate.subgroup(subgrp.model.m, B = 2,
method = "boot")
expect_is(subgrp.val, "subgroup_validated")
invisible(capture.output(print(subgrp.val, digits = 2)))
expect_warning(subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
foldid = sample(1:5, nrow(x.m), replace = TRUE)))
expect_is(subgrp.model.m, "subgroup_fitted")
expect_warning(subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
match.id = as.factor(match.id),
loss = "cox_loss_lasso",
foldid = sample(1:5, nrow(x.m), replace = TRUE)))
expect_is(subgrp.model.m, "subgroup_fitted")
expect_error(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso",
nfolds = 2) )
expect_error(fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
propensity.func = prop.func,
match.id = rep(1, length(y.m)), # provide bad match.id (only 1 level)
loss = "sq_loss_lasso",
nfolds = 3) )
expect_error(fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
match.id = as.factor(match.id),
loss = "cox_loss_lasso",
nfolds = 2) )
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
loss = "sq_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = y.m,
trt = trt.m,
match.id = as.factor(match.id),
penalty.factor = rep(1, ncol(x.m) + 1),
loss = "sq_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
match.id = as.factor(match.id),
loss = "cox_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
# subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
# trt = trt.m,
# match.id = as.factor(match.id),
# loss = "cox_loss_gbm", n.trees = 5, n.cores = 1)
#
# expect_is(subgrp.model.m, "subgroup_fitted")
subgrp.model.m <- fit.subgroup(x = x.m, y = Surv(y.time.to.event.m, status.m),
trt = trt.m,
augment.func = function(x,y,trt) {return(rep(1, NROW(x)))},
match.id = as.factor(match.id),
penalty.factor = rep(1, ncol(x.m) + 1),
loss = "cox_loss_lasso")
expect_is(subgrp.model.m, "subgroup_fitted")
## parallel
subgrp.val <- validate.subgroup(subgrp.model.m, B = 10,
parallel = TRUE,
method = "training")
subgrp.val <- validate.subgroup(subgrp.model.m, B = 10,
parallel = TRUE,
method = "boot")
}
})
test_that("test fit.subgroup for continuous outcomes and multiple trts and various losses", {
set.seed(123)
n.obs <- 100
n.vars <- 5
x <- matrix(rnorm(n.obs * n.vars, sd = 3), n.obs, n.vars)
# simulated non-randomized treatment with multiple levels
xbetat_1 <- 0.15 + 0.5 * x[,1] - 0.25 * x[,5]
xbetat_2 <- 0.15 - 0.5 * x[,2] + 0.25 * x[,3]
trt.1.prob <- exp(xbetat_1) / (1 + exp(xbetat_1) + exp(xbetat_2))
trt.2.prob <- exp(xbetat_2) / (1 + exp(xbetat_1) + exp(xbetat_2))
trt.3.prob <- 1 - (trt.1.prob + trt.2.prob)
prob.mat <- cbind(trt.1.prob, trt.2.prob, trt.3.prob)
trt <- apply(prob.mat, 1, function(rr) rmultinom(1, 1, prob = rr))
trt <- apply(trt, 2, function(rr) which(rr == 1))
# simulate response
delta1 <- 8 * (0.5 + x[,2] - x[,3] )
delta2 <- 4 * (-0.5 + x[,1] - x[,5] + x[,4] )
xbeta <- x[,1] - 2 * x[,2] - 3 * x[,5]
xbeta <- xbeta + delta1 * ((trt == 1) - (trt == 3) ) + delta2 * ((trt == 2) - (trt == 3) )
# 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)
# use multinomial logistic regression model with lasso penalty for propensity
propensity.multinom.lasso <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- propens[,match(levels(trt), colnames(propens))]
probs
}
# use multinomial logistic regression model with lasso penalty for propensity
propensity.multinom.lasso.bad <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- propens[,match(levels(trt), colnames(propens))][,1:2]
probs
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "owl_logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
expect_warning(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "owl_logistic_loss_lasso",
method = "a_learning",
nfolds = 3))
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "owl_logistic_flip_loss_lasso",
nfolds = 3))
## test matching without prop score
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
match.id = rep(1:10, each = 10),
loss = "sq_loss_lasso",
nfolds = 3)
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
summ <- summarize.subgroups(subgrp.model)
expect_is(summ, "data.frame")
invisible(capture.output(print(summ, digits = 2, p.value = 0.25)))
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso.bad,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) )
# provide incorrect reference trt
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
reference.trt = "something_wrong",
loss = "sq_loss_lasso",
nfolds = 3) )
# no prop func
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
larger.outcome.better = FALSE,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
}
# use multinomial logistic regression model with lasso penalty for propensity
propensity.multinom.lasso <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- propens[cbind(1:nrow(propens), match(levels(trt)[trt], colnames(propens)))]
probs
}
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
invisible(capture.output(print(subgrp.model, digits = 2)))
invisible(capture.output(summary(subgrp.model)))
summ <- summarize.subgroups(subgrp.model)
expect_is(summ, "data.frame")
invisible(capture.output(print(summ, digits = 2, p.value = 0.25)))
if (Sys.info()[[1]] != "windows")
{
subgrp.model <- fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = propensity.multinom.lasso,
reference.trt = 3,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y,
trt = as.factor(trt),
reference.trt = "2",
larger.outcome.better = FALSE,
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = y,
trt = trt,
propensity.func = prop.func,
loss = "sq_loss_lasso_gam",
nfolds = 3))
# different outcomes
subgrp.model <- fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "cox_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
subgrp.model <- fit.subgroup(x = x, y = y.binary,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "logistic_loss_lasso",
nfolds = 3) # option for cv.glmnet
expect_is(subgrp.model, "subgroup_fitted")
expect_error(fit.subgroup(x = x, y = Surv(y.time.to.event, status),
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "sq_loss_lasso",
nfolds = 3) )
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso,
loss = "cox_loss_lasso",
nfolds = 3) )
propensity.multinom.lasso.array <- function(x, trt)
{
if (!is.factor(trt)) trt <- as.factor(trt)
gfit <- cv.glmnet(y = trt, x = x, family = "multinomial")
# predict returns a matrix of probabilities:
# one column for each treatment level
propens <- drop(predict(gfit, newx = x, type = "response", s = "lambda.min",
nfolds = 3, alpha = 0))
# return the probability corresponding to the
# treatment that was observed
probs <- array(dim = c(dim(propens), 2, 4))
probs[is.na(probs)] <- 0.5
probs
}
expect_error(fit.subgroup(x = x, y = y,
trt = as.factor(trt),
propensity.func = propensity.multinom.lasso.array,
loss = "cox_loss_lasso",
nfolds = 3) )
}
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.