### test-iidCox.R --- 
#----------------------------------------------------------------------
## author: Brice Ozenne
## created: maj 18 2017 (09:23) 
## Version: 
## last-updated: Sep 17 2022 (07:00) 
##           By: Thomas Alexander Gerds
##     Update #: 125
#----------------------------------------------------------------------
## 
### Commentary: 
## Check internal consistency of the iid decomposition (e.g. after last event time, IF cumhazard vs. IF hazard)
## Compare the iid decomposition (lambda and beta) obtained with iidCox and timereg
##
### Change Log:
#----------------------------------------------------------------------
## 
### Code:

## * Settings
library(riskRegression)
library(testthat)
library(survival)
library(rms)
library(timereg)
library(data.table)

context("function iidCox")

## * Internal consistency
cat("[iidCox] internal consistency \n")

## ** Data
set.seed(10)
dt <- sampleData(20, outcome = "survival")[,.(time,event,X1,X2)]

timeFirstEvent <- dt[event==1,min(time)]
timeLastEvent <- dt[event==1,max(time)]
timeLastObs <- dt[,max(time)]
vec.times <-  timeFirstEvent+1:5
dt.newdata <- dt[1:4]

## ** Model
coxph.fit <- coxph(Surv(time,event==1)~ X1+X2,
                   data=dt, x = TRUE, y = TRUE)
cph.fit <- coxph(Surv(time,event==1)~ X1+X2,
                 data=dt, x = TRUE, y = TRUE)

## ** Tests
iid.coxph <-  iidCox(coxph.fit, return.object = FALSE)

test_that("[iidCox] consistency coxph,cph",{
    iid.cph <-  iidCox(cph.fit, return.object = FALSE)
    expect_equal(ignore_attr=TRUE,iid.cph, iid.coxph)
})

test_that("[iidCox] consistency with manual specification of newdata and times",{
    iid.test <- iidCox(coxph.fit, newdata = dt, tau.hazard = iid.coxph$time[[1]], return.object = FALSE)
    expect_equal(ignore_attr=TRUE,iid.coxph, iid.test)

    new.order <- sample.int(length(iid.coxph$time[[1]]))
    iid.test <- iidCox(coxph.fit, newdata = dt, tau.hazard = iid.coxph$time[[1]][new.order], return.object = FALSE)
    expect_equal(ignore_attr=TRUE,iid.coxph$IFhazard[[1]], iid.test$IFhazard[[1]][,colnames(iid.coxph$IFhazard[[1]])])
    expect_equal(ignore_attr=TRUE,iid.coxph$IFcumhazard[[1]], iid.test$IFcumhazard[[1]][,colnames(iid.coxph$IFcumhazard[[1]])])
    expect_equal(ignore_attr=TRUE,iid.coxph$IFbeta, iid.coxph$IFbeta)
})

test_that("[iidCox] cumsum(iid hazard) = iid cumhazard",{
    M.test <- t(apply(iid.coxph$IFhazard[[1]], 1, cumsum))
    expect_equal(ignore_attr=TRUE,M.test, iid.coxph$IFcumhazard[[1]])
})

