Running ‘test-AUC.R’ [11s/12s]
  Running ‘test-BinomialRegression.R’ [14s/16s]
  Running ‘test-BrierScore.R’ [13s/15s]
  Running ‘test-CauseSpecificCoxRegression.R’ [11s/12s]
  Running ‘test-FineGrayRegression.R’ [10s/11s]
  Running ‘test-ate.R’ [17s/18s]
  Running ‘test-boxplotScore.R’ [10s/10s]
  Running ‘test-iidCox.R’ [14s/14s]
  Running ‘test-phreg.R’ [10s/10s]
  Running ‘test-predictCSC.R’ [17s/19s]
  Running ‘test-predictCox.R’ [12s/13s]
  Running ‘test-predictCoxPL.R’ [11s/11s]
  Running ‘test-predictRisk-TD.R’ [13s/13s]
  Running ‘test-predictRisk.R’ [10s/10s]
  Running ‘test-score-bootstrap.R’ [36s/36s]
  Running ‘test-score.R’ [16s/16s]
  Running ‘test-summary-score.R’ [12s/12s]
Running the tests in ‘tests/test-BinomialRegression.R’ failed.
Complete output:
  > library(testthat)
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(timereg)
  Loading required package: survival
  > library(prodlim)
  > library(survival)
  > context("Binomial regression")
  > data(Melanoma)
  > # {{{ "Absolute risk regression"
  > test_that("Absolute risk regression",{
  +     set.seed(17)
  +     d <- sampleData(300,outcome="competing.risks")
  +     f1 <- ARR(Hist(time,event)~X1+X2+strata(X6),data=d,cause=1)
  +     f2 <- timereg::comp.risk(Event(time,event)~ const(X1)+const(X2)+X6,data=d,cause=1,model="rcif")
  +     f3 <- timereg::comp.risk(Event(time,event)~ const(X1)+const(X2)+X6,data=d,cause=1,model="rcif2")
  +     f4 <- ARR(Hist(time,event)~X1+X2+strata(X8),data=d,cause=1)    
  +     f5 <- timereg::comp.risk(Event(time,event)~ const(X1)+const(X2)+X8,data=d,cause=1,model="rcif")
  +     f6 <- timereg::comp.risk(Event(time,event)~ const(X1)+const(X2)+X8,data=d,cause=1,model="rcif2")
  +     expect_equal(as.numeric(f1$timeConstantEffects$coef),c(f3$gamma))
  +     expect_equal(as.numeric(f4$timeConstantEffects$coef),c(f6$gamma))
  +     d[,X4:=factor(X4)]
  +     system.time(A <- ARR(Hist(time,event)~X1+X3+X4,data=d,cause=1))
  +     ## system.time(B <- timereg::comp.risk(Event(time,event)~ const(factor(X1))+ const(factor(X3)),data=d,cause=1,model="rcif"))
  +     system.time(B <- timereg::comp.risk(Event(time,event)~ const(X1)+ const(X3)+const(X4),data=d,cause=1,model="rcif2"))
  +     ## head(A$timeVaryingEffects$coef)
  +     ## head(B$cum)
  +     expect_equal(as.numeric(A$timeConstantEffects$coef),c(B$gamma))    
  + })
  Test passed 🥇
  > 
  > # }}}
  > # {{{ "Logistic risk regression"
  > test_that("Logistic risk regression",{
  +     set.seed(17)
  +     d <- prodlim::SimCompRisk(100)
  +     a <- LRR(Hist(time,event)~X1+X2,data=d,cause=1)
  +     b <- timereg::comp.risk(Event(time,event)~ const(X1)+const(X2),data=d,cause=1,model="logistic2")
  +     expect_equal(as.numeric(a$timeConstantEffects$coef),c(b$gamma))
  +     A <- LRR(Hist(time,event)~X1+strata(X2),data=d,cause=1)
  +     B <- timereg::comp.risk(Event(time,event)~ const(X1)+X2,data=d,cause=1,model="logistic2")
  +     ## head(A$timeVaryingEffects$coef)
  +     ## head(B$cum)
  +     expect_equal(as.numeric(A$timeConstantEffects$coef),c(B$gamma))    
  + })
  Test passed 🥳
  > 
  > # }}}
  > # {{{ "Censoring model"
  > test_that("Censoring model",{
  +     expect_success(ARR(Hist(time,status)~thick+strata(invasion)+epicel,data=Melanoma,cens.model="cox",cause=1))
  +     expect_success(ARR(Hist(time,status)~thick+strata(invasion)+epicel,data=Melanoma,cens.model="cox",cens.formula=~thick+strat(invasion)+epicel,cause=1))
  +     expect_success(ARR(Hist(time,status)~sex+strata(invasion)+epicel,data=Melanoma,cens.model="cox",cens.formula=~logthick+age,cause=1))
  + })
  ── Failure (???): Censoring model ──────────────────────────────────────────────
  no expectation used.
  
  ── Failure (???): Censoring model ──────────────────────────────────────────────
  no expectation used.
  
  ── Failure (???): Censoring model ──────────────────────────────────────────────
  no expectation used.
  
  Error: Test failed
  Execution halted
Running the tests in ‘tests/test-CauseSpecificCoxRegression.R’ failed.
Complete output:
  > library(testthat)
  > context("Cause-specific Cox regression")
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(rms)
  Loading required package: Hmisc
  Loading required package: lattice
  Loading required package: survival
  Loading required package: Formula
  Loading required package: ggplot2
  
  Attaching package: 'Hmisc'
  
  The following object is masked from 'package:testthat':
  
      describe
  
  The following objects are masked from 'package:base':
  
      format.pval, units
  
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > 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(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(c(pb),c(pa),tolerance=0.1)
  +     expect_equal(c(pb),c(pc),tolerance=0.00000001)
  + #     pd <- predictEventProb(c,times=c(0,10,100,1000,2000),newdata=nd,cause=1)
  + #     expect_equal(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(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)
  + })
  Test passed 🌈
  > 
  > 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(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(unname(ah),unname(bh),tolerance = 1e-8)
  + })
  ── Error (???): predictSurv ────────────────────────────────────────────────────
  Error in `predictSurvProb(f, newdata = d[c(17, 88, 3), ], times = c(0, 
      1, 8.423, 100, 1000))`: could not find function "predictSurvProb"
  
  Error: Test failed
  Execution halted
