### predictCoxPL.R ---
#----------------------------------------------------------------------
## author: Brice Ozenne
## created: sep 4 2017 (16:43)
## Version:
## last-updated: aug 28 2018 (10:42)
## By: Brice Ozenne
## Update #: 151
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
## * predictCoxPL (documentation)
#' @title Computation of survival probabilities from Cox regression models using the product limit estimator.
#' @description Same as predictCox except that the survival is estimated using the product limit estimator.
#' @name predictCoxPL
#'
#' @inheritParams predictCox
#' @param ... additional arguments to be passed to \code{\link{predictCox}}.
#'
## * predictCoxPL (code)
#' @rdname predictCoxPL
#' @examples
#' library(survival)
#'
#' #### generate data ####
#' set.seed(10)
#' d <- sampleData(40,outcome="survival")
#' nd <- sampleData(4,outcome="survival")
#' d$time <- round(d$time,1)
#'
#' #### Cox model ####
#' fit <- coxph(Surv(time,event)~ X1 + X2 + X6,
#' data=d, ties="breslow", x = TRUE, y = TRUE)
#'
#' ## exponential approximation
#' predictCox(fit, newdata = d, times = 1:5)
#'
#' ## product limit
#' predictCoxPL(fit, newdata = d, times = 1:5)
#'
#' #### stratified Cox model ####
#' fitS <- coxph(Surv(time,event)~ X1 + strata(X2) + X6,
#' data=d, ties="breslow", x = TRUE, y = TRUE)
#'
#' ## exponential approximation
#' predictCox(fitS, newdata = d, times = 1:5)
#'
#' ## product limit
#' predictCoxPL(fitS, newdata = d, times = 1:5)
#'
#' #### fully stratified Cox model ####
#' fitS <- coxph(Surv(time,event)~ 1,
#' data=d, ties="breslow", x = TRUE, y = TRUE)
#'
#' ## product limit
#' GS <- survfit(Surv(time,event)~1, data = d)
#' range(predictCoxPL(fitS)$survival - GS$surv)
#'
#' fitS <- coxph(Surv(time,event)~ strata(X2),
#' data=d, ties="breslow", x = TRUE, y = TRUE)
#'
#' ## product limit
#' GS <- survfit(Surv(time,event)~X2, data = d)
#' range(predictCoxPL(fitS)$survival - GS$surv)
#'
## * predictCoxPL (code)
#' @rdname predictCoxPL
#' @export
predictCoxPL <- function(object,
times,
newdata = NULL,
type = c("cumhazard","survival"),
keep.strata = TRUE,
keep.infoVar = FALSE,
...){
stop <- NULL ## [:CRANcheck:] data.table
## ** normalize arguments
object.modelFrame <- coxModelFrame(object)
if(!is.null(newdata)){
if(is.data.table(newdata)){
newdata <- copy(newdata)
}else{
newdata <- as.data.table(newdata)
}
}else if (missing(times)) {
times <- numeric(0)
}
if("survival" %in% type == FALSE){
stop("The argument \'type\' must contain \"survival\" \n",
"use the function predictCox otherwise \n")
}
## ** run original prediction
original.res <- predictCox(object = object,
newdata = newdata,
times = times,
type = type,
keep.strata = TRUE,
keep.infoVar = TRUE,
...)
infoVar <- original.res$infoVar
if(keep.infoVar==FALSE){
original.res$infoVar <- NULL
}
## ** compute survival
if(is.null(newdata)){
hazard.res <- as.data.table(predictCox(object = object,
type = "hazard",
keep.strata = TRUE,
...))
if(!is.null(hazard.res$strata)){
Ustrata <- levels(hazard.res$strata)
n.Ustrata <- length(Ustrata)
for(iS in 1:n.Ustrata){ ## iS <- 1
iIndex.all <- which(hazard.res$strata==Ustrata[iS])
iIndex.times <- which(original.res$strata==Ustrata[iS])
vec.surv <- c(1,cumprod(1-hazard.res$hazard[iIndex.all]))
index.jump <- prodlim::sindex(eval.times = original.res$times[iIndex.times],
jump.times = c(0,hazard.res$times[iIndex.all]))
original.res$survival[iIndex.times] <- vec.surv[index.jump]
}
if(keep.strata == FALSE){
original.res$strata <- NULL
}
}else{
vec.surv <- c(1,cumprod(1-hazard.res$hazard))
index.jump <- prodlim::sindex(eval.times = original.res$times,
jump.times = c(0,hazard.res$times))
original.res$survival <- vec.surv[index.jump]
}
if("hazard" %in% type == FALSE){
original.res$hazard <- NULL
}
}else if(infoVar$is.strata){
object.strata <- coxStrata(object,
data = NULL,
strata.vars = infoVar$stratavars)
new.strata <- original.res$strata
if(keep.strata == FALSE){
original.res$strata <- NULL
}
Ustrata <- unique(new.strata)
n.Ustrata <- length(Ustrata)
for(iStrata in 1:n.Ustrata){ # iStrata <- 1
indexStrata.object <- which(object.strata==Ustrata[iStrata])
indexStrata.newdata <- which(new.strata==Ustrata[iStrata])
all.times <- object.modelFrame[indexStrata.object,.SD$stop]
all.times <- sort(unique(all.times[all.times <= max(times)]))
if(length(all.times)>0){
res.tempo <- predictCox(object,
newdata = newdata[indexStrata.newdata,],
times = all.times,
type = "hazard")
lastEventTime.tempo <- res.tempo$lastEventTime[which(infoVar$strata.levels==Ustrata[iStrata])]
if(0 %in% all.times){
hazard.tempo <- cbind(res.tempo$hazard,NA)
all.timesA <- c(all.times,lastEventTime.tempo + 1e-10)
}else{
hazard.tempo <- cbind(0,res.tempo$hazard,NA)
all.timesA <- c(0,all.times,lastEventTime.tempo + 1e-10)
}
if(original.res$diag){
index.jump <- prodlim::sindex(eval.times = times[indexStrata.newdata], jump.times = all.timesA)
original.res$survival[indexStrata.newdata,] <- cbind(sapply(1:length(indexStrata.newdata), function(iP){ # iP <- 16
prod(1-hazard.tempo[iP,1:index.jump[iP]])
}))
}else{
index.jump <- prodlim::sindex(eval.times = times, jump.times = all.timesA)
survival.PL <- t(apply(hazard.tempo, 1, function(x){cumprod(1-x)}))
original.res$survival[indexStrata.newdata,] <- survival.PL[,index.jump,drop=FALSE]
}
}
}
}else{
all.times <- unique(sort(object.modelFrame[stop <= max(times), stop]))
if(length(all.times)>0){
res.tempo <- predictCox(object,
newdata = newdata,
times = all.times,
type = "hazard")
if(0 %in% all.times){
hazard.tempo <- cbind(res.tempo$hazard,NA)
all.timesA <- c(all.times,res.tempo$lastEventTime+1e-10)
}else{
hazard.tempo <- cbind(0,res.tempo$hazard,NA)
all.timesA <- c(0,all.times,res.tempo$lastEventTime+1e-10)
}
index.jump <- prodlim::sindex(eval.times = times, jump.times = all.timesA)
if(original.res$diag){
original.res$survival <- cbind(sapply(1:NROW(newdata), function(iP){ # iP <- 16
prod(1-hazard.tempo[iP,1:index.jump[iP]])
}))
}else{
survival.PL <- t(apply(hazard.tempo, 1, function(x){cumprod(1-x)}))
original.res$survival <- survival.PL[,index.jump,drop=FALSE]
}
}
}
## ** export
return(original.res)
}
#----------------------------------------------------------------------
### predictCoxPL.R ends here
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.