Running ‘test-BinomialRegression.R’ [0s/1s]
  Running ‘test-BrierScore.R’ [9s/14s]
  Running ‘test-CauseSpecificCoxRegression.R’ [4s/6s]
  Running ‘test-FineGrayRegression.R’ [3s/5s]
  Running ‘test-ate.R’ [17s/39s]
  Running ‘test-boxplotScore.R’ [3s/4s]
  Running ‘test-confidenceBand_vs_timereg.R’ [6s/7s]
  Running ‘test-iidCox.R’ [3s/5s]
  Running ‘test-phreg.R’ [0s/1s]
  Running ‘test-predictCSC_vs_mstate.R’ [4s/5s]
Running the tests in ‘tests/test-predictCSC_vs_mstate.R’ failed.
Complete output:
  > ### test-predictCSC_vs_mstate.R --- 
  > #----------------------------------------------------------------------
  > ## author: Brice Ozenne
  > ## created: maj 18 2017 (09:23) 
  > ## Version: 
  > ## last-updated: jun  8 2017 (20:32) 
  > ##           By: Brice Ozenne
  > ##     Update #: 137
  > #----------------------------------------------------------------------
  > ## 
  > ### Commentary: 
  > ## Compare predictCSC with mstate
  > ### Change Log:
  > #----------------------------------------------------------------------
  > ## 
  > ### Code:
  > 
  > library(riskRegression)
  Loading required package: data.table
  Loading required package: ggplot2
  Loading required package: prodlim
  > library(testthat)
  > library(mstate)
  Loading required package: survival
  > tmat <- trans.comprisk(2, names = c("0", "1", "2"))
  > library(survival)
  > 
  > 
  > 
  > # {{{ 1- compare predict.CSC and mstate on small examples
  > 
  > ## 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)
  > 
  > dfS <- rbind(cbind(df1, grp = 1, X2 = 0),
  +              cbind(rbind(df1,df1),grp = 2, X2 = 0),
  +              cbind(df1, grp = 3, X2 = 0)
  +              )
  > 
  > ## 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)
  > 
  > # {{{ 1a- no risk factor no strata
  > 
  > test_that("predict.CSC/predictCox (1a,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, productLimit = FALSE)  
  +     pred.exp2 <- predict(CSC.exp, times = 1:10, newdata = data.frame(NA), cause = 1, productLimit = TRUE) 
  +  
  +     CSC.prodlim <- CSC(Hist(time,event) ~ 1, data = df1, survtype = "survival", method = "breslow")
  +     pred.prodlim <- predict(CSC.prodlim, times = 1:10, newdata = data.frame(NA), cause = 1, productLimit = TRUE)  
  +     expect_equal(pred.prodlim,pred.exp2)
  +     
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ 1, data = df1, survtype = "survival", method = "efron")
  +     pred.prodlimE <- predict(CSC.prodlimE, times = 1:10, newdata = data.frame(NA), cause = 1, productLimit = 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_that("predict.CSC/predictCox (1a,df2) - 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, productLimit = FALSE)  
  +     pred.exp2 <- predict(CSC.exp, times = 1:10, newdata = data.frame(NA), cause = 1, productLimit = TRUE)  
  + 
  +     CSC.prodlim <- CSC(Hist(time,event) ~ 1, data = df2, survtype = "survival", method = "breslow")
  +     pred.prodlim <- predict(CSC.prodlim, times = 1:10, newdata = data.frame(NA), cause = 1, productLimit = TRUE)  
  +     expect_equal(pred.exp2, pred.prodlim)
  + 
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ 1, data = df2, survtype = "survival", method = "efron")
  +     pred.prodlimE <- predict(CSC.prodlimE, times = 1:10, newdata = data.frame(NA), cause = 1, productLimit = 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_that("predict.CSC (1a) - compare to 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
  +         dL <- 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,
  +                            productLimit = TRUE)
  +         
  +         expect_equal(pred.probtrans[,"pstate2"],
  +                      as.double(pred.RR$absRisk)
  +                      )
  +     }
  + })
  > 
  > test_that("predict.CSC(1a) - 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, productLimit = TRUE)  
  +         haz2.prodlim <- predict(CSC.exp, times = seqTimes, newdata = newdata, cause = 2, productLimit = TRUE)  
  + 
  +         expect_equal(surv.prodlim+as.double(haz1.prodlim$absRisk)+as.double(haz2.prodlim$absRisk),rep(1,nTimes))
  +     }
  + })
  > # }}}
  > 
  > # {{{ 1b- no risk factor strata
  > 
  > test_that("predict.CSC(1b) - compare to manual estimation",{
  +     # don't work:
  +     # CSC.exp <- CSC(Hist(time,event) ~ strata(grp), data = dfS)
  +     # X2 factive variable
  +     CSC.exp <- CSC(Hist(time,event) ~ strata(grp) + X2, data = dfS, method = "breslow")
  +     pred.exp <- predict(CSC.exp, times = 1:10, newdata = data.frame(grp = 1:3, X2 = 0), cause = 1, productLimit = FALSE)  
  +     pred.exp2 <- predict(CSC.exp, times = 1:10, newdata = data.frame(grp = 1:3, X2 = 0), cause = 1, productLimit = TRUE)  
  + 
  +     CSC.prodlim <- CSC(Hist(time,event) ~ strata(grp) + X2, data = dfS, survtype = "survival", method = "breslow")
  +     pred.prodlim <- predict(CSC.prodlim, times = 1:10, newdata = data.frame(grp = 1:3, X2 = 0), cause = 1, productLimit = TRUE)  
  +     expect_equal(pred.exp2,pred.prodlim)
  +     
  +     CSC.prodlimE <- CSC(Hist(time,event) ~ strata(grp) + X2, data = dfS, survtype = "survival", method = "efron")
  +     pred.prodlimE <- predict(CSC.prodlimE, times = 1:10, newdata = data.frame(grp = 1:3, X2 = 0), cause = 1, productLimit = TRUE)  
  + 
  +     ## baseline
  +     lambda1 <- lambda2 <- 1/seq(20,2,by = -2)
  +     lambda2E.1 <- as.data.table(predictCox(CSC.prodlimE$models[[1]], type = "hazard"))
  +     lambda2E.2 <- as.data.table(predictCox(CSC.prodlimE$models[[2]], type = "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.2[strata=="grp=1"][["hazard"]]))[1:10]
  +     expect_equal(as.double(pred.prodlimE$absRisk[pred.prodlimE$strata=="grp=1"]),
  +                  cumsum(lambda1*survival))
  +     survival <- c(1,cumprod(1-lambda2E.2[strata=="grp=2"][["hazard"]]))[1:10]
  +     expect_equal(as.double(pred.prodlimE$absRisk[pred.prodlimE$strata=="grp=2"]),
  +                  cumsum(lambda2E.1[strata=="grp=2"][["hazard"]]*survival))    
  +     expect_equal(as.double(pred.prodlimE$absRisk[1,]),
  +                  as.double(pred.prodlimE$absRisk[3,]))
  + })
  Error: Test failed: 'predict.CSC(1b) - compare to manual estimation'
  * as.double(pred.exp$absRisk[1, ]) not equal to cumsum(lambda1 * survival).
  10/10 mismatches (average diff: NaN)
  [1] NA - 0.050 == NA
  [2] NA - 0.100 == NA
  [3] NA - 0.151 == NA
  [4] NA - 0.202 == NA
  [5] NA - 0.254 == NA
  [6] NA - 0.306 == NA
  [7] NA - 0.360 == NA
  [8] NA - 0.415 == NA
  [9] NA - 0.475 == NA
  ...
  * as.double(pred.prodlim$absRisk[1, ]) not equal to cumsum(lambda1 * survival).
  10/10 mismatches (average diff: NaN)
  [1] NA - 0.05 == NA
  [2] NA - 0.10 == NA
  [3] NA - 0.15 == NA
  [4] NA - 0.20 == NA
  [5] NA - 0.25 == NA
  [6] NA - 0.30 == NA
  [7] NA - 0.35 == NA
  [8] NA - 0.40 == NA
  [9] NA - 0.45 == NA
  ...
  * as.double(pred.prodlimE$absRisk[pred.prodlimE$strata == "grp=1"]) not equal to cumsum(lambda1 * survival).
  10/10 mismatches (average diff: NaN)
  [1] NA - 0.0500 == NA
  [2] NA - 0.0999 == NA
  [3] NA - 0.1495 == NA
  [4] NA - 0.1990 == NA
  [5] NA - 0.2481 == NA
  [6] NA - 0.2967 == NA
  [7] NA - 0.3447 == NA
  [8] NA - 0.3916 == NA
  [
  Execution halted