R/sfn.R

In scientific-computing-solutions/eventPrediction: Event Prediction in Clinical Trials with Time-to-Event Outcomes

#Code for the predict from parameters Survival functions/
#functions to be used in the integral when calculating the
#expected event times and times at risk


##' A Class containing the Survival Function (for a single arm) 
##' used in the integral to calculate event times
##' in predict from parameters. It should not be created by the user
##' but is created for the user (see sfns slot of \code{AnalysisResult})
##' 
##' @slot sfn A function sfn(x) is the function to be included in the events.integrate 
##' procedure when using predict from parameters. 
##' Specifically sfn(x) =  1 - P(had event [i.e. not dropout] by time x) and when no drop outs 
##' this is exactly the survival function.
##' If using dropouts this is not quite the survival function
##' However, the output of LatexSurvivalFn(object) is the survival function. 
##' @slot SurvivalFunction The actual survival function (if drop outs/finite followup are used then this will
##' not equal sfn). If finite follow up is used then S(x) = 0 for x > followup
##' @slot pdf The pdf function associated with the survival function. If finite follow up is used
##' then pdf(x) should  = 0 for x > followup 
##' @slot nullf Logical, TRUE if the object represents NULL (i.e. a survival function for 
##' a second arm in a single arm study)
##' @slot lambda The rate parameters for the arm of trial.
##' In a trial with lag, this is the rate parameters for time > T
##' @slot lambdaot If a lag was used then the rate parameters 
##' for time < T otherwise NA
##' @slot shape The Weibull shape parameter 
##' @slot LagT The lagtime for the survival function (0 for no lag) 
##' @slot followup The follow up time for each subject 
##' (Inf for studies with no fixed followup)
##' @slot dropout.shape The Weibull shape parameter for the drop out hazard function 
##' @slot dropout.lambda The rate parameter for the drop out hazard function = 0 if no dropout 
##' 
##' @export  
setClass("Sfn", 
          slots=list(sfn="function",
                     SurvivalFunction="function",
                     pdf="function",
                     nullf="logical",
                     lambda="numeric",
                     lambdaot="numeric",
                     shape="numeric",
                     LagT="numeric",
                     followup="numeric",
                     dropout.shape="numeric",
                     dropout.lambda="numeric"
                )  
)

