Nothing
### test-predictCox.R ---
#----------------------------------------------------------------------
## author: Brice Ozenne
## created: sep 4 2017 (10:38)
## Version:
## last-updated: mar 7 2023 (18:31)
## By: Brice Ozenne
## Update #: 177
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
## * Settings
library(riskRegression)
library(testthat)
library(rms)
library(survival)
library(data.table)
library(timereg); n.sim <- 500;
context("function predictCox")
## * [predictCox] Baseline hazard (no strata)
cat("[predictCox] Estimation of the baseline hazard (no strata) \n")
## ** Data
data(Melanoma, package = "riskRegression")
## ** Model
fit.coxph <- coxph(Surv(time,status == 1) ~ thick + invasion + ici, data = Melanoma, y = TRUE, x = TRUE)
fit.cph <- cph(Surv(time,status == 1) ~ thick + invasion + ici, data = Melanoma, y = TRUE, x = TRUE)
## ** Compare to survival::basehaz
test_that("baseline hazard (no strata): compare to survival::basehaz",{
## vs basehaz
expect_equal(ignore_attr=TRUE,predictCox(fit.coxph, centered = FALSE)$cumhazard,
survival::basehaz(fit.coxph, centered = FALSE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predictCox(fit.coxph, centered = TRUE)$cumhazard,
survival::basehaz(fit.coxph, centered = TRUE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predictCox(fit.cph)$cumhazard,
survival::basehaz(fit.cph)$hazard, tolerance = 1e-8)
## consistency cph coxph
## possible differences due to different fit - coef(fit.coxph)-coef(fit.cph)
expect_equal(ignore_attr=TRUE,predictCox(fit.cph, centered = FALSE),
predictCox(fit.coxph, centered = FALSE),
tolerance = 100*max(abs(coef(fit.coxph)-coef(fit.cph))))
## note centered = TRUE will no give same results as coxph do not center all variables
})
## ** Number of events
test_that("baseline hazard (no strata): number of events",{
## find all unique event times
GS.alltime <- sort(unique(Melanoma$time))
RR.coxph <- predictCox(fit.coxph)
expect_equal(ignore_attr=TRUE,RR.coxph$times, GS.alltime)
RR.cph <- predictCox(fit.cph)
expect_equal(ignore_attr=TRUE,RR.cph$times, GS.alltime)
})
## * [predictCox] Baseline hazard (strata)
cat("[predictCox] Estimation of the baseline hazard (strata) \n")
## ** Data
data(Melanoma, package = "riskRegression")
## ** Model
fitS.coxph <- coxph(Surv(time,status == 1) ~ thick + strata(invasion) + strata(ici), data = Melanoma, y = TRUE, x = TRUE)
fitS.cph <- cph(Surv(time,status == 1) ~ thick + strat(invasion) + strat(ici), data = Melanoma, y = TRUE, x = TRUE)
## ** Compare to survival::basehaz
test_that("baseline hazard (strata): compare to survival::basehaz",{
## vs basehaz
expect_equal(ignore_attr=TRUE,predictCox(fitS.coxph, centered = FALSE)$cumhazard,
basehaz(fitS.coxph, centered = FALSE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predictCox(fitS.coxph, centered = TRUE)$cumhazard,
basehaz(fitS.coxph, centered = TRUE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predictCox(fitS.cph)$cumhazard,
basehaz(fitS.cph)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predictCox(fitS.coxph, keep.strata = TRUE)$strata,
basehaz(fitS.coxph)$strata)
expect_equal(ignore_attr=TRUE,predictCox(fitS.cph, keep.strata = TRUE)$strata,
basehaz(fitS.cph)$strata)
## consistency cph coxph
## different ordering of the strata
e.coxph <- as.data.table(predictCox(fitS.coxph))
e.cph <- as.data.table(predictCox(fitS.cph))
levels(e.coxph$strata)
levels(e.cph$strata)
})
## ** Number of events
test_that("baseline hazard (strata): number of events",{
GS.alltime <- tapply(Melanoma$time, interaction(Melanoma$ici,Melanoma$invasion),
function(x){sort(unique(x))})
RR.coxph <- predictCox(fitS.coxph)
test.alltime <- tapply(RR.coxph$times, RR.coxph$strata, function(x){x})
expect_equal(ignore_attr=TRUE,unname(GS.alltime),unname(test.alltime))
GS.alltime <- tapply(Melanoma$time, interaction(Melanoma$invasion,Melanoma$ici),
function(x){sort(unique(x))})
RR.cph <- predictCox(fitS.cph)
test.alltime <- tapply(RR.cph$times, RR.cph$strata, function(x){x})
expect_equal(ignore_attr=TRUE,unname(GS.alltime),unname(test.alltime))
})
## * [predictCox] Baseline hazard with time varying covariates (no strata)
cat("[predictCox] Estimation of the baseline hazard (time varying cov, no strata) \n")
## ** Data
## example from help(coxph)
dt.TV <- list(start=c(1,2,5,2,1,7,3,4,8,8),
stop=c(2,3,6,7,8,9,9,9,14,17),
event=c(1,1,1,1,1,1,1,0,0,0),
x=c(1,0,0,1,0,1,1,1,0,0))
## ** Model
fit.coxphTV <- coxph(Surv(start, stop, event) ~ x, data = dt.TV, x = TRUE, y = TRUE)
fit.cphTV <- cph(Surv(start, stop, event) ~ x, data = dt.TV, x = TRUE, y = TRUE)
## ** Compare to survival::basehaz
test_that("baseline hazard (no strata, time varying): compare to survival::basehaz",{
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fit.coxphTV, centered = FALSE)$cumhazard),
basehaz(fit.coxphTV, centered = FALSE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fit.coxphTV, centered = TRUE)$cumhazard),
basehaz(fit.coxphTV, centered = TRUE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fit.cphTV)$cumhazard),
basehaz(fit.cphTV)$hazard, tolerance = 1e-8)
})
## * [predictCox] Baseline hazard with time varying covariates (strata)
cat("[predictCox] Estimation of the baseline hazard (time varying cov, strata) \n")
## ** Data
set.seed(10)
dtS.TV <- rbind(cbind(as.data.table(dt.TV),S = 1),
cbind(as.data.table(dt.TV),S = 2))
dtS.TV[, randomS := rbinom(.N,size = 1, prob = 1/2)]
## ** Model
fitS1.coxphTV <- coxph(Surv(start, stop, event) ~ strata(S) + x, data = dtS.TV, x = TRUE, y = TRUE)
fitS1.cphTV <- cph(Surv(start, stop, event) ~ strat(S) + x, data = dtS.TV, x = TRUE, y = TRUE)
fitS2.coxphTV <- coxph(Surv(start, stop, event) ~ strata(randomS) + x, data = dtS.TV, x = TRUE, y = TRUE)
fitS2.cphTV <- cph(Surv(start, stop, event) ~ strat(randomS) + x, data = dtS.TV, x = TRUE, y = TRUE)
## ** Compare to survival::basehaz
test_that("baseline hazard (strata, time varying): compare to survival::basehaz",{
## strata defined by S
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fitS1.coxphTV, centered = FALSE)$cumhazard),
basehaz(fitS1.coxphTV, centered = FALSE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fitS1.coxphTV, centered = TRUE)$cumhazard),
basehaz(fitS1.coxphTV, centered = TRUE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fitS1.cphTV)$cumhazard),
basehaz(fitS1.cphTV)$hazard, tolerance = 1e-8)
## strata defined by randomS
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fitS2.coxphTV, centered = FALSE)$cumhazard),
basehaz(fitS2.coxphTV, centered = FALSE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fitS2.coxphTV, centered = TRUE)$cumhazard),
basehaz(fitS2.coxphTV, centered = TRUE)$hazard, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,suppressWarnings(predictCox(fitS2.cphTV)$cumhazard),
basehaz(fitS2.cphTV)$hazard, tolerance = 1e-8)
})
## * [predictCox] Predictions,se,band,average.iid (no strata, continuous variables)
cat("[predictCox] Predictions (no strata, continuous) \n")
## ** Data
set.seed(10)
dt <- sampleData(5e1, outcome = "survival")[,.(time,event,X1,X2,X6)]
dt[,X1:=as.numeric(as.character(X1))]
dt[,X2:=as.numeric(as.character(X2))]
dt[ , X16 := X1*X6]
## sorted dataset
dt.sort <- copy(dt)
setkeyv(dt.sort,c("time"))
## ** Model
e.coxph <- coxph(Surv(time, event) ~ X1*X6, data = dt, y = TRUE, x = TRUE)
e.cph <- cph(Surv(time, event) ~ X1*X6, data = dt, y = TRUE, x = TRUE)
e.coxph_sort <- coxph(Surv(time, event) ~ X1*X6, data = dt.sort, y = TRUE, x = TRUE)
e.timereg <- cox.aalen(Surv(time, event) ~ prop(X1) + prop(X6) + prop(X1*X6),
data = dt, resample.iid = TRUE, max.timepoint.sim=NULL)
## ** Consistency between hazard/cumhazard/survival
test_that("[predictCox] - consistency of hazard/cumhazard/survival",{
predRR <- predictCox(e.coxph, type = c("hazard","cumhazard","survival"), times = sort(dt$time), newdata = dt)
expect_equal(ignore_attr=TRUE,predRR$hazard[,-1], t(apply(predRR$cumhazard,1,diff)), tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predRR$survival, exp(-predRR$cumhazard), tolerance = 1e-8)
})
## ** One time
## *** Extract information
predGS <- predict(e.timereg, newdata = dt, times = 10)
predRR1 <- predictCox(e.coxph, newdata = dt, times = 10,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR2 <- predictCox(e.coxph_sort, newdata = dt, times = 10,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR1.none <- confint(predRR1, survival.transform = "none", seed = 10,
n.sim = n.sim)
predRR1.loglog <- confint(predRR1, survival.transform = "loglog", seed = 10,
n.sim = n.sim)
## *** Test vs. cph
test_that("[predictCox] compare survival and survival.se coxph/cph (1 fixed time)",{
res.cph <- predictCox(e.cph, newdata = dt, times = 10,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1),
as.data.table(res.cph),
tolerance = 1e-3)
})
## *** Test vs. timereg
test_that("[predictCox] compare survival and survival.se to timereg (1 fixed time)",{
## punctual estimate
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predGS$S0))
expect_equal(ignore_attr=TRUE,as.double(predRR2$survival), as.double(predGS$S0))
## standard error
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predGS$se.S0))
expect_equal(ignore_attr=TRUE,as.double(predRR2$survival.se), as.double(predGS$se.S0))
})
## *** Test vs. known result
test_that("[confint.predictCox] compare to known values (1 fixed time, no transformation)",{
## confidence intervals/band
GS <- data.table("observation" = c(6, 7),
"times" = c(10, 10),
"survival" = c(0.19073, 0.016148),
"survival.se" = c(0.088812, 0.023259),
"survival.lower" = c(0.016662, 0),
"survival.upper" = c(0.364797, 0.061734),
"survival.quantileBand" = c(2.021964, 2.03706),
"survival.lowerBand" = c(0.011156, 0),
"survival.upperBand" = c(0.370304, 0.063527))
## butils::object2script(as.data.table(predRR1.none)[6:7,names(GS),with=FALSE], digit = 6)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1.none)[6:7,names(GS),with=FALSE], GS, tolerance = 1e-4, scale = 1)
})
test_that("[confint.predictCox] compare to known values (1 fixed time, log log transformation)",{
GS <- data.table("observation" = c(6, 7),
"times" = c(10, 10),
"survival" = c(0.19073, 0.01615),
"survival.se" = c(0.088812, 0.023259),
"survival.lower" = c(0.05646, 0.00028),
"survival.upper" = c(0.38475, 0.12474),
"survival.quantileBand" = c(2.02196, 2.03706),
"survival.lowerBand" = c(0.05368, 0.00022),
"survival.upperBand" = c(0.39115, 0.13183))
## butils::object2script(as.data.table(predRR1.loglog)[6:7,names(GS),with=FALSE], digit = 5)
expect_equal(ignore_attr=TRUE,
as.data.table(predRR1.loglog)[6:7,names(GS),with=FALSE],
GS,
tolerance = 1e-4,
scale = 1)
})
## ** At event times
## *** Extract information
vec.time <- sort(dt$time[1:10])
set.seed(10)
predGS <- predict(e.timereg, newdata = dt, times = vec.time)
predRR1 <- predictCox(e.coxph, newdata = dt, times = vec.time,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR2 <- predictCox(e.coxph_sort, newdata = dt, times = vec.time,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR1.none <- confint(predRR1, survival.transform = "none", seed = 10, n.sim = n.sim)
predRR1.loglog <- confint(predRR1, survival.transform = "loglog", seed = 10, n.sim = n.sim)
## *** Test vs. cph
test_that("[predictCox] compare survival and survival.se coxph/cph (eventtimes)",{
res.cph <- predictCox(e.cph, newdata = dt, times = vec.time,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1),
as.data.table(res.cph),
tolerance = 1e-3)
})
## *** Test vs. timereg
test_that("[predictCox] compare survival and survival.se to timereg (eventtimes)",{
## punctual estimate
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predGS$S0))
expect_equal(ignore_attr=TRUE,as.double(predRR2$survival), as.double(predGS$S0))
## standard error
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predGS$se.S0))
expect_equal(ignore_attr=TRUE,as.double(predRR2$survival.se), as.double(predGS$se.S0))
})
## *** Test vs. known result
test_that("[confint.predictCox] compare to known values (eventtimes, no transformation)",{
## confidence intervals/band
GS <- data.table("observation" = c(6, 7),
"times" = c(0.170031, 0.170031),
"survival" = c(0.993101, 0.982909),
"survival.se" = c(0.007437, 0.018876),
"survival.lower" = c(0.978526, 0.945913),
"survival.upper" = c(1, 1),
"survival.quantileBand" = c(2.666239, 2.66367),
"survival.lowerBand" = c(0.973273, 0.93263),
"survival.upperBand" = c(1, 1))
## butils::object2script(as.data.table(predRR1.none)[6:7,names(GS),with=FALSE], digit = 6)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1.none)[6:7,names(GS),with=FALSE], GS, tolerance = 1e-4, scale = 1)
})
test_that("[confint.predictCox] compare to known values (eventtimes, log log transformation)",{
GS <- data.table("observation" = c(6, 7),
"times" = c(0.17003, 0.17003),
"survival" = c(0.9931, 0.98291),
"survival.se" = c(0.007437, 0.018876),
"survival.lower" = c(0.94395, 0.85811),
"survival.upper" = c(0.99917, 0.99806),
"survival.quantileBand" = c(2.66624, 2.66367),
"survival.lowerBand" = c(0.88353, 0.71525),
"survival.upperBand" = c(0.99961, 0.99911))
## butils::object2script(as.data.table(predRR1.loglog)[6:7,names(GS),with=FALSE], digit = 5)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1.loglog)[6:7,names(GS),with=FALSE],
GS, tolerance = 1e-4, scale = 1)
})
## ** After last event
test_that("[predictCox] after the last event",{
predRR1 <- predictCox(e.coxph, newdata = dt, times = 1e8,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR1 <- confint(predRR1, n.sim = n.sim)
expect_true(all(is.na(predRR1$survival)))
expect_true(all(is.na(predRR1$survival.se)))
expect_true(all(is.na(predRR1$survival.lower)))
expect_true(all(is.na(predRR1$survival.upper)))
expect_true(all(is.na(predRR1$survival.lowerBand)))
expect_true(all(is.na(predRR1$survival.upperBand)))
vec.time <- max(dt$time) + c(-1e-1,0,1e-1)
predRR <- predictCox(e.coxph, type = c("hazard","cumhazard","survival"),
times = vec.time, newdata = dt)
M.true <- matrix(c(FALSE, FALSE, TRUE), nrow = NROW(dt), ncol = 3, byrow = TRUE)
expect_equal(ignore_attr=TRUE,is.na(predRR$hazard), M.true)
expect_equal(ignore_attr=TRUE,is.na(predRR$cumhazard), M.true)
expect_equal(ignore_attr=TRUE,is.na(predRR$survival), M.true)
})
test_that("Prediction - last event censored",{
dt.lastC <- copy(dt)
Utimes <- sort(unique(dt$time))
n.Utimes <- length(Utimes)
dt.lastC[time==max(time), event := 0]
fit <- coxph(Surv(time, event == 1) ~ X1 + X2 + X6, data = dt.lastC, y = TRUE, x = TRUE)
predictRR <- predictCox(fit, newdata = dt.lastC, times = tail(Utimes, 2))
survLast <- predictRR$survival[,2]
survLastM1 <- predictRR$survival[,1]
expect_true(all(survLast-survLastM1==0))
})
test_that("Prediction - last event death",{
dt.lastD <- copy(dt)
Utimes <- sort(unique(dt$time))
n.Utimes <- length(Utimes)
dt.lastD[time==max(time), event := 1]
fit <- coxph(Surv(time, event == 1) ~ X1 + X2 + X6, data = dt.lastD, y = TRUE, x = TRUE)
predictRR <- predictCox(fit, newdata = dt.lastD[1], times = tail(Utimes, 2))
survLast <- predictRR$survival[,2]
survLastM1 <- predictRR$survival[,1]
expect_true(all(survLast<survLastM1))
})
## ** Before first event
test_that("[predictCox] before the first event",{
predRR1 <- predictCox(e.coxph, newdata = dt, times = 1e-8,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR1 <- confint(predRR1, n.sim = n.sim)
expect_true(all(predRR1$survival==1))
expect_true(all(predRR1$cumhazard==0))
expect_true(all(predRR1$survival.se==0))
expect_true(all(predRR1$survival.lower==1))
expect_true(all(predRR1$survival.upper==1))
expect_true(all(predRR1$survival.lowerBand==1))
expect_true(all(predRR1$survival.upperBand==1))
})
## ** Sorted vs. unsorted times
test_that("[predictCox] - sorted vs. unsorted times (no strata)",{
vec.time <- dt$time
index.sort <- order(vec.time)
vec.time_sort <- dt$time[index.sort]
predRR1 <- predictCox(e.coxph, newdata = dt, times = vec.time, se = TRUE)
predRR2 <- predictCox(e.coxph, newdata = dt, times = vec.time_sort, se = TRUE)
expect_equal(ignore_attr=TRUE,predRR1$survival[,index.sort],
predRR2$survival)
expect_equal(ignore_attr=TRUE,predRR1$survival.se[,index.sort],
predRR2$survival.se)
expect_equal(ignore_attr=TRUE,predRR1$time[index.sort],
predRR2$time)
})
## ** iid.average
test_that("[predictCox] - fast iid average (no strata)",{
## simple average
predRR.av <- predictCox(e.coxph, times = dt$time[1:5], average.iid = TRUE, newdata = dt,
type = c("hazard","cumhazard","survival"))
predRR.GS <- predictCox(e.coxph, times = dt$time[1:5], iid = TRUE, newdata = dt,
type = c("hazard","cumhazard","survival"))
expect_equal(ignore_attr=TRUE,apply(predRR.GS$hazard.iid, MARGIN = 1:2,mean),
predRR.av$hazard.average.iid, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,apply(predRR.GS$cumhazard.iid, MARGIN = 1:2,mean),
predRR.av$cumhazard.average.iid, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,apply(predRR.GS$survival.iid, MARGIN = 1:2,mean),
predRR.av$survival.average.iid, tolerance = 1e-8)
## weighted average
fT <- TRUE
attr(fT, "factor") <- list(matrix(5, nrow = NROW(dt), ncol = 5),
matrix(1:NROW(dt), nrow = NROW(dt), ncol = 5)
)
predRR.av2 <- predictCox(e.coxph, times = sort(dt$time[1:5]), average.iid = fT, newdata = dt,
type = c("hazard","cumhazard","survival"))
calcGS <- function(iid){
apply(iid, MARGIN = 1:2, function(iCol){
mean(iCol * attr(fT, "factor")[[2]][,1])
})
}
expect_equal(ignore_attr=TRUE,predRR.av$hazard.average.iid[,order(dt$time[1:5])]*5,
predRR.av2$hazard.average.iid[[1]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predRR.av$cumhazard.average.iid[,order(dt$time[1:5])]*5,
predRR.av2$cumhazard.average.iid[[1]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predRR.av$survival.average.iid[,order(dt$time[1:5])]*5,
predRR.av2$survival.average.iid[[1]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,calcGS(predRR.GS$hazard.iid)[,order(dt$time[1:5])],
predRR.av2$hazard.average.iid[[2]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,calcGS(predRR.GS$cumhazard.iid)[,order(dt$time[1:5])],
predRR.av2$cumhazard.average.iid[[2]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,calcGS(predRR.GS$survival.iid)[,order(dt$time[1:5])],
predRR.av2$survival.average.iid[[2]],
tolerance = 1e-8)
})
## * [predictCox] Predictions,se,band,average.iid (no strata, categorical variable)
cat("[predictCox] Predictions (no strata, categorical) \n")
## ** Data
set.seed(10)
dt <- sampleData(5e1, outcome = "survival")[,.(time,event,X1,X2,X6)]
dt[,X1:=as.numeric(as.character(X1))]
dt[,X2:=as.numeric(as.character(X2))]
dt[ , Xcat2 := as.factor(paste0(X1,X2))]
## sorted dataset
dt.sort <- copy(dt)
setkeyv(dt.sort,c("time"))
## ** Model
e.coxph <- coxph(Surv(time, event) ~ Xcat2 + X6 , data = dt, y = TRUE, x = TRUE)
e.cph <- cph(Surv(time, event) ~ Xcat2 + X6 , data = dt, y = TRUE, x = TRUE)
e.timereg <- cox.aalen(Surv(time, event) ~ prop(Xcat2) + prop(X6),
data = dt, resample.iid = TRUE, max.timepoint.sim=NULL)
## ** One time and event times
test_that("[predictCox] compare survival and survival.se to timereg/cph (categorical variable)",{
## fixed time
predGS <- predict(e.timereg, newdata = dt, times = 10)
predRR1 <- predictCox(e.coxph, newdata = dt, times = 10, se = TRUE)
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predGS$S0))
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predGS$se.S0))
predCPH <- predictCox(e.cph, newdata = dt, times = 10, se = TRUE)
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predCPH$survival), tolerance = 1e-3)
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predCPH$survival.se), tolerance = 1e-3)
## event time
predGS <- predict(e.timereg, newdata = dt, times = sort(dt$time))
predRR1 <- predictCox(e.coxph, newdata = dt, times = sort(dt$time), se = TRUE)
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predGS$S0))
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predGS$se.S0))
predCPH <- predictCox(e.cph, newdata = dt, times = sort(dt$time), se = TRUE)
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predCPH$survival), tolerance = 1e-3)
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predCPH$survival.se), tolerance = 1e-3)
})
## * [predictCox] Predictions,se,band,average.iid (strata)
cat("[predictCox] Predictions (strata) \n")
## ** Data
set.seed(10)
dtStrata <- copy(dt)
dtStrata[, strata := rbinom(n = .N, size = 2, prob = c(1/3,1/2))] # strata
dtStrata.sort <- copy(dtStrata)
setkeyv(dtStrata.sort, c("strata", "time"))
## ** Model
eS.timereg <- cox.aalen(Surv(time, event) ~ strata(strata)-1 + prop(X1) + prop(X6),
data = dtStrata,
resample.iid = TRUE, max.timepoint.sim=NULL)
eS.cph <- cph(Surv(time, event) ~ strat(strata) + X1 + X6,
data = dtStrata, y = TRUE, x = TRUE)
eS.coxph <- coxph(Surv(time, event) ~ strata(strata) + X1 + X6,
data = dtStrata, y = TRUE, x = TRUE)
## ** Reject incorrect strata
test_that("[predictCox] - incorrect strata",{
dt2 <- copy(dt)
dt2$strata <- "5616"
expect_error(suppressWarnings(predictCox(eS.coxph, times = dt2$time, newdata = dt2)))
})
## ** 1 fixed time
## *** Extract information
set.seed(10)
predGS <- predict(eS.timereg, newdata = dtStrata, times = 4)
predRR1 <- predictCox(eS.coxph, newdata = dtStrata, times = 4,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR1.none <- confint(predRR1, survival.transform = "none", seed = 10,
n.sim = n.sim)
predRR1.loglog <- confint(predRR1, survival.transform = "loglog", seed = 10,
n.sim = n.sim)
## *** Test vs. cph
test_that("[predictCox] compare survival and survival.se coxph/cph (1 fixed time, strata)",{
set.seed(10)
res.cph <- predictCox(eS.cph, newdata = dtStrata, times = 4,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1),
as.data.table(res.cph),
tolerance = 1e-3)
})
## *** Test vs. timereg
test_that("[predictCox] compare survival and survival.se to timereg (1 fixed time, strata)",{
## punctual estimate
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predGS$S0))
## standard error
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predGS$se.S0))
})
## *** Test vs. known result
test_that("[confint.predictCox] compare to known values (1 fixed time, no transformation, strata)",{
## confidence intervals/band
GS <- data.table("observation" = c(6, 7),
"times" = c(4, 4),
"survival" = c(0.780713, 0.519771),
"survival.se" = c(0.082657, 0.12765),
"survival.lower" = c(0.618708, 0.269583),
"survival.upper" = c(0.942717, 0.76996),
"survival.quantileBand" = c(1.855252, 1.944986),
"survival.lowerBand" = c(0.627364, 0.271494),
"survival.upperBand" = c(0.934062, 0.768048))
## butils::object2script(as.data.table(predRR1.none)[6:7,names(GS),with=FALSE], digit = 6)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1.none)[6:7,names(GS),with=FALSE], GS, tolerance = 1e-4, scale = 1)
})
test_that("[confint.predictCox] compare to known values (1 fixed time, log log transformation, strata)",{
GS <- data.table("observation" = c(6, 7),
"times" = c(4, 4),
"survival" = c(0.78071, 0.51977),
"survival.se" = c(0.082657, 0.12765),
"survival.lower" = c(0.56416, 0.25526),
"survival.upper" = c(0.89848, 0.73082),
"survival.quantileBand" = c(1.85525, 1.94499),
"survival.lowerBand" = c(0.57849, 0.25722),
"survival.upperBand" = c(0.89408, 0.72953))
## butils::object2script(as.data.table(predRR1.loglog)[6:7,names(GS),with=FALSE], digit = 5)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1.loglog)[6:7,names(GS),with=FALSE],
GS, tolerance = 1e-4, scale = 1)
})
## ** At event times (in one strata)
## *** Extract information
vec.time <- sort(dtStrata$time)[1:10]
set.seed(10)
predGS <- predict(eS.timereg, newdata = dtStrata, times = vec.time)
predRR1 <- predictCox(eS.coxph, newdata = dtStrata, times = vec.time,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR1.none <- confint(predRR1, survival.transform = "none", seed = 10,
n.sim = n.sim)
predRR1.loglog <- confint(predRR1, survival.transform = "loglog", seed = 10,
n.sim = n.sim)
## *** Test vs. cph
test_that("[predictCox] compare survival and survival.se coxph/cph (eventtime, strata)",{
set.seed(10)
res.cph <- predictCox(eS.cph, newdata = dtStrata, times = vec.time,
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1),
as.data.table(res.cph),
tolerance = 1e-3)
})
## *** Test vs. timereg
test_that("[predictCox] compare survival and survival.se to timereg (eventtimes, strata)",{
## punctual estimate
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival), as.double(predGS$S0))
## standard error
expect_equal(ignore_attr=TRUE,as.double(predRR1$survival.se), as.double(predGS$se.S0))
})
## *** Test vs. known result
test_that("[confint.predictCox] compare to known values (eventtimes, no transformation, strata)",{
## confidence intervals/band
GS <- data.table("observation" = c(35, 36),
"times" = c(0.877076, 0.877076),
"survival" = c(0.994832, 0.992565),
"survival.se" = c(0.00369, 0.010196),
"survival.lower" = c(0.987599, 0.97258),
"survival.upper" = c(1, 1),
"survival.quantileBand" = c(2.276138, 2.331861),
"survival.lowerBand" = c(0.986432, 0.968788),
"survival.upperBand" = c(1, 1))
## butils::object2script(as.data.table(predRR1.none)[135:136,names(GS),with=FALSE], digit = 6)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1.none)[135:136,names(GS),with=FALSE], GS, tolerance = 1e-4, scale = 1)
})
test_that("[confint.predictCox] compare to known values (eventtimes, log log transformation, strata)",{
GS <- data.table("observation" = c(35, 36),
"times" = c(0.87708, 0.87708),
"survival" = c(0.99483, 0.99256),
"survival.se" = c(0.00369, 0.010196),
"survival.lower" = c(0.97914, 0.89511),
"survival.upper" = c(0.99873, 0.9995),
"survival.quantileBand" = c(2.27614, 2.33186),
"survival.lowerBand" = c(0.97391, 0.83119),
"survival.upperBand" = c(0.99898, 0.9997))
## butils::object2script(as.data.table(predRR1.loglog)[135:136,names(GS),with=FALSE], digit = 5)
expect_equal(ignore_attr=TRUE,as.data.table(predRR1.loglog)[135:136,names(GS),with=FALSE],
GS, tolerance = 1e-4, scale = 1)
})
## ** after last event
test_that("[predictCox] after the last event (strata)",{
lastevent <- dtStrata[, max(time), by = "strata"]
laststrata <- lastevent[V1==max(V1),strata]
predRR1 <- predictCox(eS.coxph, newdata = dtStrata, times = max(lastevent[["V1"]]),
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
predRR1 <- confint(predRR1, n.sim = n.sim)
expect_true(all(is.na(predRR1$survival[dtStrata$strata!=laststrata,])))
expect_true(all(is.na(predRR1$cumhazard[dtStrata$strata!=laststrata,])))
expect_true(all(!is.na(predRR1$survival[dtStrata$strata==laststrata,])))
expect_true(all(!is.na(predRR1$cumhazard[dtStrata$strata==laststrata,])))
expect_true(all(is.na(predRR1$survival.se[dtStrata$strata!=laststrata,])))
expect_true(all(is.na(predRR1$survival.lower[dtStrata$strata!=laststrata,])))
expect_true(all(is.na(predRR1$survival.upper[dtStrata$strata!=laststrata,])))
expect_true(all(is.na(predRR1$survival.lowerBand[dtStrata$strata!=laststrata,])))
expect_true(all(is.na(predRR1$survival.upperBand[dtStrata$strata!=laststrata,])))
})
## ** before first event
test_that("[predictCox] before the first event (strata)",{
firstevent <- dtStrata[, min(time), by = "strata"]
firststrata <- firstevent[V1==min(V1),strata]
predRR1 <- predictCox(eS.coxph, newdata = dtStrata, times = min(firstevent[["V1"]]))
expect_true(all(predRR1$survival[dtStrata$strata!=firststrata,]==1))
expect_true(all(predRR1$cumhazard[dtStrata$strata!=firststrata,]==0))
expect_true(all(predRR1$survival[dtStrata$strata==firststrata,]<1))
expect_true(all(predRR1$cumhazard[dtStrata$strata==firststrata,]>0))
})
## ** iid.average
test_that("[predictCox] - iid average",{
## eS.coxph <- coxph(Surv(time, event) ~ strata(strata), data = dtStrata, x = TRUE)
## seqTime <- c(0,sort(dtStrata$time)[1:5])
seqTime <- c(0,dtStrata$time[1:5])
## simple average
predRR.av <- predictCox(eS.coxph, times = seqTime, average.iid = TRUE, newdata = dtStrata,
type = c("hazard","cumhazard","survival"))
predRR.GS <- predictCox(eS.coxph, times = seqTime, iid = TRUE, newdata = dtStrata,
type = c("hazard","cumhazard","survival"))
expect_equal(ignore_attr=TRUE,apply(predRR.GS$hazard.iid, MARGIN = 1:2,mean),
predRR.av$hazard.average.iid, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,apply(predRR.GS$cumhazard.iid, MARGIN = 1:2,mean),
predRR.av$cumhazard.average.iid, tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,apply(predRR.GS$survival.iid, MARGIN = 1:2,mean),
predRR.av$survival.average.iid, tolerance = 1e-8)
## weighted average
fT <- TRUE
attr(fT, "factor") <- list(matrix(5, nrow = NROW(dtStrata), ncol = 6),
matrix(1:NROW(dtStrata), nrow = NROW(dtStrata), ncol = 6)
)
predRR.av2 <- predictCox(eS.coxph, times = sort(seqTime), average.iid = fT, newdata = dtStrata,
type = c("hazard","cumhazard","survival"))
calcGS <- function(iid){
apply(iid, MARGIN = 1:2, function(iCol){
mean(iCol * attr(fT, "factor")[[2]][,1])
})
}
expect_equal(ignore_attr=TRUE,predRR.av$hazard.average.iid[,order(seqTime)]*5,
predRR.av2$hazard.average.iid[[1]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predRR.av$cumhazard.average.iid[,order(seqTime)]*5,
predRR.av2$cumhazard.average.iid[[1]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predRR.av$survival.average.iid[,order(seqTime)]*5,
predRR.av2$survival.average.iid[[1]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,calcGS(predRR.GS$hazard.iid)[,order(seqTime)],
predRR.av2$hazard.average.iid[[2]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,calcGS(predRR.GS$cumhazard.iid)[,order(seqTime)],
predRR.av2$cumhazard.average.iid[[2]],
tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,calcGS(predRR.GS$survival.iid)[,order(seqTime)],
predRR.av2$survival.average.iid[[2]],
tolerance = 1e-8)
})
## * [predictCox] SE/CI check against manual computation
cat("[predictCox] SE/CI check against manual computation \n")
## from confint.predictCox
## ** Data
set.seed(10)
dt <- sampleData(40,outcome="survival")
## ** Model
fit <- coxph(Surv(time,event)~X1 + strata(X2) + X6,
data=dt, ties="breslow", x = TRUE, y = TRUE)
fit.pred <- predictCox(fit, newdata=dt[1:3], times=c(3,8), type = "survival",
se = TRUE, iid = TRUE, band = TRUE, confint = FALSE)
confint.pred1 <- confint(fit.pred, survival.transform = "none", n.sim = n.sim)
confint.pred2 <- confint(fit.pred, survival.transform = "loglog", n.sim = n.sim)
## ** standard errors
test_that("[predictCox] consistency iid se", {
expect_equal(ignore_attr=TRUE,fit.pred$survival.se[1,],
sqrt(colSums(fit.pred$survival.iid[,,1]^2))
)
})
## ** confidence intervals / bands computed on the original scale
test_that("[confint.predictCox] manual computation on ci", {
## orignial scale
expect_equal(ignore_attr=TRUE,confint.pred1$survival.lower,
fit.pred$survival + qnorm(0.025) * fit.pred$survival.se)
expect_equal(ignore_attr=TRUE,as.double(confint.pred1$survival.upper),
pmin(1,fit.pred$survival + qnorm(0.975) * fit.pred$survival.se))
## loglog scale
newse <- fit.pred$survival.se/(-fit.pred$survival*log(fit.pred$survival))
expect_equal(ignore_attr=TRUE,confint.pred2$survival.lower,
exp(-exp(log(-log(fit.pred$survival)) + qnorm(0.975) * newse)))
expect_equal(ignore_attr=TRUE,confint.pred2$survival.upper,
exp(-exp(log(-log(fit.pred$survival)) + qnorm(0.025) * newse)))
})
## * [predictCox] Diag argument
cat("[predictCox] Argument \'diag\' \n")
set.seed(10)
dt <- sampleData(5e1, outcome = "survival")[,.(time,event,X1,X2,X6)]
test_that("[predictCox] diag no strata", {
e.coxph <- coxph(Surv(time, event) ~ X1*X6, data = dt, y = TRUE, x = TRUE)
GS <- predictCox(e.coxph, newdata = dt, times = dt$time, se = TRUE, iid = TRUE, average.iid = TRUE)
test <- predictCox(e.coxph, newdata = dt, times = dt$time,
se = TRUE, iid = TRUE, average.iid = TRUE, diag = TRUE)
test2 <- predictCox(e.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = TRUE, diag = TRUE)
## estimates
expect_equal(ignore_attr=TRUE,dt$time, as.double(test$time))
expect_equal(ignore_attr=TRUE,diag(GS$cumhazard), as.double(test$cumhazard))
expect_equal(ignore_attr=TRUE,diag(GS$survival), as.double(test$survival))
## se
expect_equal(ignore_attr=TRUE,diag(GS$cumhazard.se), test$cumhazard.se[,1])
expect_equal(ignore_attr=TRUE,diag(GS$survival.se), test$survival.se[,1])
## iid
GS.iid.diag <- do.call(cbind,lapply(1:NROW(dt),
function(iN){GS$survival.iid[,iN,iN]}))
expect_equal(ignore_attr=TRUE,GS.iid.diag, test$survival.iid[,1,])
## average.iid
expect_equal(ignore_attr=TRUE,rowMeans(GS.iid.diag), test2$survival.average.iid[,1])
expect_equal(ignore_attr=TRUE,test$survival.average.iid, test2$survival.average.iid)
## average.iid with factor - diag=FALSE
average.iid <- TRUE
attr(average.iid,"factor") <- list(matrix(1:length(dt$time), nrow = NROW(dt), ncol = length(dt$time), byrow = TRUE),
matrix(1:NROW(dt), nrow = NROW(dt), ncol = length(dt$time)))
test3 <- predictCox(e.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = average.iid, diag = FALSE)
expect_equal(ignore_attr=TRUE,rowMultiply_cpp(GS$survival.average.iid, scale = 1:length(dt$time)), test3$survival.average.iid[[1]])
expect_equal(ignore_attr=TRUE,apply(GS$survival.iid, 1:2, function(x){sum(x * (1:length(dt$time)))/length(x)}),
test3$survival.average.iid[[2]])
## average.iid with factor - diag=FALSE, time varying factor
average.iid <- TRUE
attr(average.iid,"factor") <- list(matrix(rnorm(NROW(dt)*length(dt$time)), nrow = NROW(dt), ncol = length(dt$time)))
test4 <- predictCox(e.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = average.iid, diag = FALSE)
## iObs <- 1
expect_equal(ignore_attr=TRUE,do.call(rbind,lapply(1:NROW(dt), function(iObs){rowMeans(GS$survival.iid[iObs,,] * t(attr(average.iid,"factor")[[1]]))})),
test4$survival.average.iid[[1]])
## average.iid with factor - diag=TRUE
average.iid <- TRUE
attr(average.iid,"factor") <- list(matrix(5, nrow = NROW(dt), ncol = 1, byrow = TRUE),
matrix(1:NROW(dt), nrow = NROW(dt), ncol = 1))
test5 <- predictCox(e.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = average.iid, diag = TRUE)
expect_equal(ignore_attr=TRUE,5*test$survival.average.iid, test5$survival.average.iid[[1]])
expect_equal(ignore_attr=TRUE,rowMeans(rowMultiply_cpp(GS.iid.diag, 1:length(dt$time))),
test5$survival.average.iid[[2]][,1])
})
test_that("[predictCox] diag strata", {
eS.coxph <- coxph(Surv(time, event) ~ strata(X1) + X6, data = dt, y = TRUE, x = TRUE)
GS <- predictCox(eS.coxph, newdata = dt, times = dt$time, se = FALSE, iid = TRUE, average.iid = TRUE)
test <- predictCox(eS.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = TRUE, average.iid = TRUE, diag = TRUE)
test2 <- predictCox(eS.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = TRUE, diag = TRUE)
## estimate
expect_equal(ignore_attr=TRUE,dt$time, as.double(test$time))
expect_equal(ignore_attr=TRUE,diag(GS$cumhazard), as.double(test$cumhazard))
expect_equal(ignore_attr=TRUE,diag(GS$survival), as.double(test$survival))
## iid
GS.iid.diag <- do.call(cbind,lapply(1:NROW(dt),
function(iN){GS$survival.iid[,iN,iN]}))
expect_equal(ignore_attr=TRUE,GS.iid.diag, test$survival.iid[,1,])
## average.iid
expect_equal(ignore_attr=TRUE,rowMeans(GS.iid.diag), test2$survival.average.iid[,1])
expect_equal(ignore_attr=TRUE,test$survival.average.iid, test2$survival.average.iid)
## average.iid with factor - diag=FALSE
average.iid <- TRUE
attr(average.iid,"factor") <- list(matrix(1:length(dt$time), nrow = NROW(dt), ncol = length(dt$time), byrow = TRUE),
matrix(1:NROW(dt), nrow = NROW(dt), ncol = length(dt$time)))
test3 <- predictCox(eS.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = average.iid, diag = FALSE)
expect_equal(ignore_attr=TRUE,rowMultiply_cpp(GS$survival.average.iid, scale = 1:length(dt$time)), test3$survival.average.iid[[1]])
expect_equal(ignore_attr=TRUE,apply(GS$survival.iid, 1:2, function(x){sum(x * (1:length(dt$time)))/length(x)}),
test3$survival.average.iid[[2]])
## average.iid with factor - diag=FALSE, time varying factor
average.iid <- TRUE
attr(average.iid,"factor") <- list(matrix(rnorm(NROW(dt)*length(dt$time)), nrow = NROW(dt), ncol = length(dt$time)))
test4 <- predictCox(eS.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = average.iid, diag = FALSE)
expect_equal(ignore_attr=TRUE,do.call(rbind,lapply(1:NROW(dt), function(iObs){rowMeans(GS$survival.iid[iObs,,] * t(attr(average.iid,"factor")[[1]]))})),
test4$survival.average.iid[[1]])
## average.iid with factor - diag=TRUE
average.iid <- TRUE
attr(average.iid,"factor") <- list(matrix(5, nrow = NROW(dt), ncol = 1, byrow = TRUE),
matrix(1:NROW(dt), nrow = NROW(dt), ncol = 1))
test5 <- predictCox(eS.coxph, newdata = dt, times = dt$time,
se = FALSE, iid = FALSE, average.iid = average.iid, diag = TRUE)
expect_equal(ignore_attr=TRUE,5*test$survival.average.iid, test5$survival.average.iid[[1]])
expect_equal(ignore_attr=TRUE,rowMeans(rowMultiply_cpp(GS.iid.diag, 1:length(dt$time))),
test5$survival.average.iid[[2]][,1])
})
## * [predictCox] Miscellaneous
## ** Confidence bands vs timereg
cat("[predictCox] Confidence band vs timereg \n")
## *** Data
set.seed(10)
dt <- sampleData(1e2, outcome = "survival")
newdata <- dt[1:10,]
dtStrata <- data.frame(time=c(4,3,1,1,2,2,3),
status=c(1,1,1,0,1,1,0),
x=c(0,2,1,1,1,0,0),
sex=c(0,0,0,0,1,1,1))
## *** Model
e.timereg <- cox.aalen(Surv(time, event) ~ prop(X1) + prop(X2), data = dt, max.timepoint.sim=NULL)
e.coxph <- coxph(Surv(time, event) ~ X1 + X2, data = dt, x = TRUE, y = TRUE)
vec.times <- e.timereg$time.sim.resolution
## *** Compute quantile for confidence bands
resTimereg <- list()
for(i in 1:NROW(newdata)){ # i <- 1
set.seed(10)
resTimereg[[i]] <- predict.timereg(e.timereg,
newdata = newdata[i,,drop=FALSE],
times = vec.times,
resample.iid = 1,
n.sim = n.sim)
}
## *** Tests
test_that("[predictCox] Quantile for the confidence band of the cumhazard", {
predRR <- predictCox(e.coxph,
newdata = newdata,
times = vec.times,
se = TRUE,
iid = TRUE,
band = TRUE,
confint = FALSE,
type = "cumhazard")
## compatibility with timereg
ref <- unlist(lapply(resTimereg,"[[", "unif.band"))
set.seed(10)
iid2cpp <- array(NA, dim(predRR$cumhazard.iid))
for(iC in 1:dim(predRR$cumhazard.iid)[3]){ ## iC <- 1
iid2cpp[,,iC] <- rowScale_cpp(predRR$cumhazard.iid[,,iC],sqrt(diag(crossprod(predRR$cumhazard.iid[,,iC]))))
}
pred.band2 <- riskRegression:::quantileProcess_cpp(nSample = dim(predRR$cumhazard.iid)[1],
nContrast = dim(predRR$cumhazard.iid)[3],
nSim = n.sim,
iid = aperm(iid2cpp, c(2,1,3)),
alternative = 3,
global = FALSE,
confLevel = 0.95)
expect_equal(ignore_attr=TRUE,pred.band2,ref)
## note confint is removing the first column since the standard error is 0
set.seed(10)
pred.band2.no0 <- riskRegression:::quantileProcess_cpp(nSample = dim(predRR$cumhazard.iid)[1],
nContrast = dim(predRR$cumhazard.iid)[3],
nSim = n.sim,
iid = aperm(iid2cpp[,-1,], c(2,1,3)),
alternative = 3,
global = FALSE,
confLevel = 0.95)
## should not set transform to NA because at time 0 se=0 so the log-transform fails
pred.confint <- confint(predRR, n.sim = n.sim, seed = 10,
cumhazard.transform = "none")
expect_equal(ignore_attr=TRUE,pred.confint$cumhazard.quantileBand, pred.band2.no0)
expect_equal(ignore_attr=TRUE,pred.confint$cumhazard.quantileBand, ref)
})
## *** Display
predRR <- predictCox(e.coxph,
newdata = newdata[1],
times = vec.times,
se = TRUE,
band = TRUE,
type = c("cumhazard","survival")
)
## plotRR <- autoplot(predRR, type = "survival", band = TRUE, ci = TRUE, plot = FALSE)
## dev.new()
## plotTR <- plot.predict.timereg(resTimereg[[1]])
## dev.new()
## plotRR$plot + coord_cartesian(ylim = c(0,1))
## graphics.off()
## *** With strata
## Fit a stratified model
eS.coxph <- coxph(Surv(time, status) ~ x + strata(sex),
data = dtStrata, x = TRUE, y = TRUE)
eS.pred <- predictCox(eS.coxph, newdata = dtStrata, times = 1:4,
se = TRUE, iid = TRUE, band = TRUE)
eS.confint <- confint(eS.pred, seed = 10)
# eS.confint$survival.quantileBand
## ** Store.iid argument
cat("[predictCox] Check same result store.iid=minimal vs. full \n")
## *** Data
set.seed(10)
d <- sampleData(50, outcome = "survival")
setkey(d,time)
## *** no strata
m.coxph <- coxph(Surv(time, event) ~ X1+X6, data = d, y = TRUE, x = TRUE)
seqTime <- c(1e-16,4:10,d$time[1:10],1e6)
newdata <- d
## system.time(
## res1 <- predictCox(m.coxph, times = seqTime, newdata = newdata, store.iid = "minimal", se = TRUE, iid = FALSE)
## )
## system.time(
## res3 <- predictCox(m.coxph, times = seqTime, newdata = newdata, store.iid = "full", se = TRUE, iid = FALSE)
## )
test_that("[predictCox] store.iid = minimal vs. full - no strata", {
res1 <- predictCox(m.coxph, times = seqTime, newdata = newdata,
type = c("cumhazard", "survival"),
store.iid = "minimal", se = TRUE, iid = TRUE, average.iid = TRUE)
res2 <- predictCox(m.coxph, times = seqTime, newdata = newdata,
type = c("cumhazard", "survival"),
store.iid = "full", se = TRUE, iid = TRUE)
expect_equal(ignore_attr=TRUE,res1$cumhazard.se,res2$cumhazard.se)
expect_equal(ignore_attr=TRUE,res1$survival.se,res2$survival.se)
expect_equal(ignore_attr=TRUE,res1$cumhazard.iid,res2$cumhazard.iid)
expect_equal(ignore_attr=TRUE,res1$survival.iid,res2$survival.iid)
expect_equal(ignore_attr=TRUE,res1$cumhazard.average.iid, apply(res2$cumhazard.iid,1:2,mean))
expect_equal(ignore_attr=TRUE,res1$survival.average.iid, apply(res2$survival.iid,1:2,mean))
## pre-store
m2.coxph <- iidCox(m.coxph, store.iid = "minimal")
res1bis <- predictCox(m2.coxph, times = seqTime, newdata = newdata,
type = c("cumhazard", "survival"),
se = TRUE, iid = TRUE, average.iid = TRUE)
expect_equal(ignore_attr=TRUE,res1bis$cumhazard.se,res2$cumhazard.se)
expect_equal(ignore_attr=TRUE,res1bis$survival.se,res2$survival.se)
expect_equal(ignore_attr=TRUE,res1bis$cumhazard.iid,res2$cumhazard.iid)
expect_equal(ignore_attr=TRUE,res1bis$survival.iid,res2$survival.iid)
expect_equal(ignore_attr=TRUE,res1bis$cumhazard.average.iid, apply(res2$cumhazard.iid,1:2,mean))
expect_equal(ignore_attr=TRUE,res1bis$survival.average.iid, apply(res2$survival.iid,1:2,mean))
})
## *** strata
m.coxph <- coxph(Surv(time, event) ~ strata(X1)+X6, data = d, y = TRUE, x = TRUE)
seqTime <- c(1e-16,4:10,d$time[1:10],1e6)
newdata <- d
test_that("[predictCox] store.iid = minimal vs. full - strata", {
newdata <- rbind(d[1],d[1])
res1 <- predictCox(m.coxph, times = seqTime, newdata = newdata,
type = c("cumhazard", "survival"),
store.iid = "minimal", se = TRUE, iid = TRUE, average.iid = TRUE)
res2 <- predictCox(m.coxph, times = seqTime, newdata = newdata,
type = c("cumhazard", "survival"),
store.iid = "full", se = TRUE, iid = TRUE)
expect_equal(ignore_attr=TRUE,res1$cumhazard.se,res2$cumhazard.se)
expect_equal(ignore_attr=TRUE,res1$survival.se,res2$survival.se)
expect_equal(ignore_attr=TRUE,res1$cumhazard.iid,res2$cumhazard.iid)
expect_equal(ignore_attr=TRUE,res1$survival.iid,res2$survival.iid)
expect_equal(ignore_attr=TRUE,res1$cumhazard.average.iid, apply(res2$cumhazard.iid,1:2,mean))
expect_equal(ignore_attr=TRUE,res1$survival.average.iid, apply(res2$survival.iid,1:2,mean))
## pre store
m2.coxph <- iidCox(m.coxph, store.iid = "minimal")
res1bis <- predictCox(m2.coxph, times = seqTime, newdata = newdata,
type = c("cumhazard", "survival"),
se = TRUE, iid = TRUE, average.iid = TRUE)
expect_equal(ignore_attr=TRUE,res1bis$cumhazard.se,res2$cumhazard.se)
expect_equal(ignore_attr=TRUE,res1bis$survival.se,res2$survival.se)
expect_equal(ignore_attr=TRUE,res1bis$cumhazard.iid,res2$cumhazard.iid)
expect_equal(ignore_attr=TRUE,res1bis$survival.iid,res2$survival.iid)
expect_equal(ignore_attr=TRUE,res1bis$cumhazard.average.iid, apply(res2$cumhazard.iid,1:2,mean))
expect_equal(ignore_attr=TRUE,res1bis$survival.average.iid, apply(res2$survival.iid,1:2,mean))
})
## ** Weigthed cox
cat("[predictCox] Does not handle weights \n")
## *** Data
set.seed(10)
data(Melanoma, package = "riskRegression")
wdata <- runif(nrow(Melanoma), 0, 1)
times1 <- unique(Melanoma$time)
## *** Test
test_that("[predictCox] - weights",{
fitW.coxph <- coxph(Surv(time,status == 1) ~ thick*age, data = Melanoma, weights = wdata, y = TRUE, x = TRUE)
fitW.cph <- cph(Surv(time,status == 1) ~ thick*age, data = Melanoma, y = TRUE, x = TRUE, weights = wdata)
expect_error(resW <- predictCox(fitW.coxph, times = Melanoma$time, newdata = Melanoma))
expect_error(resW <- predictCox(fitW.cph, times = Melanoma$time, newdata = Melanoma))
})
## ** Deal with negative timepoints
cat("[predictCox] dealing with negative timepoints or NA \n")
data(Melanoma, package = "riskRegression")
fit.coxph <- coxph(Surv(time,status == 1) ~ thick*age, data = Melanoma, y = TRUE, x = TRUE)
test_that("Deal with negative/NA time points",{
expect_equal(ignore_attr=TRUE,unname(predictCox(fit.coxph, times = -1, newdata = Melanoma)$survival),
matrix(1,nrow = NROW(Melanoma), ncol = 1))
expect_error(predictCox(fit.coxph, times = c(1,2,NA), newdata = Melanoma))
})
# }}}
## ** Deal with no event
set.seed(10)
dt <- sampleData(1e2, outcome = "survival")
dt$event <- 0
e.coxph <- coxph(Surv(time,event) ~ 1, data = dt, x = TRUE, y = TRUE)
test_that("Deal with no event",{
expect_true(all(predictCox(e.coxph)$cumhazard==0))
expect_true(all(predictCox(e.coxph)$survival==1))
expect_true(all(predictCox(e.coxph, newdata = dt, times = 1:5)$cumhazard==0))
expect_true(all(predictCox(e.coxph, newdata = dt, times = 1:5)$survival==1))
})
## ** Fully stratified model and NAs
set.seed(10)
d <- sampleData(1e2, outcome = "survival")
setkeyv(d, "time")
d[X1==0,event := c(event[1:(.N-1)],0)]
d[X1==1,event := c(event[1:(.N-1)],1)]
tau <- c(d[,max(time),by="X1"][[2]],10000)
test_that("After last event - fully stratified model",{
## one strata variable
X <- unique(d[,"X1",drop=FALSE])
e.coxph <- coxph(Surv(time,event) ~ strata(X1), data = d, x = TRUE, y = TRUE)
test <- as.data.table(predictCox(e.coxph, times = tau, newdata = X, se = TRUE))
expect_equal(ignore_attr=TRUE,test[strata == 1 & times == 10000,survival], test[strata == 1 & times == d[X1==1,max(time)],survival], tolerance = 1e-6)
expect_equal(ignore_attr=TRUE,test[strata == 1 & times == 10000,survival.se], test[strata == 1 & times == d[X1==1,max(time)],survival.se], tolerance = 1e-6)
expect_true(is.na(test[strata == 0 & times == 10000,survival]))
expect_true(is.na(test[strata == 0 & times == 10000,survival.se]))
test2 <- as.data.table(predictCoxPL(e.coxph, times = tau, newdata = X, se = TRUE))
expect_equal(ignore_attr=TRUE,test2[strata == 1 & times == 10000,survival], 0, tolerance = 1e-6)
expect_equal(ignore_attr=TRUE,test2[strata == 1 & times == 10000,survival.se], test2[strata == 1 & times == d[X1==1,max(time)],survival.se], tolerance = 1e-6)
expect_true(is.na(test2[strata == 0 & times == 10000,survival]))
expect_true(is.na(test2[strata == 0 & times == 10000,survival.se]))
## two strata variables
X <- unique(d[,c("X1","X2"),drop=FALSE])
e2.coxph <- coxph(Surv(time,event) ~ strata(X1)+strata(X2), data = d, x = TRUE, y = TRUE)
test <- as.data.table(predictCox(e2.coxph, times = tau, newdata = X, se = TRUE))
expect_equal(ignore_attr=TRUE,test[strata == "1, 0" & times == 10000,survival], test[strata == "1, 0" & times == d[X1==1 & X2 == 0,max(time)],survival])
expect_equal(ignore_attr=TRUE,test[strata == "1, 0" & times == 10000,survival.se], test[strata == "1, 0" & times == d[X1==1 & X2 == 0,max(time)],survival.se])
expect_equal(ignore_attr=TRUE,test[strata == "1, 0" & times == 10000,survival], test[strata == "1, 0" & times == d[X1==1 & X2 == 0,max(time)],survival])
expect_equal(ignore_attr=TRUE,test[strata == "1, 0" & times == 10000,survival.se], test[strata == "1, 0" & times == d[X1==1 & X2 == 0,max(time)],survival.se])
expect_equal(ignore_attr=TRUE,test[strata == "1, 1" & times == 10000,survival], test[strata == "1, 1" & times == d[X1==1 & X2 == 0,max(time)],survival])
expect_equal(ignore_attr=TRUE,test[strata == "1, 1" & times == 10000,survival.se], test[strata == "1, 1" & times == d[X1==1 & X2 == 0,max(time)],survival.se])
expect_true(all(is.na(test[strata %in% c("0, 0","0, 1") & times == 10000,survival])))
expect_true(all(is.na(test[strata %in% c("0, 0","0, 1") & times == 10000,survival.se])))
test2 <- as.data.table(predictCoxPL(e2.coxph, times = tau, newdata = X, se = TRUE))
expect_true(all(test2[strata %in% c("1, 0","1, 1") & times == 10000,survival]==0))
expect_true(all(!is.na(test2[strata %in% c("1, 0","1, 1") & times == 10000,survival.se])))
expect_true(all(is.na(test2[strata %in% c("0, 0","0, 1") & times == 10000,survival])))
expect_true(all(is.na(test2[strata %in% c("0, 0","0, 1") & times == 10000,survival.se])))
})
## * [predictCox] Previous Bug
cat("[predictCox] Previous bug \n")
## ** Some coef are NA
set.seed(10)
dt <- sampleData(5e2, outcome = "survival")
e.coxph <- coxph(Surv(time, event) ~ X1+ X6 , data = dt, y = TRUE, x = TRUE)
e.coxph$coefficients[] <- as.numeric(NA)
test_that("Return error when coef contains NA", {
expect_error(capture.output(predictCox(e.coxph, newdata = dt, times = 1)))
})
## ** average.iid
test_that("Cox - output of average.iid should not depend on other arguments", {
set.seed(10)
d <- sampleData(70,outcome="survival")
d[, X1 := paste0("T",rbinom(.N, size = 2, prob = c(0.51)))]
fit <- coxph(Surv(time,event)~X1 + strata(X2) + X6,
data=d, ties="breslow", x = TRUE, y = TRUE)
out1 <- predictCox(fit, newdata = d[1:5], times = 1:3, average.iid = TRUE)
out2 <- predictCox(fit, newdata = d[1:5], times = 1:3, se = TRUE, average.iid = TRUE)
expect_equal(ignore_attr=TRUE,out1$survival.average.iid,out2$survival.average.iid, tolerance = 1e-8)
})
## ** Incorrect calculation of the standard error with atanh (i.e. se/(1+b^2) instead of se/(1-b^2))
## from: Paul Blanche <pabl@sund.ku.dk>
## subject: suspected error in riskRegression
## date: Tue, 30 Jul 2019 11:42:14 +0200
test_that("Standard error after atanh transformation", {
set.seed(10)
n <- 1e2
x <- rnorm(n)
y <- rnorm(n)
e.test <- cor.test(x,y)
rho <- e.test$estimate
rho.se <- (1-rho^2)
e.trans <- transformCIBP(estimate = 0.5, se = cbind(1), null = 0,
type = "atanh", ci = TRUE, conf.level = 0.95, alternative = "two.sided",
min.value = NULL, max.value = NULL, band = FALSE,
p.value = TRUE, seed = 10, df = Inf)
expect_equal(ignore_attr=TRUE,abs(atanh(0.5)-atanh(e.trans$lower)),abs(atanh(0.5)-atanh(e.trans$upper)), tolerance = 1e-6)
expect_equal(ignore_attr=TRUE,as.double((atanh(0.5)-atanh(e.trans$lower)))/qnorm(0.975),1/(1-0.5^2), tolerance = 1e-6)
})
## ** se/iid should not depend on the ordering of the argument times
test_that("Cox - iid/se should not depend on other arguments", {
set.seed(10)
d <- sampleData(70,outcome="survival")
d[, X1 := paste0("T",rbinom(.N, size = 2, prob = c(0.51)))]
fit <- coxph(Surv(time,event)~X1 + strata(X2) + X6,
data=d, ties="breslow", x = TRUE, y = TRUE)
seqTau <- abs(rnorm(10))
out1 <- predictCox(fit, newdata = d[1:5], times = seqTau,
se = TRUE, iid = TRUE, average.iid = TRUE)
out2 <- predictCox(fit, newdata = d[1:5], times = sort(seqTau),
se = TRUE, iid = TRUE, average.iid = TRUE)
out3 <- predictCox(fit, newdata = d[1:5], times = seqTau,
average.iid = TRUE)
out4 <- predictCox(fit, newdata = d[1:5], times = sort(seqTau),
average.iid = TRUE)
out5 <- predictCox(fit, newdata = d[1:5], times = seqTau,
se = TRUE, iid = TRUE, average.iid = TRUE, store.iid = "minimal")
out6 <- predictCox(fit, newdata = d[1:5], times = sort(seqTau),
se = TRUE, iid = TRUE, average.iid = TRUE, store.iid = "minimal")
expect_equal(ignore_attr=TRUE,out1$survival.iid[,order(seqTau),],out2$survival.iid)
expect_equal(ignore_attr=TRUE,out1$survival.se[,order(seqTau)],out2$survival.se)
expect_equal(ignore_attr=TRUE,out1$survival.average.iid[,order(seqTau)],out2$survival.average.iid)
expect_equal(ignore_attr=TRUE,out2$survival.average.iid,out4$survival.average.iid)
expect_equal(ignore_attr=TRUE,out3$survival.average.iid[,order(seqTau)],out4$survival.average.iid)
expect_equal(ignore_attr=TRUE,out5$survival.iid[,order(seqTau),],out6$survival.iid)
expect_equal(ignore_attr=TRUE,out5$survival.se[,order(seqTau)],out6$survival.se)
expect_equal(ignore_attr=TRUE,out5$survival.average.iid[,order(seqTau)],out6$survival.average.iid)
expect_equal(ignore_attr=TRUE,out2$survival.iid,out6$survival.iid)
expect_equal(ignore_attr=TRUE,out2$survival.se,out6$survival.se)
expect_equal(ignore_attr=TRUE,out2$survival.average.iid,out6$survival.average.iid)
})
## ** iidCox is working when stopped early
test_that("iidCox - stopped early", {
n <- 500
set.seed(10)
dt <- sampleData(n, outcome="survival")
seqJump <- sort(dt[dt$event==1,time])
m.cox <- coxph(Surv(time,event)~ X1*X6+strata(X2),data=dt, x = TRUE, y = TRUE)
m2.cox <- iidCox(m.cox, tau.max = seqJump[5], store.iid = "minimal")
m3.cox <- iidCox(m.cox, tau.max = seqJump[5], store.iid = "full")
GS <- predictCox(m.cox, newdata = dt, times = seqJump[1:5], average.iid = TRUE)
test1 <- predictCox(m2.cox, newdata = dt, times = seqJump[1:5], average.iid = TRUE)
test2 <- predictCox(m3.cox, newdata = dt, times = seqJump[1:5], average.iid = TRUE)
expect_equal(ignore_attr=TRUE,GS$survival.average.iid,test1$survival.average.iid)
expect_equal(ignore_attr=TRUE,GS$survival.average.iid,test2$survival.average.iid)
})
## ** predictCox with interacitons
test_that("predictCox - : operator", {
n <- 500
set.seed(10)
dt <- sampleData(n, outcome="survival")
m.cox <- coxph(Surv(time,event)~ X1*X6+strata(X2),data=dt, x = TRUE, y = TRUE)
test <- predictCox(m.cox, newdata = dt[1], times = 5:8, average.iid = TRUE)
m2.cox <- coxph(Surv(time,event)~ X1:X6+strata(X2),data=dt, x = TRUE, y = TRUE)
test2 <- predictCox(m2.cox, newdata = dt[1], times = 5:8, average.iid = TRUE) ## previously an error because of missing main term
expect_equal(as.double(test2$survival),c(0.5786438,0.5475801,0.4452037,0.3406072), tol = 1e-6)
})
#----------------------------------------------------------------------
### test-predictCox.R ends here
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.