test_that("[iidCox] iid hazard = 0 at non event times and NA after last observation",{
    iid.test <- iidCox(coxph.fit, tau.hazard = sort(c(dt$time,c(0.1,1,5,8,12,25))), return.object = FALSE)

    expect_equal(ignore_attr=TRUE,iid.test$IFhazard[[1]][,as.character(iid.coxph$time[[1]])],
                 iid.coxph$IFhazard[[1]])

    vec.noEventTime <- setdiff(iid.test$time[[1]],dt[event==1,time])
    vec.noEventTime_beforeTmax <- vec.noEventTime[vec.noEventTime<=timeLastObs]
    vec.noEventTime_afterTmax <- vec.noEventTime[vec.noEventTime>timeLastObs]
    
    expect_true(all(iid.test$IFhazard[[1]][,as.character(vec.noEventTime_beforeTmax)]==0))
    expect_true(all(is.na(iid.test$IFhazard[[1]][,as.character(vec.noEventTime_afterTmax)])))

    ## additional tests
    ## compatibility IFhazard, IFcumhazard
    M.test <- t(apply(iid.test$IFhazard[[1]], 1, cumsum))
    expect_equal(ignore_attr=TRUE,M.test, iid.test$IFcumhazard[[1]])

    ## iid cumHazard[lastEvent+] =  cumHazard[lastEvent]
    vec.time <- intersect(iid.test$time[[1]][iid.test$time[[1]]>timeLastEvent],
                          iid.test$time[[1]][iid.test$time[[1]]<=timeLastObs])
    expect_true(all(iid.test$IFcumhazard[[1]][,as.character(vec.time)]-iid.test$IFcumhazard[[1]][,as.character(timeLastEvent)]==0))
})
## ** selectJump
## test_that("[iidCox] selectJump",{
##     seqTest <- c(0,dt$time[1:10],dt$time[1:10]-1e-12,dt$time[1:10]+1e-5,1e5,1)
##     GS <-  iidCox(coxph.fit,
##                   tau.hazard = seqTest,
##                   return.object = FALSE)
##     test <- riskRegression:::selectJump(iid.coxph, times = seqTest, type = c("hazard","cumhazard"))

##     expect_equal(ignore_attr=TRUE,unname(GS$IFhazard[[1]]),unname(test$IFhazard[[1]]))
##     expect_equal(ignore_attr=TRUE,unname(GS$IFcumhazard[[1]]),unname(test$IFcumhazard[[1]]))
## })

## ** Empty strata
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))

e.coxph <- coxph(Surv(time, event) ~ strata(X1), data = dtS,
                 x = TRUE, y = TRUE)