# Constructor for \code{Sfn} object
# @param lambda The rate parameters for the arm of the trial.
# In a trial with lag, this is the rate parameters for time > T
# @param lambdaot If a lag was used then the rate parameters 
# for time < T otherwise NA
# @param lag.T The lagtime for the survival function (0 for no lag) 
# @param shape The Weibull shape parameter 
# @param followup The follow up time for each subject (Inf for studies with no fixed followup)
# @param dropout.shape The Weibull shape parameter for the drop out hazard function 
# @param dropout.lambda The rate parameter for the drop out hazard function - 0
# if no drop out
Sfn <- function(lambda,lambdaot,lag.T,shape,followup,dropout.shape=1,dropout.lambda=0){
  
  if(lag.T==0){
    lambdaot <- as.numeric(NA)
  }
  
  #survival function without dropouts or fixed follow up...
  f <- if(lag.T==0) 
         function(x){
           return(exp(-(lambda*x)^shape))
         }
       else 
         function(x){
           ifelse(x<lag.T, exp(-(lambdaot*x)^shape),exp(-(lambda*x)^shape+(lambda^shape-lambdaot^shape)*lag.T^shape))
         }
  
  #... and its corresponding hazard function 
  h <- if(lag.T==0) 
         function(x){
           return(shape*lambda^shape*x^(shape-1))
         }
       else 
         function(x){
           ifelse(x<lag.T,shape*lambdaot^shape*x^(shape-1),shape*lambda^shape*x^(shape-1))
         }
  
  ################################################################
  
  #next we calulate the function giving 1- P(having event) which is placed in the sfn slot 
  #of the object when including dropouts this function is the integrated or cumulative hazard for hazard
  #event. 
  
  #The integrand of the cumulative hazard function
  f2 <- function(x){h(x)*f(x)*exp(-(x*dropout.lambda)^dropout.shape)}
  
  #we precompute the value of the integral when there is dropout and lag
  lag.T.int <- if(lag.T != 0 && dropout.lambda != 0) integrate(f=f2,lower=0,upper=lag.T)$value else 0
  
  #we precompute the value of 1-P(having event by time=follow up) to save calculating the integral
  #multiple times
  follow.up.int <- if(dropout.lambda!=0 && !is.infinite(followup)) 
                      integrate(f=f2,lower=0,upper=followup)$value else 0 
  
  
  #include drop out and follow up if required for slot sfn (so no longer the survival function)
  dropf <- function(x){
    #First deal with case of no dropout
    if(dropout.lambda==0){
      return(ifelse(x<followup,f(x),f(followup)))
    }
        
    #next no lag but followup 
    if(lag.T==0 && !is.infinite(followup)){
      return(sapply(x, function(t){
        ifelse(t < followup,
          1-integrate(f=f2,lower=0,upper=t)$value,
          1-follow.up.int) #1-P(have event by time t > followup) = 1-P(have event by time followup)
      }))
    }
    
    #Next no lag or followup < lag (so no discontinuity)
    if(lag.T==0 || followup <= lag.T){    
      return(sapply(x,function(t){1-integrate(f=f2,lower=0,upper=min(t,followup))$value}))
    } 
    
    #finally case with lag and follow up > lag   
    sapply(x, function(t){
      ifelse(t < lag.T,
          1-integrate(f=f2,lower=0,upper=min(t,followup))$value,  
          1-integrate(f=f2,lower=lag.T,upper=min(t,followup))$value-lag.T.int #split integral into two parts   
      )
    })
      
  }
  #the function to be included in the events integral
  sfn <- dropf

  ########################################
  
  #Now we deal with the Survival function
  #Add drop outs
  S <- function(x){
    if(dropout.lambda==0){
      return(f(x))
    }
    f(x)*exp(-(x*dropout.lambda)^dropout.shape)
  }
  
  #and fixed follow ups
  #In '<=' the equal sign here is crucial as we need 
  #to use S(followup) when calculating the atrisk information
  SurvivalFunction <- function(x){ifelse(x<=followup,S(x),0)} 
  
  #######################################
  
  ##Next work on pdf
    
  dropouth <- function(x){
    if(dropout.lambda==0){
      return(0)
    }
    return((dropout.lambda*dropout.shape)*(dropout.lambda*x)^(dropout.shape-1))
  }
  
  #pdf is hazard function * survival function
  pdfwithdropout <- if(dropout.lambda==0) function(x){S(x)*h(x)} else function(x){S(x)*(h(x)+dropouth(x))}
      
  ##and fixed follow up
  pdf <- if(is.infinite(followup)) pdfwithdropout else function(x){ifelse(x<=followup,pdfwithdropout(x),0)}
  
  new("Sfn",sfn=sfn,nullf=FALSE,
      lambda=lambda,lambdaot=lambdaot,LagT=lag.T,shape=shape,followup=followup,
      dropout.lambda=dropout.lambda,dropout.shape=dropout.shape,SurvivalFunction=SurvivalFunction,pdf=pdf)
  
}

# Create a Null version of the \code{Sfn} object 
# For use as the survival function for the experimental arm
# of a single arm trial
# @return A \code{Sfn} object
NullSfn <- function(){
  f <- function(x){0}
  new("Sfn",sfn=f,nullf=TRUE,lambda=0,lambdaot=0,LagT=0,shape=0,followup=0,
      dropout.shape=1,dropout.lambda=0,SurvivalFunction=f,pdf=f)
}

##' Method to output a Latex String of the object
##' 
##' @param results Object to output 
##' @param \ldots Additional parameters to be passed into the function
##' @param lambda The symbol to be used for the arm's rate 
##' parameter, if lagged study then this is a vector of before and after lag
##' symbols. The latex backslash character needs to be escaped  
##' @param shape The symbol to be used for Weibull shape parameter
##' @return A latex string of the object (for \code{Sfn} this is the Survival function)
##' @rdname LatexSurvivalFn-methods
##' @name LatexSurvivalFn
##' @export
setGeneric( "LatexSurvivalFn", 
            def = function( results, ... )
              standardGeneric( "LatexSurvivalFn" ))


