library(testthat)
context("Cause-specific Cox regression")
library(riskRegression)
library(rms)
library(survival)
library(prodlim)
test_that("CSC vs prodlim",{
data(Melanoma)
nd <- data.frame(sex=factor(levels(Melanoma$sex)))
A <- prodlim(Hist(time,status)~1,data=Melanoma)
B <- CSC(Hist(time,status)~1,data=Melanoma)
pB <- predictRisk(B,times=sort(unique(Melanoma$time)),newdata=nd[1,,drop=FALSE],cause=1)
pA <- predictRisk(A,times=sort(unique(Melanoma$time)),newdata=nd[1,,drop=FALSE],cause=1)
expect_equal(ignore_attr=TRUE,c(pB),pA,tolerance=1e3)
a <- prodlim(Hist(time,status)~sex,data=Melanoma)
b <- CSC(Hist(time,status)~strat(sex),data=Melanoma,fitter="cph")
c <- CSC(Hist(time,status)~strat(sex),data=Melanoma,surv.type="survival",fitter="cph")
pa <- predictRisk(a,times=c(0,10,100,1000,2000),newdata=nd,cause=1)
pb <- predictRisk(b,times=c(0,10,100,1000,2000),newdata=nd,cause=1)
pc <- predictRisk(c,times=c(0,10,100,1000,2000),newdata=nd,cause=1)
expect_equal(ignore_attr=TRUE,c(pb),c(pa),tolerance=0.1)
expect_equal(ignore_attr=TRUE,c(pb),c(pc),tolerance=0.00000001)
# pd <- predictEventProb(c,times=c(0,10,100,1000,2000),newdata=nd,cause=1)
# expect_equal(ignore_attr=TRUE,c(pc),c(pd),tolerance=0.00000001)
u <- CSC(Hist(time,status)~strat(sex)+age+invasion,data=Melanoma,fitter="cph")
v <- CSC(Hist(time,status)~strat(sex)+age+invasion,data=Melanoma,surv.type="survival",fitter="cph")
pu <- predictRisk(u,times=c(0,10,100,1000,2000),newdata=Melanoma[c(17,84),],cause=1)
pv <- predictRisk(v,times=c(0,10,100,1000,2000),newdata=Melanoma[c(17,84),],cause=1)
expect_equal(ignore_attr=TRUE,c(pu),c(pv),tolerance=0.1)
## plot(a,cause=1,ylim=c(0,.5),confint=FALSE)
## lines(sort(unique(Melanoma$time)),pb[1,],lty=3)
## lines(sort(unique(Melanoma$time)),pb[2,],lty=3,col=2)
})
if (requireNamespace("pec",quietly=TRUE)){
library(pec)
test_that("predictSurv",{
set.seed(17)
d <- prodlim::SimSurv(100)
f <- coxph(Surv(time,status)~X1+X2,data=d,x=TRUE,y=TRUE)
h <- cph(Surv(time,status)~X1+X2,data=d,surv=TRUE,x=TRUE,y=TRUE)
af <- predictRisk(f,newdata=d[c(17,88,3),],times=c(0,1,8.423,100,1000))
bf <- 1-predictSurvProb(f,newdata=d[c(17,88,3),],times=c(0,1,8.423,100,1000))
expect_equal(ignore_attr=TRUE,unname(af),unname(bf),tolerance = 1e-8)
ah <- predictRisk(h,newdata=d[c(17,88,3),],times=c(0,1,8.423,100,1000))
bh <- 1-predictSurvProb(h,newdata=d[c(17,88,3),],times=c(0,1,8.423,100,1000))
colnames(bh) <- NULL
expect_equal(ignore_attr=TRUE,unname(ah),unname(bh),tolerance = 1e-8)
})
}
test_that("Cox models",{
set.seed(17)
d <- prodlim::SimCompRisk(100)
a <- CSC(Hist(time,event)~X1+X2,data=d)
A <- CSC(Hist(time,event)~X1+X2,data=d,surv.type="surv")
a1 <- coxph(Surv(time,event==1)~X1+X2,data=d)
a2 <- coxph(Surv(time,event==2)~X1+X2,data=d)
A2 <- coxph(Surv(time,event!=0)~X1+X2,data=d)
expect_equal(ignore_attr=TRUE,coef(a$models[[1]]),coef(a1),tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,coef(a$models[[2]]),coef(a2),tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,coef(A$models[[2]]),coef(A2),tolerance = 1e-8)
})
test_that("strata",{
set.seed(17)
d <- prodlim::SimCompRisk(100)
a <- CSC(Hist(time,event)~strata(X1)+X2,data=d)
A <- CSC(Hist(time,event)~strata(X1)+X2,data=d,surv.type="surv")
a1 <- coxph(Surv(time,event==1)~strata(X1)+X2,data=d)
a2 <- coxph(Surv(time,event==2)~strata(X1)+X2,data=d)
A2 <- coxph(Surv(time,event!=0)~strata(X1)+X2,data=d)
expect_equal(ignore_attr=TRUE,coef(a$models[[1]]),coef(a1),tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,coef(a$models[[2]]),coef(a2),tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,coef(A$models[[2]]),coef(A2),tolerance = 1e-8)
})
test_that("strat and strata",{
data(Melanoma)
a <- CSC(Hist(time,status)~strat(sex)+age+invasion+logthick+strat(epicel)+strat(ulcer),data=Melanoma,fitter="cph")
predictRisk(a,times=c(0,100,1000,4000),newdata=Melanoma[c(17,77,188),],cause=2)
b <- CSC(Hist(time,status)~strata(sex)+age+invasion+logthick+strata(epicel)+strata(ulcer),data=Melanoma,fitter="coxph")
pa <- predictRisk(a,times=c(0,100,1000,4000),newdata=Melanoma[c(17,77,188),],cause=2)
pb <- predictRisk(b,times=c(0,100,1000,4000),newdata=Melanoma[c(17,77,188),],cause=2)
expect_equal(ignore_attr=TRUE,pa,pb,tolerance=1e-6)
})
# test_that("CSC many character valued causes",{
# set.seed(17)
# d <- prodlim::SimCompRisk(100)
# d$event <- as.character(factor(d$event,labels=c("a","b","c")))
# m1 <- CSC(Hist(time,event)~strat(X1)+X2,data=d,fitter="cph")
# m2 <- CSC(Hist(time,event)~strata(X1)+X2,data=d,surv.type="surv",cause="b")
# expect_equal(ignore_attr=TRUE,round(coef(m1$models[[2]])[[1]],6),0.535059)
# })
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.