test_that("[iidCox] empty strata", {
    test <- iidCox(e.coxph, tau.max = 1)
    expect_equal(ignore_attr=TRUE,test$iid$time[[1]],0)
    expect_equal(ignore_attr=TRUE,test$iid$time[[2]],c(0.2,0.6,0.9))
    expect_equal(ignore_attr=TRUE,test$iid$time[[3]],0)

    expect_equal(ignore_attr=TRUE,unname(test$iid$IFhazard[[1]]), matrix(0, nrow = NROW(dtS), ncol = 1))
    GS <- cbind("0.2" = c(0, 0, 0, -0.00173611, 0.03993056, 0, 0, -0.00173611, 0, 0, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, -0.00173611, 0, 0, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, -0.00173611, 0, -0.00173611, -0.00173611, -0.00173611, -0.00173611, 0, 0, -0.00173611, 0, 0, -0.00173611, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, -0.00173611), 
                "0.6" = c(0, 0, 0, -0.00189036, 0, 0, 0, -0.00189036, 0, 0, 0, 0, -0.00189036, -0.00189036, -0.00189036, 0, -0.00189036, 0, 0, 0, 0, -0.00189036, -0.00189036, -0.00189036, 0, 0, 0, -0.00189036, -0.00189036, -0.00189036, 0, -0.00189036, 0, -0.00189036, -0.00189036, -0.00189036, -0.00189036, 0, 0, -0.00189036, 0, 0, -0.00189036, 0, 0, -0.00189036, 0.0415879, -0.00189036, 0, -0.00189036), 
                "0.9" = c(0, 0, 0, -0.00206612, 0, 0, 0, -0.00206612, 0, 0, 0, 0, -0.00206612, -0.00206612, -0.00206612, 0, -0.00206612, 0, 0, 0, 0, 0.04338843, -0.00206612, -0.00206612, 0, 0, 0, -0.00206612, -0.00206612, -0.00206612, 0, -0.00206612, 0, -0.00206612, -0.00206612, -0.00206612, -0.00206612, 0, 0, -0.00206612, 0, 0, -0.00206612, 0, 0, -0.00206612, 0, -0.00206612, 0, -0.00206612))
    expect_equal(ignore_attr=TRUE,test$iid$IFhazard[[2]], GS, tolerance = 1e-6)
    expect_equal(ignore_attr=TRUE,unname(test$iid$IFhazard[[3]]), matrix(0, nrow = NROW(dtS), ncol = 1))

    expect_equal(ignore_attr=TRUE,unname(test$iid$IFcumhazard[[1]]), matrix(0, nrow = NROW(dtS), ncol = 1))
    GS <- cbind("0.2" = c(0, 0, 0, -0.00173611, 0.03993056, 0, 0, -0.00173611, 0, 0, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, -0.00173611, 0, 0, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, -0.00173611, 0, -0.00173611, -0.00173611, -0.00173611, -0.00173611, 0, 0, -0.00173611, 0, 0, -0.00173611, 0, 0, -0.00173611, -0.00173611, -0.00173611, 0, -0.00173611), 
                "0.6" = c(0, 0, 0, -0.00362647, 0.03993056, 0, 0, -0.00362647, 0, 0, 0, 0, -0.00362647, -0.00362647, -0.00362647, 0, -0.00362647, 0, 0, 0, 0, -0.00362647, -0.00362647, -0.00362647, 0, 0, 0, -0.00362647, -0.00362647, -0.00362647, 0, -0.00362647, 0, -0.00362647, -0.00362647, -0.00362647, -0.00362647, 0, 0, -0.00362647, 0, 0, -0.00362647, 0, 0, -0.00362647, 0.03985179, -0.00362647, 0, -0.00362647), 
                "0.9" = c(0, 0, 0, -0.00569259, 0.03993056, 0, 0, -0.00569259, 0, 0, 0, 0, -0.00569259, -0.00569259, -0.00569259, 0, -0.00569259, 0, 0, 0, 0, 0.03976196, -0.00569259, -0.00569259, 0, 0, 0, -0.00569259, -0.00569259, -0.00569259, 0, -0.00569259, 0, -0.00569259, -0.00569259, -0.00569259, -0.00569259, 0, 0, -0.00569259, 0, 0, -0.00569259, 0, 0, -0.00569259, 0.03985179, -0.00569259, 0, -0.00569259))
    expect_equal(ignore_attr=TRUE,test$iid$IFcumhazard[[2]], GS, tolerance = 1e-6)
    expect_equal(ignore_attr=TRUE,unname(test$iid$IFcumhazard[[3]]), matrix(0, nrow = NROW(dtS), ncol = 1))

    expect_equal(ignore_attr=TRUE,test$iid$etime.max, c(11.4,14.8,9.2))
    ## expect_equal(ignore_attr=TRUE,test$iid$indexObs,order(dtS$time))
})

## * Compare to timereg
## ** Data
data(Melanoma, package = "riskRegression")

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]
dt[ , Xcat2 := as.factor(paste0(X1,X2))]

## sorted dataset
dt.sort <- copy(dt)
setkeyv(dt.sort,c("time")) 

## dataset with ties
dtTies <- copy(dt)[1:10,]
dtTies[, event := 1]
dtTies[7:8, X1 := 1]
dtTies[7:8, X6 := 1]
setkeyv(dtTies, c("time"))
dtTies[, timeTies := time]
dtTies[7:8, timeTies := time[1]]  # ties

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"))

## ** No strata, no interaction, continous
cat("[iidCox] compare to timereg - no strata, no interaction, continuous \n")
## *** Models
e.coxph <- coxph(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),
                          data = dt, resample.iid = TRUE, max.timepoint.sim=NULL)

eApprox.timereg <- cox.aalen(Surv(time, event) ~ prop(X1)+prop(X6),
                                data = dt, resample.iid = TRUE, max.timepoint.sim=NULL,
                                rate.sim = 0, clusters = 1:NROW(dt))