##' @name LatexSurvivalFn
##' @rdname LatexSurvivalFn-methods
##' @aliases LatexSurvivalFn,Sfn-method
##' @export
setMethod("LatexSurvivalFn",
  signature("Sfn"),
  function(results,decimalplaces,lambda,shape){
    
    #First set the function which will produce the output 
    if(results@LagT!= 0 && results@shape ==1 ){
      LatexOutput <- function(lambda,shape,followup,dropouttext){
        return(paste("$$ S(t) = \\begin{cases} \\exp(-",lambda[1],"t",dropouttext,") & \\quad \\text{if } t \\le T \\\\
                     \\exp(-",lambda[2],"t-(",lambda[1],"-",lambda[2],")T",dropouttext,")  & \\quad \\text{if }  t > T\\ \\end{cases} $$"))
      }
    }else if(results@LagT!= 0 && results@shape !=1){
      LatexOutput <- function(lambda,shape,followup,dropouttext){
        return(paste("$$ S(t) = \\begin{cases} \\exp(-(",lambda[1],"t)^",shape,dropouttext,") & \\quad \\text{if } t \\le T \\\\
                     \\exp(-(",lambda[2],"t)^",shape,"-(",lambda[1],"^",shape,
                     "-",lambda[2],"^",shape,")T^",shape,dropouttext,")  & \\quad \\text{if }  t > T\\ \\end{cases} $$"))
      }
    }else if(results@LagT == 0 && results@shape ==1){  
      LatexOutput <- function(lambda,shape,followup,dropouttext){
          fuptxt <- if(is.infinite(followup)) "" else paste("\\quad (\\text{if } t <",round(followup,decimalplaces),")") 
          return(paste("$$S(t) = \\exp(-",lambda[1],"t",dropouttext,")",fuptxt,"$$"))
      }
    }else{    
      LatexOutput <- function(lambda,shape,followup,dropouttext){
        fuptxt <- if(is.infinite(followup)) "" else paste("\\quad (\\text{if } t <",round(followup,decimalplaces),")") 
        return(paste("$$S(t) = \\exp(-(",lambda[1],"t)^",shape,dropouttext,")",fuptxt,"$$"))
      }
    }
    
    #Then produce the output
    ans <- ""
    extra <- ""
    
    dropouttext <- ""
    if(results@dropout.lambda!=0){
      
      if(results@dropout.shape==1){
        dropouttext <- paste("-",lambda[3],"t")
      }
      else{
        dropouttext <- paste("-(",lambda[3],"t)^",shape[2])
      }
    }
    
    
    if(results@LagT!=0){
      extra <- paste(",\\:",lambda[2],"=",round(results@lambda,decimalplaces))
      lambda_val <- round(results@lambdaot,decimalplaces)
    }
    else{
      lambda_val <- round(results@lambda,decimalplaces)
    }    
    
    ans <- paste(ans,LatexOutput(lambda,shape=shape[1],followup=results@followup,dropouttext),"\n")
    ans <- paste(ans,"$$",lambda[1],"=",lambda_val,extra)
    
    if(results@shape !=1) {
      ans <- paste(ans,",\\:",shape[1],"=",round(results@shape,decimalplaces))
    }
    
    if(results@dropout.lambda!=0){
      ans <- paste(ans,",\\:",lambda[3],"=",round(results@dropout.lambda,decimalplaces))
      
      if(results@dropout.shape!=1){
        ans <- paste(ans,",\\:",shape[2],"=",round(results@dropout.shape,decimalplaces))
      }
    }
    
    ans <- paste(ans,"$$\n") 
    return(ans)
  }
)

# Calculate the expected time at risk
# 
# @param SurvFn sfn object. If nullf slot of sfn object is true 
# then 0 time at risk occurs occur (this is used
# when calculating single arm trials) 
# @param B Accrual time.
# @param k Non-uniformity of accrual (numeric, 1=uniform).
# @param t A vector of prediction times for which the expected time  
# at risk is required
# @return A vector of time at risk at the given times
atrisk.integ <-  function( SurvFn, B, k, t ) {
  if(SurvFn@nullf){
    return (rep(0,length(t)))
  }
  
  #see the predict from parameters vignette for the three formula calculated here
  #and their derivations
  
  #####################
  
  #First calculate the amount of time of subjects who are still on the trial at time t
  #= integral from 0 to min(t,B) of P(still on trial at t | rec at s)P(rec at s)*time_at_risk ds
  #and time_at_risk=t-s
  myf1 <- function(s,k,t){(k*s^(k-1)/B^k)*SurvFn@SurvivalFunction(t-s)*(t-s)}
  
  i1 <- unlist(lapply(t,
                      function(x){
                        if(x==0) 0 else integrate(f=myf1,lower=0,upper=min(x,B),k=k,t=x)$value
                      })) 
  
  #Next calculate the amount of time of subjects who had an event or dropped out before time t
  #This does not include subjects who have dropped out as they were censored at the end of the follow up period
  
  #This function is given by 
  #integral from 0 to t [ W ] dt' where W = time at risk for subject who had event at time t'
  #specifically
  # W = integral from 0 to t' P(had event at time t' | rec at s)P(rec at s)(t'-s)  ds
  #notw the "at" in P(had event at time t') implies we need the pdf rather than the survival function
  
  internalf <- function(s,t,k,B){sapply(s,function(x,t,k,B){(k*x^(k-1)/B^k)*SurvFn@pdf(t-x)*(t-x)},t=t,k=k,B=B)}
  myf2 <- function(t,k,B){sapply(t,function(x,k,B){
    upper <- min(x,B)
    lower <- max(x-SurvFn@followup,0)
    lagpoint <- max(x- SurvFn@LagT,0)
    if(lower > upper){
      return(0)
    }
    
    if(lagpoint >= upper || lagpoint <= lower){
      return( integrate(internalf,lower=lower,upper=upper,t=x,k=k,B=B)$value)
    }
    
    integrate(internalf,lower=lower,upper=lagpoint,t=x,k=k,B=B)$value+
      integrate(internalf,lower=lagpoint,upper=upper,t=x,k=k,B=B)$value},k=k,B=B)
    
    
  }
  
  i2 <- unlist(lapply(t,function(x){
    if(x==0) 0 else integrate(f=myf2,lower=0,upper=x,k=k,B=B)$value  
  }))
  
  #Finally calculate the time on study for subjects who dropped out at the follow up period time before time t
  #The three terms multiplied are: 
  #P(surviving until F)
  #Time spent at risk (i.e. F)
  #P(recruited before time t-F), 
  i3 <- unlist(lapply(t,function(x){
    if(x < SurvFn@followup) 0 else SurvFn@SurvivalFunction(SurvFn@followup)*SurvFn@followup*(min(B,x-SurvFn@followup)^k)/B^k
  }))
  
  #And sum the three values together
  i1+i2+i3
}

# Calculate expected events using adaptive quadrature
# 
# This function uses R's integrate function to calculate
# the integral s^{k-1}*SurvFn(t-s) ds for given
# k, t and SurvFn. The limits of the integral are 0 to 
# min(t,B)
# 
# @param SurvFn sfn object. If nullf slot of sfn object is true 
# then 0 events occur (this is used
# when calculating single arm trials) 
# @param B Accrual time.
# @param k Non-uniformity of accrual (numeric, 1=uniform).
# @param t A vector of prediction times for which the number 
# of events is required
# @return A vector of number of events at the given times
events.integ <- function( SurvFn, B, k, t ) {
  
  if(SurvFn@nullf){
    return (rep(0,length(t)))
  }
  
  #the function to be integrated (s is the dummy variable)
  #see vignette
  myf <- function(s,k,t){s^{k-1}*SurvFn@sfn(t-s)}
  
  
  integrated_val <- unlist(lapply(t,
           function(x){
             if(x==0){
               return(0) #integrate unhappy if both limits are 0 and k < 1
             }
             upperlim <- min(x,B)
             
             #If integrand is differentiable over limits then integrate in one go
             if(SurvFn@LagT == 0 || upperlim <= SurvFn@LagT  ){
               return(integrate(f=myf,lower=0,upper=upperlim,k=k,t=x)$value)
             }
                                    
             #when there is a lag the integrand is non-differentiable
             #so splitting the limits into two differentiable parts 
             #vastly improves performance
             split.point <- upperlim - SurvFn@LagT
             return(integrate(f=myf,lower=0,upper=split.point,k=k,t=x)$value+
                    integrate(f=myf,lower=split.point,upper=upperlim,k=k,t=x)$value)
                                    
             }       
  ))
  #See vignette for how integrated_val is used
  (pmin(t,B)/B)^k - k*integrated_val/B^k
}