Running the tests in ‘tests/test-ate.R’ failed.
Complete output:
  > ## minor change compare to the hard coded version can be due to a previous bug in iidCox (the value have not been updated)
  > ## (the calculation of the iid was stopped on jump to early when using tau.max)
  > library(testthat)
  > library(data.table)
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(survival)
  > library(rms)
  Loading required package: Hmisc
  Loading required package: lattice
  Loading required package: Formula
  Loading required package: ggplot2
  
  Attaching package: 'Hmisc'
  
  The following object is masked from 'package:testthat':
  
      describe
  
  The following objects are masked from 'package:base':
  
      format.pval, units
  
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > 
  > verbose <- FALSE
  > calcIterm <- function(factor, iid, indexJump, iid.outsideI){
  +     n <- length(indexJump)
  + 
  +     if(iid.outsideI){ ## compute iid int factor
  +         ls.I <- lapply(1:n, function(iObs){ ## iObs <- 2
  +             iIID <- iid[,1,iObs]
  +             iFactor <- factor[iObs,1:indexJump[iObs]]
  +         
  +             if(indexJump[iObs]==0){
  +                 return(rep(0,n))
  +             }else {
  +                 return(iIID*sum(iFactor))
  +             }
  +         })
  +     } else { ## compute int iid * factor
  +         ls.I <- lapply(1:n, function(iObs){ ## iObs <- 2
  +             iIID <- iid[,1:indexJump[iObs],iObs]
  +             iFactor <- factor[iObs,1:indexJump[iObs]]
  +         
  +             if(indexJump[iObs]==0){
  +                 return(rep(0,n))
  +             }else if(indexJump[iObs]==1){
  +                 return(iIID*iFactor)
  +             }else{
  +                 return(rowSums(rowMultiply_cpp(iIID, iFactor)))
  +             }
  +         })
  +     }
  +     return(do.call(cbind,ls.I))
  + }
  > 
  > 
  > ## ** Ate using list of formulae
  > test_that("[ate] using list of formulae",{
  +     set.seed(10)
  +     d <- sampleData(100)
  +     expect_output(print(ate(list(Hist(time,event)~X1+X2,
  +                             Hist(time,event)~X1+X2),
  +                        treatment = "X1", data = d, times = 1:5, cause = 1,
  +                        verbose = FALSE)))
  + })
  Test passed 😸
  > 
  > ## * [ate] G-formula only
  > cat("[ate] G-formula only \n")
  [ate] G-formula only 
  > 
  > ## ** data
  > set.seed(10)
  > n <- 5e1
  > dtS <- sampleData(n,outcome="survival")
  > dtS$time <- round(dtS$time,1)
  > dtS$X1 <- factor(rbinom(n, prob = c(0.3,0.4) , size = 2), labels = paste0("T",0:2))
  > 
  > ## ** Compare to fix value / explicit computation 
  > test_that("[ate] G-formula,survival - compare to fix values / explicit computation",{
  +     seqTime <- c(1,5,10)
  +     
  +     fit <- cph(formula = Surv(time,event)~ X1+X2,data=dtS,y=TRUE,x=TRUE)
  +     
  +     ateFit <- ate(fit, data = dtS, treatment = "X1", contrasts = NULL,
  +                   times = seqTime, band = FALSE, verbose = FALSE)
  +     dt.ateFit <- as.data.table(summary(ateFit, band = TRUE, p.value = FALSE, short = 2))
  + 
  +     ## point estimate
  +     expect_equal(dt.ateFit[level == "T0" & type == "meanRisk",estimate],
  +                  c(0.05842145, 0.46346908, 0.66107744),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0" & type == "meanRisk",lower],
  +                  c(0.0000000, 0.2756147, 0.4670461),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0" & type == "meanRisk",upper],
  +                  c(0.1178148, 0.6513235, 0.8551088),
  +                  tol = 1e-6)
  + 
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "diffRisk",estimate],
  +                  c(0.0006541715, 0.0037795576, 0.0040931551),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "diffRisk",lower],
  +                  c(-0.04241346, -0.24364096, -0.26383983),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "diffRisk",upper],
  +                  c(0.0437218, 0.2512001, 0.2720261),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "diffRisk", lowerBand],
  +                  c(-0.04230655, -0.24302677, -0.26317472),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "diffRisk", upperBand],
  +                  c(0.04361489, 0.25058589, 0.27136103),
  +                  tol = 1e-6)
  + 
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "ratioRisk",estimate],
  +                  c(1.011197, 1.008155, 1.006192),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "ratioRisk",lower],
  +                  c(0.2745873, 0.4723104, 0.5997906),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "ratioRisk",upper],
  +                  c(1.747808, 1.543999, 1.412593),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "ratioRisk", lowerBand],
  +                  c(0.26976028, 0.46879896, 0.5971274),
  +                  tol = 1e-6)
  +     expect_equal(dt.ateFit[level == "T0.T1" & type == "ratioRisk", upperBand],
  +                  c(1.75263463, 1.5475109, 1.41525588),
  +                  tol = 1e-6)
  + 
  +     ## standard error
  +     ATE <- list()
  +     ATE.iid <- list()
  +     for(iT in c("T0","T1","T2")){ ## iT <- "T0"
  +         ## risk
  +         newdata0 <- copy(dtS)
  +         newdata0$X1 <- iT
  +         fit <- cph(formula = Surv(time,event)~ X1+X2,data=dtS,y=TRUE,x=TRUE)
  +         resPred <- predictCox(fit, newdata = newdata0, time = seqTime, iid = TRUE)
  +         ATE[[iT]] <- colMeans(1-resPred$survival)
  +         
  +         ATE.iid_term1 <- apply(-resPred$survival.iid,1,rowMeans)
  +         ATE.iid_term2 <- apply(1-resPred$survival, 1, function(x){x-ATE[[iT]]})/n
  +         ATE.iid[[iT]] <- t(ATE.iid_term1) + t(ATE.iid_term2)
  +         ATE.se <- sqrt(apply(ATE.iid[[iT]]^2, 2, sum))
  +         ATE.lower <- ATE[[iT]]+ qnorm(0.025) * ATE.se
  +         ATE.upper <- ATE[[iT]] + qnorm(0.975) * ATE.se
  + 
  +         expect_equal(pmax(0,ATE.lower), dt.ateFit[level == iT,lower])
  +         expect_equal(pmin(1,ATE.upper), dt.ateFit[level == iT,upper])
  +     }
  + 
  +     ## difference in risk
  +     diffATE <- ATE[["T1"]]-ATE[["T0"]]
  +     diffATE.iid <- ATE.iid[["T1"]]-ATE.iid[["T0"]]
  +     diffATE.se <- sqrt(apply(diffATE.iid^2, 2, sum))
  +     diffATE.lower <- diffATE + qnorm(0.025) * diffATE.se
  +     diffATE.upper <- diffATE + qnorm(0.975) * diffATE.se
  + 
  +     expect_equal(diffATE.lower, dt.ateFit[level == "T0.T1" & type == "diffRisk", lower])
  +     expect_equal(diffATE.lower, c(-0.04241346, -0.24364096, -0.26383983), tol = 1e-6)
  + 
  +     expect_equal(diffATE.upper, dt.ateFit[level == "T0.T1" & type == "diffRisk", upper])
  +     expect_equal(diffATE.upper, c(0.0437218, 0.2512001, 0.2720261), tol = 1e-6)
  + 
  +     ## ratio of the risks
  +     ratioATE <- ATE[["T1"]]/ATE[["T0"]]
  +     ratioATE.iid <- rowScale_cpp(ATE.iid[["T1"]],ATE[["T0"]])-rowMultiply_cpp(ATE.iid[["T0"]], ATE[["T1"]]/ATE[["T0"]]^2)
  +     expect_equal(unname(ATE.iid[["T1"]]),unname(ateFit$iid$GFORMULA$T1))
  + 
  +     ratioATE.se <- sqrt(rowSums(ratioATE.iid^2))
  +     ratioATE.se <- sqrt(colSums(ratioATE.iid^2))
  +     ratioATE.lower <- ratioATE + qnorm(0.025) * ratioATE.se
  +     ratioATE.upper <- ratioATE + qnorm(0.975) * ratioATE.se
  +     
  +     expect_equal(ratioATE.lower, dt.ateFit[level == "T0.T1" & type == "ratioRisk", lower])
  +     expect_equal(ratioATE.lower, c(0.2745873, 0.47231036, 0.59979056), tol = 1e-6)
  +     expect_equal(ratioATE.upper, dt.ateFit[level == "T0.T1" & type == "ratioRisk", upper])
  +     expect_equal(ratioATE.upper, c(1.7478076, 1.5439995, 1.41259273))
  + })
  Test passed 🌈
  > 
  > ## ** strata argument
  > test_that("[ate] G-formula,survival - strata argument",{
  +     ## Cox model
  +     e.coxph <- coxph(Surv(time, event) ~ X1, data = dtS,
  +                      x = TRUE, y = TRUE)
  +     outATE <- ate(e.coxph, data = dtS, treatment = NULL, strata = "X1", time = 1,
  +                   product.limit = FALSE, verbose = FALSE, band = TRUE)
  +     dt.outATE <- as.data.table(outATE)
  + 
  +     outPred <- predictCox(e.coxph,
  +                           newdata = dtS,
  +                           times = 1,
  +                           iid = TRUE,
  +                           se = TRUE,
  +                           keep.newdata = TRUE,
  +                           type = "survival")
  + 
  +     GS.se <- as.data.table(outPred)[,.SD[1],by = "X1"]
  +     setkeyv(GS.se, cols = "X1")
  + 
  +     expect_equal(1-GS.se$survival, dt.outATE[type=="meanRisk",estimate])
  +     expect_equal(GS.se$survival.se, dt.outATE[type=="meanRisk",se])
  + 
  +     ## Stratified Cox model
  +     e.coxph <- coxph(Surv(time, event) ~ strata(X1), data = dtS,
  +                      x = TRUE, y = TRUE)
  +     outATE <- ate(e.coxph, data = dtS, treatment = NULL, strata = "X1", time = 3,
  +                   verbose = FALSE, band = TRUE)
  +     dt.outATE <- as.data.table(outATE)
  +     
  +     outPred <- predictCox(e.coxph,
  +                           newdata = dtS,
  +                           times = 3,
  +                           se = TRUE,
  +                           keep.newdata = TRUE,
  +                           type = "survival")
  +     
  +     GS.se <- as.data.table(outPred)[,.SD[1],by = "X1"]
  +     setkeyv(GS.se, cols = "X1")
  + 
  +     expect_equal(1-GS.se$survival, dt.outATE[type == "meanRisk",estimate])
  +     expect_equal(GS.se$survival.se, dt.outATE[type == "meanRisk",se])
  +     expect_equal(GS.se$survival.se, dt.outATE[type == "meanRisk",se])
  + })    
  Test passed 😸
  > 
  > 
  > ## * [ate] Logistic regression
  > cat("[ate] Logistic regression \n")
  [ate] Logistic regression 
  > 
  > n <- 100
  > set.seed(10)
  > dtB <- sampleData(n, outcome="binary")
  > ## dtB <- data.table(X1 = rbinom(n,size =1,prob=0.5),
  >                   ## X2 = rbinom(n,size =1,prob=0.5))
  > ## dtB$Y <- rbinom(n, size = 1, prob = expit(dtB$X1+dtB$X2))
  > ## dtB$X1 <- as.factor(dtB$X1)
  > ## dtB$X2 <- as.factor(dtB$X2)
  > 
  > test_that("[ate] logistic regression - compare to lava",{
  +     fitY <- glm(formula = Y ~ X1 + X2, data=dtB, family = "binomial")
  +     fitT <- glm(formula = X1 ~ X2, data=dtB, family = "binomial")
  + 
  +     ## G-formula
  +     e0.ate <- ate(fitY, data = dtB, treatment = "X1",
  +                  se = TRUE, iid = TRUE, B = 0, verbose = FALSE)
  +     e.ate <- ate(fitY, data = dtB, treatment = "X1",
  +                  times = 5, 
  +                  se = TRUE, iid = TRUE, B = 0, verbose = FALSE)
  +     dt.ate <- as.data.table(e.ate)
  +     
  +     e.lava <- estimate(fitY, function(p, data){
  +         a <- p["(Intercept)"] ; b <- p["X11"] ; c <- p["X21"] ;
  +         R.X11 <- expit(a + b + c * (data[["X2"]]=="1"))
  +         R.X10 <- expit(a + c * (data[["X2"]]=="1"))
  +         list(risk0=R.X10,risk1=R.X11,riskdiff=R.X11-R.X10)},
  +         average=TRUE)
  + 
  +     expect_equal(unname(e.lava$coef), dt.ate[type %in% c("meanRisk","diffRisk"),estimate])
  +     expect_equal(unname(e.lava$vcov[1:2,1:2]), unname(crossprod(do.call(cbind,e.ate$iid$GFORMULA))))
  +     expect_equal(unname(sqrt(diag(e.lava$vcov))), c(dt.ate[type %in% c("meanRisk","diffRisk"),se]))
  + 
  +     ## AIPTW
  +     e.ate2 <- ate(fitY,
  +                   treatment = fitT,
  +                   data = dtB, 
  +                   times = 5, 
  +                   se = TRUE, iid = TRUE, B = 0, verbose = FALSE
  +                   )
  +     dt.ate2 <- as.data.table(e.ate2)
  +     
  +     dtB$Y0 <- 0
  +     dtB$Y1 <- 1
  + 
  +     ## iid ate
  +     iPredT <- predict(fitT, type = "response")
  +     dtBC <- rbind(cbind(dtB[,.(Y,X2)], X1 = factor(0, levels = levels(dtB$X1)), X1test = dtB$X1=="0", pi = 1-iPredT),
  +                   cbind(dtB[,.(Y,X2)], X1 = factor(1, levels = levels(dtB$X1)), X1test = dtB$X1=="1", pi = iPredT))
  +     dtBC$r <- predict(fitY, newdata = dtBC, type = "response")
  +     dtBC[, ate := Y*X1test/pi + r*(1-X1test/pi)]
  +     dtBC[, ate.iid := (ate-mean(ate))/.N, by = "X1"]
  + 
  +     ## iid outcome
  +     iid.risk <- attr(predictRisk(fitY, newdata = dtBC, iid = TRUE),"iid")
  +     nuisanceY.iid <- rowMultiply_cpp(iid.risk, scale = (1-dtBC$X1test/dtBC$pi))
  +     dtBC$AnuisanceY.iid <- c(rowMeans(nuisanceY.iid[,dtBC$X1=="0"]),rowMeans(nuisanceY.iid[,dtBC$X1=="1"]))
  + 
  +     ## iid treatment
  +     iid.pi <- attr(predictRisk(fitT, newdata = dtBC, iid = TRUE),"iid")
  +     nuisanceT.iid <- rowMultiply_cpp(iid.pi, scale = (-1)^(dtBC$X1=="1")*dtBC$X1test*(dtBC$Y-dtBC$r)/dtBC$pi^2)
  +     dtBC$AnuisanceT.iid <- c(rowMeans(nuisanceT.iid[,dtBC$X1=="0"]),rowMeans(nuisanceT.iid[,dtBC$X1=="1"]))
  +     
  +     ## global
  +     expect_equal(dtBC[, mean(ate), by = "X1"][[2]],
  +                  dt.ate2[type=="meanRisk",estimate])
  +     expect_equal(dtBC[, sqrt(sum( (ate.iid+AnuisanceY.iid+AnuisanceT.iid)^2 )), by = "X1"][[2]],
  +                  dt.ate2[type=="meanRisk",se])
  +     ## butils:::object2script(dt.ate2)
  +     ## GS <- data.table("type" = c("ate", "ate", "diffAte", "ratioAte"), 
  +     ##                  "level" = c("0", "1", "0.1", "0.1"), 
  +     ##                  "time" = c(5, 5, 5, 5), 
  +     ##                  "value" = c(0.4167808, 0.7341667, 0.3173858, 1.7615174), 
  +     ##                  "se" = c(0.05227194, 0.13439533, 0.14429419, 0.39122437), 
  +     ##                  "lower" = c(0.31432969, 0.47075666, 0.03457442, 0.99473171), 
  +     ##                  "upper" = c(0.5192320, 0.9975767, 0.6001973, 2.5283031), 
  +     ##                  "estimator" = c("AIPTW", "AIPTW", "AIPTW", "AIPTW"))
  +     
  + })
  ── Error (???): [ate] logistic regression - compare to lava ────────────────────
  Error in `estimate(fitY, function(p, data) {
      a <- p["(Intercept)"]
      b <- p["X11"]
      c <- p["X21"]
      R.X11 <- expit(a + b + c * (data[["X2"]] == "1"))
      R.X10 <- expit(a + c * (data[["X2"]] == "1"))
      list(risk0 = R.X10, risk1 = R.X11, riskdiff = R.X11 - R.X10)
  }, average = TRUE)`: could not find function "estimate"
  
  Error: Test failed
  Execution halted