## *** Extract information
## compute the IF of the baseline hazard and then for the survival
IF.coxph <- iidCox(e.coxph, keep.times = FALSE, store.iid = "full", return.object = FALSE)
IF.coxph_sort <- iidCox(e.coxph_sort, keep.times = FALSE, store.iid = "full", return.object = FALSE)
## store "sufficient statistics" for the IF of the baseline hazard
## and then compute the IF for the survival
IFminimal.coxph <- iidCox(e.coxph, keep.times = FALSE, store.iid = "minimal", return.object = FALSE)
## same as before but using an approximation

IFlambda_GS <- t(as.data.table(e.timereg$B.iid))
 
## *** Tests
test_that("[iidCox] beta - no strata, no interaction, continuous",{
    expect_equal(ignore_attr=TRUE,unname(IF.coxph$IFbeta),e.timereg$gamma.iid)
    expect_equal(ignore_attr=TRUE,unname(IF.coxph$IFbeta[order(dt$time),]),unname(IF.coxph_sort$IFbeta))
    expect_equal(ignore_attr=TRUE,unname(IFminimal.coxph$IFbeta),e.timereg$gamma.iid)
})

test_that("[iidCox] lambda - no strata, no interaction, continuous",{
    expect_equal(ignore_attr=TRUE,as.double(IF.coxph$IFcumhazard[[1]]), as.double(IFlambda_GS[,-1]))
    expect_equal(ignore_attr=TRUE,IF.coxph$IFcumhazard[[1]][order(dt$time),], IF.coxph_sort$IFcumhazard[[1]])
    expect_equal(ignore_attr=TRUE,IFminimal.coxph$IFcumhazard[[1]], NULL)
})
  

## ** no strata, interactions, continous
cat("[iidCox] compare to timereg - no strata, interaction, continuous \n")

## *** Models
e.coxph <- coxph(Surv(time, event) ~ X1*X6, data = dt, 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)

## *** Extract information
IF.coxph <- iidCox(e.coxph, keep.times = FALSE, return.object = FALSE)
IFlambda_GS <- t(as.data.table(e.timereg$B.iid))

## *** Tests
test_that("[iidCox] beta - no strata, interactions, continuous",{
    expect_equal(ignore_attr=TRUE,unname(IF.coxph$IFbeta),e.timereg$gamma.iid)
})
test_that("[iidCox] lambda - no strata, interactions, continuous",{
    expect_equal(ignore_attr=TRUE,as.double(IF.coxph$IFcumhazard[[1]]), as.double(IFlambda_GS[,-1]))
})


## ** no covariate
cat("[iidCox] compare to timereg - no covariate \n")
## cox.aalen do not work without covariable
## cox.aalen(Surv(eventtime, event) ~ 1, data = dt, resample.iid = TRUE, max.timepoint.sim=NULL)
## Error in model.frame.default(formula = Surv(survs$stop, survs$status) ~  : 
##   type (NULL) incorrect pour la variable 'Z'

test_that("[iidCox] beta - no covariate",{
    e.coxph <- coxph(Surv(time, event) ~ 1, data = dt, y = TRUE, x = TRUE)
    IF.coxph <- iidCox(e.coxph, keep.times = FALSE, return.object = FALSE)
    expect_true(all(is.na(IF.coxph$IFbeta)))
})

## ** no strata, no interaction, with a categorical variable
cat("[iidCox] compare to timereg - no strata, no interaction, categorical \n")

## *** Models
e.coxph <- coxph(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)

## coef(e.coxph)
## coef(e.timereg)

## *** Extract information
IFlambda_GS <- t(as.data.table(e.timereg$B.iid))
IF.coxph <- iidCox(e.coxph, keep.times = FALSE, return.object = FALSE)

## *** Tests
test_that("[iidCox] beta - no strata, interactions, categorical",{
    expect_equal(ignore_attr=TRUE,unname(IF.coxph$IFbeta),e.timereg$gamma.iid)
})
  
test_that("[iidCox] lambda - no strata, interactions, categorical",{
    expect_equal(ignore_attr=TRUE,as.double(IF.coxph$IFcumhazard[[1]]),
                 as.double(IFlambda_GS[,-1]))
})



