R/accrual.R

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

#This file contains the generator functions for recruiting extra subjects 
#into a predict from data simulation. It also contains the code for
#estimating k in the accrual cdf G(t) t^k/B^k from recruitment data.

##' @include eventData.R fromDataInternal.R
NULL

##' A class to allow additional subjects 
##' to be recruited onto the study
##' @slot f A function with 1 argument (the number of subjects to be recruited)
##' which returns a vector of recruitment times (of the required length)
##' @slot model A string, name of the accrual procedure used  
##' @slot text The summary text to be output when printing the accrual procedure
##' for example "a poisson process with rate..."
##' @seealso \code{\link{show,AccrualGenerator-method}}
##' @export
setClass("AccrualGenerator",
         slots=list(f="function",
                    model="character",
                    text="character"
         )
)

##' @name show
##' @rdname show-methods
##' @aliases show,AccrualGenerator-method
##' @export
setMethod("show","AccrualGenerator",function(object){
  cat(paste("Accrual model:",object@model,"\n"))
  cat(paste("Model Description:\n"))
  cat(object@text)  
})

##' Create A Poisson process accrual generator
##' @param start.date The start date for the recruitment process
##' @param rate The rate of the required Poisson process.
##' @return An \code{AccrualGenerator} object which can 
##' be used when recruiting subjects with a Possion process
##' with given \code{start.date} and \code{rate}
##' @export 
Generate.PoissonAccrual <- function(start.date,rate){
  if(rate <= 0) stop("rate must be positive")
  start.date <- FixDates(start.date)
  f <- function(N){
    start.date + ceiling( cumsum( rexp( N, rate=rate ) ) )
  }
  
  r <- if(rate < 0.01) r else round(rate,2)
  
  text <- paste("a poisson process with rate=", 
                 r, " and start date=", start.date,".",sep="")
  
  new("AccrualGenerator",f=f,
      model="Poisson process",text=text)    
    
}

##' Create An AccrualGenerator using a power law recruitment
##' 
##' Subjects are accrued according to the c.d.f 
##' \code{G(t)=t^k/B^k} where \code{k} is a parameter, 
##' \code{t} is the time and \code{B} is the recruitment period. 
##'  
##' See the predict from data vignette for further details  
##' @param start.date The start of the subject accrual period
##' @param end.date The date the last subject is accrued so \code{B} = 
##' \code{end.date} - \code{start.date} unless \code{rec.start.date} is used see 
##' below
##' @param k The non-uniformity accrual parameter
##' @param deterministic Logical, if FALSE then the recruitment times 
##' are non-stochastically chosen so that their cumulative distribution function is \code{G(t)}
##' otherwise they are generated by sampling random variables with a cdf \code{G(t)}
##' @param rec.start.date If this argument is used the subjects are still recruited between 
##' start.date and end.date but they follow the cdf \code{G(t)=(t^k-L^k)/(B^k-L^k)} where
##' \code{t} is in \code{[L,B]} and \code{B = end.date - rec.start.date} and \code{L = start.date- rec.start.date}.
##' @return An AccrualGenerator object which generates the required subject accruals
##' @export
Generate.Accrual <- function(start.date,end.date,k,deterministic=FALSE,rec.start.date=NULL){
  start.date <- FixDates(start.date)
  end.date <- FixDates(end.date)
  rec.start.date <- if(is.null(rec.start.date)) start.date else FixDates(rec.start.date)
  
  # Length of recruitment period
  B <- as.integer(end.date - rec.start.date)
  L <- as.integer(start.date - rec.start.date)
  
  if(L < 0 || k <= 0|| B <= L){
    stop("Invalid arguments")
  }
  
  if(deterministic){
    f <- function(N){
      as.Date(((1:N/N)*(B^k-L^k) + L^k)^(1/k)-L,origin=start.date)
    }
    stext <- "non-stochastic"
  }
  else{
    f <- function(N){
      ans <- ((B^k-L^k)*runif(N) + L^k)^(1/k)-L
      ans[order(ans)]
      as.Date(ans,origin=start.date)
    }
    stext <- "stochastic"
  }  
  
  text <- paste("a ",stext  ," allocation following G(t)=t^k/B^k, where k=", 
  round(k,2), " and B=",B," days is the recruitment period [", rec.start.date, ", ", end.date, "].",sep="")
  
  if(L!=0){
    text <- paste(text," New subjects will be recruited after ",start.date,".",sep="") 
  }
  
  new("AccrualGenerator",f=f,
      model="Power law allocation",text=text) 
  
} 

##' Generic function to estimate non-uniformity accrual 
##' parameter \code{k}
##' 
##' Note The estimate produced assumes recruitment has completed.
##' @docType methods
##' @param object An \code{EventData} object
##' @param ... Additional arguments to be passed to specific method
##' @rdname estimateAccrualParameter-methods  
##' @export
setGeneric("estimateAccrualParameter",function(object,...) standardGeneric("estimateAccrualParameter"))



##' @param start.date The startdate for recruitment, by default the date 
##' the first subject is recruited. If subjects are recruited on start.date, it is
##' assumed they are recruited on the following day to avoid log(0) in the 
##' maximum likelihood estimate.
##' @param end.date The end date for recruitment, by default the date 
##' the last subject is recruited. 
##' @return A maximum likelihood estimate for k, see vignette 
##' for further details    
##' @rdname estimateAccrualParameter-methods 
##' @aliases estimateAccrualParameter,EventData-method
##' @export
setMethod("estimateAccrualParameter", "EventData",
  function(object, start.date=min(object@subject.data$rand.date),
           end.date=max(object@subject.data$rand.date)){
    message("Note this estimate assumes recruitment has completed")          
    
    start.date <- FixDates(start.date)
    end.date <- FixDates(end.date)
    
    B <- as.numeric(end.date-start.date)
    if(B <= 0) stop("end.date must be after start.date")
    
    recs <- (object@subject.data$rand.date-start.date)
    recs <- ifelse(recs==0,0.5,recs)
    
    if(any(recs< 0 | recs > B)) stop("Subjects recruited outside of recruitment period")
    
    return(MLestimateK(B,as.numeric(recs)))
    
  }
)