Running the tests in ‘tests/test-phreg.R’ failed.
Complete output:
  > ### test-phreg.R --- 
  > #----------------------------------------------------------------------
  > ## author: Brice Ozenne
  > ## created: feb 28 2017 (09:52) 
  > ## Version: 
  > ## last-updated: Dec  7 2020 (09:23) 
  > ##           By: Thomas Alexander Gerds
  > ##     Update #: 22
  > #----------------------------------------------------------------------
  > ## 
  > ### Commentary: 
  > ## 
  > ### Change Log:
  > #----------------------------------------------------------------------
  > ## 
  > ### Code:
  > 
  > ## * Settings
  > library(testthat)
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(rms)
  Loading required package: Hmisc
  Loading required package: lattice
  Loading required package: survival
  Loading required package: Formula
  Loading required package: ggplot2
  
  Attaching package: 'Hmisc'
  
  The following object is masked from 'package:testthat':
  
      describe
  
  The following objects are masked from 'package:base':
  
      format.pval, units
  
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > library(survival)
  > library(prodlim)
  > library(mets)
  Loading required package: timereg
  Loading required package: lava
  
  Attaching package: 'lava'
  
  The following object is masked from 'package:ggplot2':
  
      vars
  
  The following object is masked from 'package:testthat':
  
      compare
  
  mets version 1.2.9
  
  Attaching package: 'mets'
  
  The following object is masked from 'package:rms':
  
      bplot
  
  > 
  > context("Compatibility with the phreg function from the mets package")
  > 
  > ## * Survival
  > 
  > ## ** Data
  > set.seed(10)
  > d <- sampleData(5e1, outcome = "survival")[,.(eventtime,event,X1,X2,X3,X4,X6)]
  > d[ , Xcat2 := as.factor(paste0(X1,X2))]
  > d[, entry:= 0]
  > d[, id := 1:.N]
  > d[, group := factor(rbinom(.N,2,0.5))]
  > 
  > ## ** Tests
  > test_that("survival - no strata", {
  +     m.phreg <- phreg(Surv(entry,eventtime,event)~X1+X2,data=d)
  +     m.coxph <- coxph(Surv(eventtime,event)~X1+X2,data=d, x = TRUE, y = TRUE)
  +     expect_equal(predictCox(m.phreg),predictCox(m.coxph), tol = 1e-8)
  + 
  +     pred.phreg <- predictCox(m.phreg, newdata = d, time = 1:5, se = TRUE)
  +     pred.coxph <- predictCox(m.coxph, newdata = d, time = 1:5, se = TRUE)
  +     expect_equal(pred.phreg$survival, pred.coxph$survival, tol = 1e-8)
  + })
  ── Failure (???): survival - no strata ─────────────────────────────────────────
  predictCox(m.phreg) not equal to predictCox(m.coxph).
  Component "cumhazard": Mean relative difference: 0.3224839
  Component "survival": Mean relative difference: 0.1827965
  
  Error: Test failed
  Execution halted