## ** Ties
cat("[iidCox] compare to timereg - ties \n")

e.coxph_Efron <- coxph(Surv(timeTies, event) ~ X1+X6, data = dtTies, y = TRUE, x = TRUE,
                 ties = "efron")
e.coxph_Breslow <- coxph(Surv(timeTies, event) ~ X1+X6, data = dtTies, y = TRUE, x = TRUE,
                         ties = "breslow")
e.timereg <- cox.aalen(Surv(timeTies, event) ~ prop(X1)+prop(X6), data = dtTies, resample.iid = TRUE, max.timepoint.sim=NULL)

## different coef for timereg - cannot compare
coef(e.coxph_Breslow)
coef(e.coxph_Efron)
coef(e.timereg)[,"Coef."]

## *** Extract information
IF.coxph_Efron <- iidCox(e.coxph_Efron, keep.times = FALSE, return.object = FALSE)
IF.coxph_Breslow <- iidCox(e.coxph_Breslow, keep.times = FALSE, return.object = FALSE)

## *** Tests
## ??? Gold standard


## ** Strata, no interaction, continuous
cat("[iidCox] compare to timereg - strata \n")
## *** Model
e.timereg <- cox.aalen(Surv(time, event) ~ strata(strata)-1 + prop(X1) + prop(X6), data = dtStrata, 
                       resample.iid = TRUE, max.timepoint.sim=NULL)
e.coxph <- coxph(Surv(time, event) ~ strata(strata) + X1 + X6, data = dtStrata, y = TRUE, x = TRUE)

## *** Extract information
IF.coxph <- iidCox(e.coxph, return.object = FALSE)

## *** Tests
test_that("[iidCox] beta - strata",{
    expect_equal(ignore_attr=TRUE,unname(IF.coxph$IFbeta),e.timereg$gamma.iid)
})
  
test_that("[iidCox] lambda - strata",{

    name.strata <- unique(dtStrata$strata)
    n.strata <- length(name.strata)
    
    for(iStrata in 1:n.strata){ ## iStrata <- 1
        IF.GS <- do.call(rbind,
                         lapply(e.timereg$B.iid,function(x){x[,iStrata]})
                         )
        colnames(IF.GS) <- e.timereg$time.sim.resolution
      
        checkTimes <- intersect(e.timereg$time.sim.resolution,
                                IF.coxph$time[[iStrata]])
      
        term1 <- IF.coxph$IFcumhazard[[iStrata]][,which(IF.coxph$time[[iStrata]] %in% checkTimes),drop = FALSE]
        term2 <- IF.GS[,which(e.timereg$time.sim.resolution %in% checkTimes)]
        diff <- term1 - term2
        expect_true(all(abs(na.omit(as.double(diff)))<1e-10))
    }
    
  })



## ** Melanoma data
cat("[iidCox] compare to timereg - Melanoma \n")
data(Melanoma, package = "riskRegression")

test_that("[iidCox] Compare to timereg on Melanoma dta",{
    e.timereg <- cox.aalen(Surv(time,status==1)~prop(sex), data = Melanoma)
    e.coxph <- coxph(Surv(time,status==1)~sex, data=Melanoma, x=TRUE, y=TRUE)

    RR.iid <- iidCox(e.coxph, return.object = FALSE)
    timereg.iidLambda <- t(as.data.table(e.timereg$B.iid))

    expect_equal(ignore_attr=TRUE,unname(RR.iid$IFbeta),e.timereg$gamma.iid)
    expect_equal(ignore_attr=TRUE,as.double(RR.iid$IFcumhazard[[1]]), ## as.double(RR.iid$IFcumhazard[[1]])[c(1:4,206:209,411)]
                 as.double(timereg.iidLambda[,-1])) ## as.double(timereg.iidLambda[,-1])[c(1:4,206:209,411)]
})
#----------------------------------------------------------------------
### test-iidCox.R ends here