tests/UnitTestsPlots.R

In discSurv: Discrete Time Survival Analysis

library(discSurv)

############################
# Calibration
# Simulation with true model

set.seed(12)
sampleSize <- 250
xTrain <- rnorm(n=sampleSize)
ObsMat <- t(sapply(X=1:sampleSize, FUN=function(i) {
  
  # Censoring time
  censT <- sample(1:5, size=1, replace=FALSE,
                  prob=c(4/10, 3/10, 2/10, 1/10, 0))
  #censT <- 5
  
  survT <- 1
  defProb <- exp(xTrain[i]+0.1*survT)/(1+exp(xTrain[i]+0.1*survT))
  randNumb <- rbinom(n=1, size=1, prob=defProb)
  while( randNumb!=1 & survT < 5 ){
    survT <- survT + 1
    defProb <- exp(xTrain[i]+0.1*survT)/(1+exp(xTrain[i]+0.1*survT))
    randNumb <- rbinom(n=1, size=1, prob=defProb)
  }
  
  if( survT <= censT ){
    return(c(event=1, obsTime=survT))
  } else{
    return(c(event=0, obsTime=censT))
  }
  
} ))
datShort <- data.frame(ObsMat, x=xTrain)
datLong <- dataLong(dataShort=datShort, eventColumn = "event",
                    timeColumn = "obsTime", timeAsFactor=TRUE)

# Test data
sampleSize <- 250
xTest <- rnorm(n=sampleSize)
ObsMatTest <- t(sapply(X=1:sampleSize, FUN=function(i) {
  
  # Censoring time
  censT <- sample(1:5, size=1, replace=FALSE,
                  prob=c(4/10, 3/10, 2/10, 1/10, 0))
  #censT <- 5
  
  survT <- 1
  defProb <- exp(xTest[i]+0.1*survT)/(1+exp(xTest[i]+0.1*survT))
  randNumb <- rbinom(n=1, size=1, prob=defProb)
  while( randNumb!=1 & survT < 5 ){
    survT <- survT + 1
    defProb <- exp(xTest[i]+0.1*survT)/(1+exp(xTest[i]+0.1*survT))
    randNumb <- rbinom(n=1, size=1, prob=defProb)
  }
  
  if( survT <= censT ){
    return(c(event=1, obsTime=survT))
  } else{
    return(c(event=0, obsTime=censT))
  }
  
} ))
datShortTest <- data.frame(ObsMatTest, x=xTest)
datLongTest <- dataLong(dataShort=datShortTest, eventColumn = "event",
                        timeColumn = "obsTime", timeAsFactor=TRUE)

# Estimate model
estMod <- glm(formula=y ~ timeInt + x, family=binomial(), data=datLong)
estHazLongTest <- predict(estMod, newdata=datLongTest, type="response")
trueHazLongTest <- exp(datLongTest$x+0.1*as.numeric(datLongTest$timeInt))/
  (1+exp(datLongTest$x+0.1*as.numeric(datLongTest$timeInt)))

# Calibration plot
calPlot(estHazards=estHazLongTest, testDataLong=datLongTest)

# Comparison with true hazard rates
mean(abs(estHazLongTest-trueHazLongTest)) <= 0.1

# Estimate recalibration model
estLinPredLongTest <- predict(estMod, newdata=datLongTest, type="link")
m1 <- estRecal(testLinPred=estLinPredLongTest, testDataLong=datLongTest)
summary(m1)
# Test 1: alpha=0 and beta=0
attr(m1, "Test1")
# Test 2: alpha=0 | beta=1
attr(m1, "Test2")
# Test 3: beta=1 | alpha
attr(m1, "Test3")