Running the tests in ‘tests/test-predictCSC.R’ failed.
Complete output:
  > ### test-predictCSC.R --- 
  > #----------------------------------------------------------------------
  > ## author: Brice Ozenne
  > ## created: maj 18 2017 (09:23) 
  > ## Version: 
  > ## last-updated: Dec  7 2020 (09:18) 
  > ##           By: Thomas Alexander Gerds
  > ##     Update #: 303
  > #----------------------------------------------------------------------
  > ## 
  > ### Commentary: 
  > ##
  > ### Change Log:
  > #----------------------------------------------------------------------
  > ## 
  > ### Code:
  > 
  > ## * Settings
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(testthat)
  > library(data.table)
  > library(prodlim)
  > library(rms)
  Loading required package: Hmisc
  Loading required package: lattice
  Loading required package: survival
  Loading required package: Formula
  Loading required package: ggplot2
  
  Attaching package: 'Hmisc'
  
  The following object is masked from 'package:testthat':
  
      describe
  
  The following objects are masked from 'package:base':
  
      format.pval, units
  
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > library(survival)
  > library(timereg); nsim.band <- 100
  > 
  > ## * [predictCSC] predict.CSC vs. manual calculation and mstate on small examples
  > cat("[predictCSC] predict.CSC vs. manual calculation on small examples \n")
  [predictCSC] predict.CSC vs. manual calculation on small examples 
  > ## ** Data
  > ## simulatenous events
  > df1 <- data.frame(time = rep(1:10,2),
  +                  event = c(rep(1,10),rep(2,10)),
  +                  X1 = 0:1
  +                  )
  > df1$event1 <- as.numeric(df1$event == 1)
  > df1$event2 <- as.numeric(df1$event == 2)
  > 
  > set.seed(11)
  > dfS <- rbind(cbind(df1, grp = 1, X2 = rbinom(20,1,.4)),
  +              cbind(rbind(df1,df1), grp = 2, X2 = rbinom(20,1,.4)),
  +              cbind(df1, grp = 3, X2 = rbinom(20,1,.4))
  +              )
  > 
  > ## distinct events
  > df2 <- data.frame(time = c(1:10,-0.1+1:10),
  +                   event = c(rep(1,10),rep(2,10)),
  +                   X1 = 0:1
  +                   )
  > df2$event1 <- as.numeric(df2$event == 1)
  > df2$event2 <- as.numeric(df2$event == 2)
  > 
  > ## ** No risk factor no strata
  > test_that("[predictCSC] (no strata, data=df1): compare to manual estimation",{
  + 
  +     CSC.exp <- CSC(Hist(time,event) ~ 1, data = df1)
  +     pred.exp <- predict(CSC.exp, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = FALSE)  
  +     pred.exp2 <- predict(CSC.exp, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = TRUE) 
  +     CSC.prodlim <- CSC(Hist(time,event) ~ 1, data = df1, surv.type = "survival", method = "breslow")
  +     pred.prodlim <- predict(CSC.prodlim, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = TRUE)  
  +     expect_equal(pred.prodlim,pred.exp2)
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ 1, data = df1, surv.type = "survival", method = "efron")
  +     pred.prodlimE <- predict(CSC.prodlimE, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = TRUE)  
  + 
  +     ## test baseline hazard
  +     lambda1 <- lambda2 <- 1/seq(20,2,by = -2)
  +     lambda2E <- predictCox(CSC.prodlimE$models[[2]], type = "hazard")$hazard
  +     expect_equal(predictCox(CSC.exp$models[[1]], type = "hazard")$hazard,
  +                  lambda1)
  +     expect_equal(predictCox(CSC.exp$models[[2]], type = "hazard")$hazard,
  +                  lambda2)
  +     expect_equal(predictCox(CSC.prodlim$models[[1]], type = "hazard")$hazard,
  +                  lambda1)
  +     expect_equal(predictCox(CSC.prodlim$models[[2]], type = "hazard")$hazard,
  +                  lambda1+lambda2)
  +     expect_equal(basehaz(CSC.prodlimE$models[[2]])$hazard,
  +                  cumsum(lambda2E)
  +                  )
  + 
  +     ## test absolute risk
  +     survival <- c(1,exp(cumsum(-lambda1-lambda2)))[1:10]
  +     expect_equal(as.double(pred.exp$absRisk),
  +                  cumsum(lambda1*survival))
  +     survival <- c(1,cumprod(1-(lambda1+lambda2)))[1:10]
  +     expect_equal(as.double(pred.prodlim$absRisk),
  +                  cumsum(lambda1*survival))
  +     survival <- c(1,cumprod(1-lambda2E))[1:10]
  +     expect_equal(as.double(pred.prodlimE$absRisk),
  +                  cumsum(lambda1*survival))
  + })
  Test passed 🎉
  > 
  > test_that("[predictCSC] (no strata, data=df1): compare to manual estimation",{
  +     CSC.exp <- CSC(Hist(time,event) ~ 1, data = df2)
  +     pred.exp <- predict(CSC.exp, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = FALSE)  
  +     pred.exp2 <- predict(CSC.exp, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = TRUE)  
  +     CSC.prodlim <- CSC(Hist(time,event) ~ 1, data = df2, surv.type = "survival", method = "breslow")
  +     pred.prodlim <- predict(CSC.prodlim, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = TRUE)  
  +     expect_equal(pred.exp2, pred.prodlim)
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ 1, data = df2, surv.type = "survival", method = "efron")
  +     pred.prodlimE <- predict(CSC.prodlimE, times = 1:10, newdata = data.frame(NA), cause = 1, product.limit = TRUE)  
  +     ## test baseline hazard
  +     lambda1 <- 1/seq(19,1,by = -2)
  +     lambda2 <- 1/seq(20,2,by = -2)
  +     lambda2E <- as.data.table(predictCox(CSC.prodlimE$models[[2]], type = "hazard"))
  +     expect_equal(predictCox(CSC.exp$models[[1]], times = 1:10, type = "hazard")$hazard,
  +                  lambda1)
  +     expect_equal(predictCox(CSC.exp$models[[2]], times = -0.1 + 1:10, type = "hazard")$hazard,
  +                  lambda2)
  +     expect_equal(predictCox(CSC.prodlim$models[[1]], times = 1:10, type = "hazard")$hazard,
  +                  lambda1)
  +     expect_equal(predictCox(CSC.prodlim$models[[2]], type = "hazard")$hazard,
  +                  as.double(rbind(lambda2,lambda1)))
  +     expect_equal(basehaz(CSC.prodlimE$models[[2]])$hazard,
  +                  cumsum(lambda2E$hazard)
  +                  )
  +     ## test absolute risk
  +     vec.lambda <- cumsum(as.double(rbind(lambda2,lambda1)))
  +     survival <- exp(-matrix(vec.lambda, ncol = 2, byrow = TRUE)[,1])
  +     expect_equal(as.double(pred.exp$absRisk),
  +                  cumsum(lambda1*survival))
  +     vec.lambda <- as.double(rbind(lambda2,lambda1))
  +     survival <- matrix(cumprod(1-vec.lambda), ncol = 2, byrow = TRUE)[,1]
  +     expect_equal(as.double(pred.prodlim$absRisk),
  +                  cumsum(lambda1*survival))
  +     vec.lambda <-lambda2E$hazard
  +     survival <- matrix(cumprod(1-vec.lambda), ncol = 2, byrow = TRUE)[,1]
  +     expect_equal(as.double(pred.prodlimE$absRisk),
  +                  cumsum(lambda1*survival))
  + })
  Test passed 😸
  > 
  > test_that("[predictCSC]: compare to mstate",{
  +     if (!requireNamespace("mstate",quietly=TRUE)){
  +         message("Package mstate not installed. Skip this test of predictCSC.")
  +     }else{
  +         library(mstate)
  +         for(iData in 1:2){# iData <- 1
  +             df <- switch(as.character(iData),
  +                          "1"=df1,
  +                          "2"=df2)
  +             ## fit CSC
  +             CSC.exp <- CSC(Hist(time,event) ~ 1, data = df)
  +             ## fit mstate
  +             tmat <- mstate::trans.comprisk(2, names = c("0", "1", "2"))
  +             dL <- mstate::msprep(time = c(NA, "time", "time"),
  +                          status = c(NA,"event1", "event2"),
  +                          data = df, keep = c("X1"),
  +                          trans = tmat)
  +             dL.exp <- expand.covs(dL,  c("X1"))
  +             e.coxph <- coxph(Surv(time, status) ~ strata(trans),
  +                              data = dL.exp)
  +             ## compare
  +             newdata.L <- data.frame(trans = c(1, 2), strata = c(1, 2))
  +             newdata <- data.frame(NA)
  +             pred.msfit <- msfit(e.coxph, newdata = newdata.L, trans = tmat)
  +             suppressWarnings(
  +                 pred.probtrans <- probtrans(pred.msfit,0)[[1]]
  +             )
  +             pred.RR <- predict(CSC.exp,
  +                                times = pred.probtrans[,"time"],
  +                                newdata = newdata,
  +                                cause = 1,
  +                                se = FALSE,
  +                                product.limit = TRUE)
  +             expect_equal(pred.probtrans[,"pstate2"],
  +                          as.double(pred.RR$absRisk)
  +                          )
  +         }
  +     }
  + })
  Test passed 😀
  > 
  > 
  > test_that("[predict.CSC] check absolute risks add up to one",{
  +     for(iData in 1:2){
  +         df <- switch(as.character(iData),
  +                      "1"=df1,
  +                      "2"=df2)
  +         ## fit CSC
  +         CSC.exp <- CSC(Hist(time,event) ~ 1, data = df, method = "breslow")
  +         ## compute all transition probabilites
  +         seqTimes <- sort(unique(df$time))
  +         nTimes <- length(seqTimes)
  +         newdata <- data.frame(NA)
  +         hazAll <- predictCox(CSC.exp$models[[1]], times = seqTimes, newdata = newdata, type = "hazard")$hazard+predictCox(CSC.exp$models[[2]], times = seqTimes, newdata = newdata, type = "hazard")$hazard
  +         surv.prodlim <- cumprod(1-as.double(hazAll))
  +         haz1.prodlim <- predict(CSC.exp, times = seqTimes, newdata = newdata, cause = 1, product.limit = TRUE)  
  +         haz2.prodlim <- predict(CSC.exp, times = seqTimes, newdata = newdata, cause = 2, product.limit = TRUE)  
  +         expect_equal(surv.prodlim+as.double(haz1.prodlim$absRisk)+as.double(haz2.prodlim$absRisk),rep(1,nTimes))
  +     }
  + })
  Test passed 🎉
  > 
  > ## ** No risk factor strata
  > test_that("[predict.CSC] (strata): compare to manual estimation",{
  +     ## as.data.table(dfS)[,.N, by = "grp"]
  +     CSC.exp <- CSC(Hist(time,event) ~ strata(grp), data = dfS, method = "breslow")    
  +     pred.exp <- predict(CSC.exp, times = 1:10, newdata = data.frame(grp = 1:3), cause = 1, product.limit = FALSE)
  +     pred.exp2 <- predict(CSC.exp, times = 1:10, newdata = data.frame(grp = 1:3), cause = 1, product.limit = TRUE)  
  + 
  +     CSC.prodlim <- CSC(Hist(time,event) ~ strata(grp), data = dfS, surv.type = "survival", method = "breslow")
  +     pred.prodlim <- predict(CSC.prodlim, times = 1:10, newdata = data.frame(grp = 1:3), cause = 1, product.limit = TRUE)  
  +     expect_equal(pred.exp2,pred.prodlim)
  +     
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ strata(grp), data = dfS, surv.type = "survival", method = "efron")
  +     pred.prodlimE <- predict(CSC.prodlimE, times = 1:10, newdata = data.frame(grp = 1:3), cause = 1, product.limit = TRUE)  
  + 
  +     ## baseline
  +     baseline1 <- as.data.table(predictCox(CSC.exp$models[[1]], type = "hazard"))
  +     lambda1 <- baseline1[strata=="grp=1",hazard]
  +     expect_equal(lambda1, 1/seq(20,2,by = -2))
  + 
  +     baseline2 <- as.data.table(predictCox(CSC.exp$models[[1]], type = "hazard"))
  +     lambda2 <- baseline2[strata=="grp=1",hazard]
  +     expect_equal(lambda2, 1/seq(20,2,by = -2))
  +         
  +     baseline1E <- as.data.table(predictCox(CSC.prodlimE$models[[1]], type = "hazard"))
  +     lambda1E.1 <- baseline1E[strata=="grp=1",hazard]
  +     lambda1E.2 <- baseline1E[strata=="grp=2",hazard]
  +     expect_equal(lambda1E.1, 1/seq(20,2,by = -2))
  + 
  +     baseline2E <- as.data.table(predictCox(CSC.prodlimE$models[[2]], type = "hazard"))
  +     lambda2E.1 <- baseline2E[strata=="grp=1",hazard]
  +     lambda2E.2 <- baseline2E[strata=="grp=2",hazard]
  +     ## 
  +         
  +     ## test absolute risk
  +     survival <- c(1,exp(cumsum(-lambda1-lambda2)))[1:10]
  +     expect_equal(as.double(pred.exp$absRisk[1,]),
  +                  cumsum(lambda1*survival))
  +     expect_equal(as.double(pred.exp$absRisk[1,]),
  +                  as.double(pred.exp$absRisk[2,]))
  +     expect_equal(as.double(pred.exp$absRisk[1,]),
  +                  as.double(pred.exp$absRisk[3,]))
  + 
  +     survival <- c(1,cumprod(1-(lambda1+lambda2)))[1:10]
  +     expect_equal(as.double(pred.prodlim$absRisk[1,]),
  +                  cumsum(lambda1*survival))
  +     expect_equal(as.double(pred.prodlim$absRisk[1,]),
  +                  as.double(pred.prodlim$absRisk[2,]))
  +     expect_equal(as.double(pred.prodlim$absRisk[1,]),
  +                  as.double(pred.prodlim$absRisk[3,]))
  + 
  +     survival <- c(1,cumprod(1-lambda2E.1))[1:10]
  +     expect_equal(as.double(pred.prodlimE$absRisk[pred.prodlimE$strata=="grp=1"]),
  +                  cumsum(lambda1E.1*survival))
  +     survival <- c(1,cumprod(1-lambda2E.2))[1:10]
  +     expect_equal(as.double(pred.prodlimE$absRisk[pred.prodlimE$strata=="grp=2"]),
  +                  cumsum(lambda1E.2*survival))    
  +     expect_equal(as.double(pred.prodlimE$absRisk[1,]),
  +                  as.double(pred.prodlimE$absRisk[3,]))
  + })
  Test passed 🌈
  > 
  > ## ** risk factor no strata
  > test_that("[predict.CSC] (risk factor, strata): compare to manual estimation",{
  +     CSC.exp <- CSC(Hist(time,event) ~ X1, data = df1)
  +     CSC.prodlim <- CSC(Hist(time,event) ~ X1, data = df1, surv.type = "survival", method = "breslow")
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ X1, data = df1, surv.type = "survival", method = "efron")
  + 
  +     calcCIF <- function(CSC, X){
  +         lambda01 <- as.data.table(predictCox(CSC$models[[1]], centered = FALSE, type = "hazard"))
  +         lambda02 <- as.data.table(predictCox(CSC$models[[2]], centered = FALSE, type = "hazard"))
  +         eXb.1 <- as.double(exp(X%*%coef(CSC)[[1]]))
  +         eXb.2 <- as.double(exp(X%*%coef(CSC)[[2]]))
  +         if(CSC$surv.type == "survival"){
  +             survival <- cumprod(1-eXb.2*lambda02$hazard)
  +         }else{           
  +             Lambda01 <- cumsum(lambda01$hazard)
  +             Lambda02 <- cumsum(lambda02$hazard)
  +             survival <- exp(-eXb.1*Lambda01-eXb.2*Lambda02)
  +         }
  +         survival <- c(1,survival[-length(survival)])
  +         absRisk <- cumsum(eXb.1*lambda01$hazard*survival)
  +         return(absRisk)
  +     }
  +     for(iX1 in 0:1){
  +         newdata <- data.frame(X1=iX1)
  +         pred.exp <- predict(CSC.exp, times = 1:10, newdata = newdata, cause = 1, product.limit = FALSE)  
  +         expect_equal(as.double(pred.exp$absRisk),
  +                      calcCIF(CSC.exp, X = as.matrix(newdata))
  +                      )
  +         pred.prodlim <- predict(CSC.prodlim, times = 1:10, newdata = newdata, cause = 1, product.limit = TRUE)
  +         expect_equal(as.double(pred.prodlim$absRisk),
  +                      calcCIF(CSC.prodlim, X = as.matrix(newdata))
  +                      )
  +         pred.prodlimE <- predict(CSC.prodlimE, times = 1:10, newdata = newdata, cause = 1, product.limit = TRUE)  
  +         expect_equal(as.double(pred.prodlimE$absRisk),
  +                      calcCIF(CSC.prodlimE, X = as.matrix(newdata))
  +                      )
  +     }
  + })
  Test passed 🎉
  > 
  > test_that("[predict.CSC] - compare to manual estimation",{
  +     CSC.exp <- CSC(Hist(time,event) ~ X1, data = df2, method = "breslow")
  + 
  +     CSC.prodlim <- CSC(Hist(time,event) ~ X1, data = df2, surv.type = "survival", method = "breslow")
  +     
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ X1, data = df2, surv.type = "survival", method = "efron")
  + 
  +     calcCIF <- function(CSC, X){
  + 
  +         lambda01 <- as.data.table(predictCox(CSC$models[[1]], centered = FALSE, type = "hazard"))
  +         lambda02 <- as.data.table(predictCox(CSC$models[[2]], centered = FALSE, type = "hazard"))
  +         eXb.1 <- as.double(exp(X%*%coef(CSC)[[1]]))
  +         eXb.2 <- as.double(exp(X%*%coef(CSC)[[2]]))
  + 
  +         if(CSC$surv.type == "survival"){
  +             survival <- cumprod(1-eXb.2*lambda02$hazard)
  +         }else{           
  +             Lambda01 <- cumsum(lambda01$hazard)
  +             Lambda02 <- cumsum(lambda02$hazard)
  +             survival <- exp(-eXb.1*Lambda01-eXb.2*Lambda02)
  +         }
  + 
  +         survival <- c(1,survival[-length(survival)])
  +         absRisk <- cumsum(eXb.1*lambda01$hazard*survival)
  +         return(absRisk)
  +         
  +     }
  +  
  +     for(iX1 in 0:1){ # iX1 <- 0
  +         newdata <- data.frame(X1=iX1)
  +         pred.exp <- predict(CSC.exp, times = sort(unique(df2$time)), newdata = newdata, cause = 1, product.limit = FALSE)  
  +         expect_equal(as.double(pred.exp$absRisk),
  +                      calcCIF(CSC.exp, X = as.matrix(newdata))
  +                      )
  +     
  +         pred.prodlim <- predict(CSC.prodlim, times = sort(unique(df2$time)), newdata = newdata, cause = 1, product.limit = TRUE)  
  +         expect_equal(as.double(pred.prodlim$absRisk),
  +                      calcCIF(CSC.prodlim, X = as.matrix(newdata))
  +                      )
  + 
  +         pred.prodlimE <- predict(CSC.prodlimE, times = sort(unique(df2$time)), newdata = newdata, cause = 1, product.limit = TRUE)  
  +         expect_equal(as.double(pred.prodlimE$absRisk),
  +                      calcCIF(CSC.prodlimE, X = as.matrix(newdata))
  +                      )
  +     }
  + })
  Test passed 🥇
  > 
  > 
  > test_that("[predict.CSC]: compare to mstate",{
  +     if (!requireNamespace("mstate",quietly=TRUE)){
  +         message("Package mstate not installed. Skip this test of predictCSC.")
  +     }else{
  +         library(mstate)
  +         for(iData in 1:2){
  +             df <- switch(as.character(iData),
  +                          "1"=df1,
  +                          "2"=df2)
  +             ## fit CSC
  +             CSC.exp <- CSC(Hist(time,event) ~ X1, data = df, method = "breslow")
  +             ## fit mstate
  +             df$event1 <- as.numeric(df$event == 1)
  +             df$event2 <- as.numeric(df$event == 2)
  +             tmat <- mstate::trans.comprisk(2, names = c("0", "1", "2"))
  +             dL <- mstate::msprep(time = c(NA, "time", "time"),
  +                          status = c(NA,"event1", "event2"),
  +                          data = df, keep = c("X1"),
  +                          trans = tmat)
  +             dL.exp <- expand.covs(dL,  c("X1"))
  +             e.coxph <- coxph(Surv(time, status) ~ X1.1 + X1.2 + strata(trans),
  +                              data = dL.exp)
  +             ## compare
  +             for(iX1 in 0:1){ # iX1 <- 0
  +                 newdata.L <- data.frame(X1.1 = c(iX1,0), X1.2 = c(0,iX1), trans = c(1, 2), strata = c(1, 2))    
  +                 newdata <- data.frame(X1 = iX1)
  +                 pred.msfit <- msfit(e.coxph, newdata = newdata.L, trans = tmat)
  +                 suppressWarnings(
  +                     pred.probtrans <- probtrans(pred.msfit,0)[[1]]
  +                 )
  +                 pred.RR <- predict(CSC.exp, times = pred.probtrans[,"time"], newdata = newdata, cause = 1, product.limit = TRUE)
  +                 expect_equal(pred.probtrans[,"pstate2"],
  +                              as.double(pred.RR$absRisk)
  +                              )
  +             }
  +         }
  +     }
  + })
  Test passed 🥳
  > 
  > test_that("[predict.CSC]: check absolute risks add up to one",{
  +     for(iData in 1:2){
  +         df <- switch(as.character(iData),
  +                      "1"=df1,
  +                      "2"=df2)
  + 
  +         ## fit CSC
  +         CSC.exp <- CSC(Hist(time,event) ~ X1, data = df, method = "breslow")
  + 
  +         ## compute all transition probabilites
  +         seqTimes <- sort(unique(df$time))
  +         nTimes <- length(seqTimes)
  +   
  +         for(iX1 in 0:1){            
  +             newdata <- data.frame(X1 = iX1)
  + 
  +             hazAll <- predictCox(CSC.exp$models[[1]], times = seqTimes, newdata = newdata, type = "hazard")$hazard+predictCox(CSC.exp$models[[2]], times = seqTimes, newdata = newdata, type = "hazard")$hazard
  +             surv.prodlim <- cumprod(1-as.double(hazAll))
  +             haz1.prodlim <- predict(CSC.exp, times = seqTimes, newdata = newdata, cause = 1, product.limit = TRUE)  
  +             haz2.prodlim <- predict(CSC.exp, times = seqTimes, newdata = newdata, cause = 2, product.limit = TRUE)  
  + 
  +             expect_equal(surv.prodlim+as.double(haz1.prodlim$absRisk)+as.double(haz2.prodlim$absRisk),rep(1,nTimes))
  +         }
  +     }
  + })
  Test passed 🎊
  > 
  > test_that("[predict.CSC]: mstate simulated data",{
  +     if (!requireNamespace("mstate",quietly=TRUE)){
  +         message("Package mstate not installed. Skip this test of predictCSC.")
  +     }else{
  +         library(mstate)
  +         ## * [predictCSC] predict.CSC vs. mstate on simulated and real data
  +         cat("[predictCSC] predict.CSC vs. mstate \n")
  +         ## ** Data
  +         set.seed(10)
  +         d <- sampleData(1e2, outcome = "competing.risks")[,.(time,event,X1,X2,X6)]
  +         d[,X1:=as.numeric(as.character(X1))]
  +         d[,X2:=as.numeric(as.character(X2))]
  +         d[ , X16 := X1*X6]
  +         d[ , Xcat2 := as.factor(paste0(X1,X2))]
  +         d <- d[event!=0]
  +         d[, event1 := as.numeric(event == 1)]
  +         d[, event2 := as.numeric(event == 2)]
  +         ## ** Mstate
  +         tmat <- mstate::trans.comprisk(2, names = c("0", "1", "2"))
  +         dL <- mstate::msprep(time = c(NA, "time", "time"),
  +                              status = c(NA,"event1", "event2"),
  +                              data = d, keep = c("X1","X2","X16","Xcat2"),
  +                              trans = tmat)
  +         dL.exp <- expand.covs(dL,  c("X1","X2","X16","Xcat2"))
  +         ## ** No covariates
  +         ## test_that("predict.CSC (no covariates): compare to mstate",{
  +         newdata <- data.frame(NA)
  +         newdata.L <- data.frame(trans = c(1, 2), strata = c(1, 2))
  +         # mstate
  +         e.coxph <- coxph(Surv(time, status) ~1 + strata(trans),
  +                          data = dL.exp)
  +         pred.msfit <- msfit(e.coxph, newdata = newdata.L, trans = tmat)
  +         pred.probtrans <- probtrans(pred.msfit,0)[[1]]
  +         ## riskRegression
  +         CSC.RR1 <- CSC(Hist(time,event)~1, data = d, method = "breslow")
  +         pred.RR1a <- predict(CSC.RR1, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR1b <- predict(CSC.RR1, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         pred.RR1c <- predict(CSC.RR1, newdata, cause = 2, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR1d <- predict(CSC.RR1, newdata, cause = 2, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         CSC.RR2 <- CSC(Hist(time,event)~1, data = d, surv.type = "survival", method = "breslow")
  +         pred.RR2a <- predict(CSC.RR2, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR2b <- predict(CSC.RR2, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         expect_equal(as.double(pred.RR1a$absRisk),pred.probtrans[,"pstate2"])
  +         expect_equal(as.double(pred.RR1b$absRisk),pred.probtrans[,"pstate2"], tol = 5e-3)
  +         expect_equal(as.double(pred.RR1c$absRisk),pred.probtrans[,"pstate3"])
  +         expect_equal(as.double(pred.RR1d$absRisk),pred.probtrans[,"pstate3"], tol = 5e-2)
  +         expect_equal(as.double(pred.RR2a$absRisk),pred.probtrans[,"pstate2"])
  +         expect_equal(as.double(pred.RR2b$absRisk),pred.probtrans[,"pstate2"], tol = 5e-3)
  +         # quantile(as.double(pred.RR1a$absRisk.se) - pred.probtrans[,"se2"])
  +         expect_equal(as.double(pred.RR1a$absRisk.se),pred.probtrans[,"se2"], tol = 5e-3)
  +         expect_equal(as.double(pred.RR1b$absRisk.se),pred.probtrans[,"se2"], tol = 5e-3)
  +         expect_equal(as.double(pred.RR1c$absRisk.se),pred.probtrans[,"se3"], tol = 5e-3)
  +         expect_equal(as.double(pred.RR1d$absRisk.se),pred.probtrans[,"se3"], tol = 5e-3)
  +         expect_equal(as.double(pred.RR2a$absRisk.se),pred.probtrans[,"se2"], tol = 5e-3)    
  +         expect_equal(as.double(pred.RR2b$absRisk.se),pred.probtrans[,"se2"], tol = 5e-3)
  +         ## ** With covariates
  +         # test_that("predict.CSC (covariates): compare to mstate",{
  +         for(iX in 0:1){
  +             newdata <- data.frame(X1 = iX, X2 = 0, X16 = 0)
  +             newdata.L <- data.frame(X1.1 = c(iX, 0), X1.2 = c(0, iX),
  +                                     X2.1 = c(0, 0), X2.2 = c(0, 0),
  +                                     X16.1 = c(0, 0), X16.2 = c(0, 0),
  +                                     trans = c(1, 2), strata = c(1, 2))
  +             # mstate
  +             e.coxph <- coxph(Surv(time, status) ~ X1.1 + X1.2 + X2.1 + X2.2 + X16.1 + X16.2 + strata(trans),
  +                              data = dL.exp)
  +             pred.msfit <- msfit(e.coxph, newdata = newdata.L, trans = tmat)
  +             suppressWarnings(
  +                 pred.probtrans <- probtrans(pred.msfit,0)[[1]]
  +             )
  +             ## riskRegression
  +             CSC.RR1 <- CSC(Hist(time,event)~X1+X2+X16, data = d, method = "breslow")
  +             pred.RR1a <- predict(CSC.RR1, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                                  keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +             pred.RR1b <- predict(CSC.RR1, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                                  keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +             pred.RR1c <- predict(CSC.RR1, newdata, cause = 2, time = pred.probtrans[,"time"],
  +                                  keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +             pred.RR1d <- predict(CSC.RR1, newdata, cause = 2, time = pred.probtrans[,"time"],
  +                                  keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +             CSC.RR2 <- CSC(Hist(time,event)~X1+X2+X16, data = d, surv.type = "survival", method = "breslow")
  +             pred.RR2a <- predict(CSC.RR2, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                                  keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +             pred.RR2b <- predict(CSC.RR2, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                                  keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +             expect_equal(as.double(pred.RR1a$absRisk),pred.probtrans[,"pstate2"])
  +             expect_equal(as.double(pred.RR1b$absRisk),pred.probtrans[,"pstate2"], tol = 5e-3)
  +             expect_equal(as.double(pred.RR1c$absRisk),pred.probtrans[,"pstate3"])
  +             expect_equal(as.double(pred.RR1d$absRisk),pred.probtrans[,"pstate3"], tol = 1e-2)
  +             expect_equal(as.double(pred.RR2a$absRisk),pred.probtrans[,"pstate2"], tol = 1e-2)
  +             expect_equal(as.double(pred.RR2b$absRisk),pred.probtrans[,"pstate2"], tol = 1e-1)
  +             #if(iX==0){
  +             # quantile(as.double(pred.RR1a$absRisk.se) - pred.probtrans[,"se2"])
  +             # quantile(as.double(pred.RR1c$absRisk.se) - pred.probtrans[,"se3"])
  +             expect_equal(as.double(pred.RR1a$absRisk.se),pred.probtrans[,"se2"], tol = 1e-2)
  +             expect_equal(as.double(pred.RR1b$absRisk.se),pred.probtrans[,"se2"], tol = 1e-2)
  +             if(iX==0){
  +                 expect_equal(as.double(pred.RR1c$absRisk.se),pred.probtrans[,"se3"], tol = 1e-2)
  +                 expect_equal(as.double(pred.RR1d$absRisk.se),pred.probtrans[,"se3"], tol = 1e-2)
  +             }
  +             expect_equal(as.double(pred.RR2a$absRisk.se),pred.probtrans[,"se2"], tol = 1e-2)    
  +             expect_equal(as.double(pred.RR2b$absRisk.se),pred.probtrans[,"se2"], tol = 1e-2)
  +             ## }else{
  +             ##     expect_equal(as.double(pred.RR1a$absRisk.se),pred.probtrans[,"se2"], tol = 1e-1)
  +             ##     expect_equal(as.double(pred.RR1b$absRisk.se),pred.probtrans[,"se2"], tol = 1e-1)
  +             ##     expect_equal(as.double(pred.RR1c$absRisk.se),pred.probtrans[,"se3"], tol = 1e-1)
  +             ##     expect_equal(as.double(pred.RR1d$absRisk.se),pred.probtrans[,"se3"], tol = 1e-1)
  +             ##     expect_equal(as.double(pred.RR2a$absRisk.se),pred.probtrans[,"se2"], tol = 1e-1)    
  +             ##     expect_equal(as.double(pred.RR2b$absRisk.se),pred.probtrans[,"se2"], tol = 1e-1)
  +             ## }
  +         }
  +         ## ** Strata 
  +         #test_that("predict.CSC (strata): compare to mstate",{
  +         newdata <- data.frame(X1 = 0, X2 = 0, X16 = 0)
  +         newdata.L <- data.frame(X1.1 = c(0, 0), X1.2 = c(0, 0),
  +                                 X2.1 = c(0, 0), X2.2 = c(0, 0),
  +                                 X16.1 = c(0, 0), X16.2 = c(0, 0),
  +                                 trans = c(1, 2), strata = c(1, 2))
  +         e.coxph <- coxph(Surv(time, status) ~ X1.1 + X1.2 + X2.1 + X2.2 + X16.1 + X16.2 + strata(trans),
  +                          data = dL.exp)
  +         pred.msfit <- msfit(e.coxph, newdata = newdata.L, trans = tmat)
  +         suppressWarnings(
  +             pred.probtrans <- probtrans(pred.msfit,0)[[1]]
  +         )
  +         ## riskRegression no strata
  +         CSC.RR1 <- CSC(Hist(time,event)~X1+X2+X16, data = d, method = "breslow")
  +         pred.RR1a <- predict(CSC.RR1, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR1b <- predict(CSC.RR1, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         pred.RR1c <- predict(CSC.RR1, newdata, cause = 2, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR1d <- predict(CSC.RR1, newdata, cause = 2, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         CSC.RR2 <- CSC(Hist(time,event)~X1+X2+X16, data = d, surv.type = "survival", method = "breslow")
  +         pred.RR2a <- predict(CSC.RR2, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR2b <- predict(CSC.RR2, newdata, cause = 1, time = pred.probtrans[,"time"],
  +                              keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         ## riskRegression strata
  +         d2 <- rbind(cbind(d,grp=1),cbind(d,grp=2))
  +         CSC.RR1_strata <- CSC(Hist(time,event)~X1+X2+X16+strata(grp), data = d2, method = "breslow")
  +         pred.RR1a_strata <- predict(CSC.RR1_strata, cbind(newdata,grp=1), cause = 1, time = pred.probtrans[,"time"],
  +                                     keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR1b_strata <- predict(CSC.RR1_strata, cbind(newdata,grp=1), cause = 1, time = pred.probtrans[,"time"],
  +                                     keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         pred.RR1c_strata <- predict(CSC.RR1_strata, cbind(newdata,grp=1), cause = 2, time = pred.probtrans[,"time"],
  +                                     keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR1d_strata <- predict(CSC.RR1_strata, cbind(newdata,grp=1), cause = 2, time = pred.probtrans[,"time"],
  +                                     keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         expect_equal(pred.RR1a_strata$absRisk,pred.RR1a$absRisk)
  +         expect_equal(pred.RR1b_strata$absRisk,pred.RR1b$absRisk)
  +         expect_equal(pred.RR1c_strata$absRisk,pred.RR1c$absRisk)
  +         expect_equal(pred.RR1d_strata$absRisk,pred.RR1d$absRisk)
  +         CSC.RR2_strata<- CSC(Hist(time,event)~X1+X2+X16+strata(grp), data = d2, surv.type = "survival", method = "breslow")
  +         pred.RR2a_strata <- predict(CSC.RR2_strata, cbind(newdata,grp=1), cause = 1, time = pred.probtrans[,"time"],
  +                                     keep.newdata = FALSE, se = TRUE, product.limit = TRUE)
  +         pred.RR2b_strata <- predict(CSC.RR2_strata, cbind(newdata,grp=1), cause = 1, time = pred.probtrans[,"time"],
  +                                     keep.newdata = FALSE, se = TRUE, product.limit = FALSE)
  +         expect_equal(pred.RR2a_strata$absRisk,pred.RR2a$absRisk)
  +         expect_equal(pred.RR2b_strata$absRisk,pred.RR2b$absRisk)
  +     }
  + })
  [predictCSC] predict.CSC vs. mstate 
  ── Error (???): [predict.CSC]: mstate simulated data ───────────────────────────
  Error in ``[.data.table`(data, , tcols)`: j (the 2nd argument inside [...]) is a single symbol but column name 'tcols' is not found. Perhaps you intended DT[, ..tcols]. This difference to data.frame is deliberate and explained in FAQ 1.1.
  Backtrace:
   1. mstate::msprep(...)
   4. data.table:::`[.data.table`(data, , tcols)
  
  Error: Test failed
  Execution halted
Running the tests in ‘tests/test-predictCox.R’ failed.
Complete output:
  > ### test-predictCox.R --- 
  > #----------------------------------------------------------------------
  > ## author: Brice Ozenne
  > ## created: sep  4 2017 (10:38) 
  > ## Version: 
  > ## last-updated: Dec  7 2020 (09:23) 
  > ##           By: Thomas Alexander Gerds
  > ##     Update #: 168
  > #----------------------------------------------------------------------
  > ## 
  > ### Commentary: 
  > ##
  > ### Change Log:
  > #----------------------------------------------------------------------
  > ## 
  > ### Code:
  > 
  > ## * Settings
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(testthat)
  > library(rms)
  Loading required package: Hmisc
  Loading required package: lattice
  Loading required package: survival
  Loading required package: Formula
  Loading required package: ggplot2
  
  Attaching package: 'Hmisc'
  
  The following object is masked from 'package:testthat':
  
      describe
  
  The following objects are masked from 'package:base':
  
      format.pval, units
  
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > 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")
  [predictCox] Estimation of the baseline hazard (no strata) 
  > 
  > ## ** 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(predictCox(fit.coxph, centered = FALSE)$cumhazard, 
  +                survival::basehaz(fit.coxph, centered = FALSE)$hazard, tolerance = 1e-8)
  +   expect_equal(predictCox(fit.coxph, centered = TRUE)$cumhazard, 
  +                survival::basehaz(fit.coxph, centered = TRUE)$hazard, tolerance = 1e-8)
  +   expect_equal(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(predictCox(fit.cph),
  +                predictCox(fit.coxph, centered = TRUE), 
  +                tolerance = 100*max(abs(coef(fit.coxph)-coef(fit.cph))))
  + })
  ── Failure (???): baseline hazard (no strata): compare to survival::basehaz ────
  predictCox(fit.cph) not equal to predictCox(fit.coxph, centered = TRUE).
  Component "cumhazard": Mean relative difference: 0.8294145
  Component "survival": Mean relative difference: 0.2117415
  
  Error: Test failed
  Execution halted
Running the tests in ‘tests/test-predictRisk-TD.R’ failed.
Complete output:
  > ### test-predictRisk-TD.R --- 
  > #----------------------------------------------------------------------
  > ## Author: Thomas Alexander Gerds
  > ## Created: Sep  2 2017 (08:01) 
  > ## Version: 
  > ## Last-Updated: Dec  5 2020 (13:20) 
  > ##           By: Thomas Alexander Gerds
  > ##     Update #: 12
  > #----------------------------------------------------------------------
  > ## 
  > ### Commentary: 
  > ## 
  > ### Change Log:
  > #----------------------------------------------------------------------
  > ## 
  > ### Code:
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(testthat)
  > library(timereg)
  Loading required package: survival
  > library(survival)
  > test_that("start stop Cox: basehaz vs timereg",{
  +     data(veteran)
  +     ## 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(bh.cox[!duplicated(bh.cox$hazard),]$time,ca$cum[-1,"time"])
  +     expect_equal(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(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(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)))
  + })
  ── Warning (???): start stop Cox: basehaz vs timereg ───────────────────────────
  data set 'veteran' not found
  Backtrace:
   1. utils::data(veteran)
  
  > 
  > 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(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)
  +     expect_success(riskRegression:::predictRisk.CSCTD(fit,newdata=d[c(17:32),],landmark=2:5,times=1))
  + })
  ── Failure (???): start stop CauseSpecificCox ──────────────────────────────────
  no expectation used.
  
  Error: Test failed
  Execution halted
Running the tests in ‘tests/test-predictRisk.R’ failed.
Complete output:
  > ### test-predictRisk.R --- 
  > #----------------------------------------------------------------------
  > ## Author: Thomas Alexander Gerds
  > ## Created: Aug 10 2017 (08:56) 
  > ## Version: 
  > ## Last-Updated: Dec  5 2020 (13:20) 
  > ##           By: Thomas Alexander Gerds
  > ##     Update #: 33
  > #----------------------------------------------------------------------
  > ## 
  > ### Commentary: 
  > ## 
  > ### Change Log:
  > #----------------------------------------------------------------------
  > ## 
  > ### Code:
  > library(riskRegression)
  riskRegression version 2020.12.08
  > library(testthat)
  > library(rms)
  Loading required package: Hmisc
  Loading required package: lattice
  Loading required package: survival
  Loading required package: Formula
  Loading required package: ggplot2
  
  Attaching package: 'Hmisc'
  
  The following object is masked from 'package:testthat':
  
      describe
  
  The following objects are masked from 'package:base':
  
      format.pval, units
  
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > library(survival)
  > library(data.table)
  > library(lava)
  
  Attaching package: 'lava'
  
  The following object is masked from 'package:ggplot2':
  
      vars
  
  The following object is masked from 'package:testthat':
  
      compare
  
  > 
  > ## * SuperLearner
  > cat("[predictRisk.SuperPredictor] \n")
  [predictRisk.SuperPredictor] 
  > test_that("wrap the SuperLearner", {
  +     if (!requireNamespace("SuperLearner",quietly=TRUE)){
  +         message("Package SuperLearner not installed. Skip this test.")
  +     }else{
  +         library(SuperLearner)
  +         d <- sampleData(139,outcome="binary")
  +         sl = SuperLearner(Y = d$Y,
  +                           X = d[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10"),with=0L],
  +                           family = binomial(),
  +                           SL.library = "SL.glm")
  +     }
  + })
  ── Skip (???): wrap the SuperLearner ───────────────────────────────────────────
  Reason: empty test
  
  > 
  > ## * rfsrc
  > cat("[predictRisk.rfsrc]  \n")
  [predictRisk.rfsrc]  
  > test_that("Additional arguments: example with imputation of missing data", {
  +     if (!requireNamespace("randomForestSRC",quietly=TRUE)){
  +         message("Package randomForestSRC not installed. Skip this test.")
  +     }else{
  +         data(pbc,package="survival")
  +         set.seed(10)
  +         forest <- rfsrc(Surv(time,status)~chol+age+sex,data=pbc,ntree=10,nsplit=10)
  +         ## no specification
  +         expect_error(Score(list(forest),formula=Hist(time,status)~1,data=pbc,conf.int=FALSE,metrics="brier",cens.model="km"))
  +         ## correct specification
  +         expect_output(print(Score(list(forest),formula=Hist(time,status)~1,data=pbc,conf.int=FALSE,metrics="brier",cens.model="km",predictRisk.args=list("rfsrc"=list(na.action="na.impute")))))
  +         ## wrong specification
  +         expect_error(print(Score(list(forest),formula=Hist(time,status)~1,data=pbc,conf.int=FALSE,metrics="brier",cens.model="km",predictRisk.args=list("randomForestSRC"=list(na.action="na.impute")))))
  +     }
  + })
  ── Error (???): Additional arguments: example with imputation of missing data ──
  Error in `rfsrc(Surv(time, status) ~ chol + age + sex, data = pbc, ntree = 10, 
      nsplit = 10)`: could not find function "rfsrc"
  
  Error: Test failed
  Execution halted