### test-predictRisk-TD.R ---
#----------------------------------------------------------------------
## Author: Thomas Alexander Gerds
## Created: Sep 2 2017 (08:01)
## Version:
## Last-Updated: Sep 17 2022 (07:01)
## By: Thomas Alexander Gerds
## Update #: 15
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
library(riskRegression)
library(testthat)
library(timereg)
library(survival)
test_that("start stop Cox: basehaz vs timereg",{
## remove ties
set.seed(17)
veteran$jtime <- veteran$time+rnorm(NROW(veteran),mean=0,sd=0.001)
cox <- coxph(Surv(jtime, status) ~ celltype+karno + age + trt, data= veteran,x=1)
ca <- cox.aalen(Surv(jtime, status) ~ prop(celltype)+prop(karno) + prop(age) + prop(trt), data= veteran)
bh.cox <- basehaz(cox,centered=FALSE)
## timereg's cum does not include censored times where the hazard does not change
expect_equal(ignore_attr=TRUE,bh.cox[!duplicated(bh.cox$hazard),]$time,ca$cum[-1,"time"])
expect_equal(ignore_attr=TRUE,bh.cox[!duplicated(bh.cox$hazard),"hazard"],ca$cum[-1,"(Intercept)"])
test.times <- c(0,4,300,8,veteran$jtime[c(17,10)])
expect_equal(ignore_attr=TRUE,1-predictRisk(cox,newdata=veteran[c(47,8,17,10),],times=test.times),
predict(ca,se.fit=FALSE,newdata=veteran[c(47,8,17,10),],times=test.times)$S0)
## time-dependent
vet2 <- survSplit(Surv(jtime, status) ~., veteran, cut=c(60, 120), episode ="timegroup")
coxTD <- coxph(Surv(tstart, jtime, status) ~ celltype+karno + age + trt, data= vet2,x=1)
caTD <- cox.aalen(Surv(tstart, jtime, status) ~ prop(celltype)+prop(karno) + prop(age) + prop(trt), data= vet2)
bh.coxTD <- basehaz(coxTD,centered=FALSE)
## timereg's cum does not include censored times where the hazard does not change
expect_equal(ignore_attr=TRUE,bh.coxTD[!duplicated(bh.coxTD$hazard),"hazard"],caTD$cum[-1,"(Intercept)"])
## with(basehaz(coxTD,centered=FALSE),plot(time,hazard,type="s"))
## lines(caTD$cum[,"time"],caTD$cum[,"(Intercept)"],col=2,type="s")
expect_gt(unlist(riskRegression:::predictRisk.coxphTD(coxTD,newdata=vet2[c(17),],landmark=60,times=48)),
unlist(riskRegression:::predictRisk.coxphTD(coxTD,newdata=vet2[c(17),],landmark=60,times=18)))
})
test_that("start stop CauseSpecificCox",{
library(data.table)
set.seed(9)
d1 <- sampleData(201,outcome="competing.risk")
d1[,start:=0]
fit1 <- CSC(list(Hist(time=time,event=event,entry=start)~X1+X4+log(X7)+X9),data=d1)
## expect_equal(ignore_attr=TRUE,a <- riskRegression:::predictRisk.CSCTD(fit1,newdata=d1[c(17:32),],landmark=2:5,times=1,cause=1),
## b <- riskRegression:::predictRisk.CauseSpecificCox(fit1,newdata=d1[c(17:32),],landmark=2:5,times=1,cause=1))
d2 <- sampleData(d1[,sum(event==0)],outcome="competing.risk")
d2[,start:=d1[event==0,time]]
d2[,time:=time+d1[event==0,time]]
d3 <- sampleData(sum(d2$event==0),outcome="competing.risk")
d3[,start:=d2[event==0,time]]
d3[,time:=time+d2[event==0,time]]
d <- rbindlist(list(d1,d2,d3))
## d[start>time,.(start,time,event)]
fit <- CSC(list(Hist(time=time,event=event,entry=start)~X1+X4+log(X7)+X9),data=d)
## NOTE: expect_error with regexp=NA tests no error
expect_error(riskRegression:::predictRisk.CSCTD(fit,newdata=d[c(17:32),],landmark=2:5,times=1), regexp = NA)
})
######################################################################
### test-predictRisk-TD.R ends here
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.