R/eventData.R

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

# This file contains the public functions associated with eventData object
# in the predict from data part of the package. Note the estimateAccrualParameter
# function can be found in the accrual.R file.

##' @include fromDataSimParam.R eventPrediction_package.R
NULL

##' @import survival methods
NULL

##' @importFrom  mvtnorm rmvnorm
NULL

##' @importFrom  scales date_format
NULL

# Function that will check the validity of EventData object
# when constructed
# @param object EventData object
checkEventData <- function(object){
  errors <- ""
  
  required.colnames <- c("has.event", "rand.date", "withdrawn", "subject", "time","censored.at.follow.up")
  
  data <- object@subject.data
  
  if(nrow(data)>0){
  
    if(any(!data$has.event %in% c(0,1))) errors <- paste(errors,"has.event must be 1 or 0" )
    if(any(!data$withdrawn %in% c(0,1))) errors <- paste(errors,"withdrawn must be 1 or 0" )
    if(any(!data$censored.at.follow.up %in% c(0,1))) errors <- paste(errors,"censored.at.follow.up must be 1 or 0" )

    if(any(!required.colnames %in% colnames(data))) errors <- paste(errors,"invalid column names")
  
    
    if(class(data$rand.date)!="Date") errors <- paste(errors,"Invalid rand.date")
    if(class(data$time)!="numeric") {
      errors <- paste(errors,"Times must be numeric")
    }
    else{
      if(all(data$time==0)) errors <- paste(errors,"time value incorrect all are 0!")  
      if(any(data$time < 0 | is.na(data$time))) {
        errors <- paste(errors,"subjects cannot have non-positive time on trial. Please check subjects ",
                      paste(data[data$time <= 0 | is.na(data$time),]$subject,collapse=", "))
      
      }
    }
  
    if(any(duplicated(data$subject))) errors <- paste(errors,"subject ID must be unique")
  
    if(all(data$has.event==0)){
      warning("No events have occurred - a model cannot be fit")
    }
    if(all(data$has.event==1)) errors <- paste(errors,"all events have occurred!")
  
    idx <- which(data$has.event==1& data$withdrawn==1)
    if(length(idx)){
      errors <- paste(errors,"subjects cannot be both withdrawn and have an event. Please check subjects ",
                    paste(data$subject[idx],collapse=", "))
    }
  
  }
  
  if(!is.numeric(object@followup)||length(object@followup)!= 1|| object@followup <= 0){
    errors <- paste(errors,"Invalid followup argument")
  }
  
  if(errors == "") TRUE else errors
}
  
##' Class representing data to use for new predictions
##' @slot subject.data a data frame with 6 columns
##' "subject", "rand.date", "has.event", "withdrawn", "censored.at.follow.up" and "time" and "site" and "event.type"
##' see vignette and \code{EventData} for further details. 
##' @slot followup A numeric value for the fixed follow up period a subject
##' is followed. This is in days. If there is no fixed followup period then 
##' Inf should be used
##' @export
setClass("EventData", 
          slots= list(subject.data = "data.frame",
                      followup="numeric"),
          validity = checkEventData)


##' Constructor for event data object
##' 
##' All dates must be in one of the following formats:
##' YYYY-MM-DD, DD/MM/YY or DD Month YYYY
##'  
##' 
##' @param data A data frame 
##' @param subject string, the column name of subject identifiers
##' @param rand.date string, the column name of randomization dates
##' @param has.event string, the column name of whether a subject had an event
##' (1) or not (0) 
##' @param withdrawn string, the column name of whether a subject has withdrawn
##' (1) or not (0)  
##' @param time Either a string, the column name of time each subject has been on the study
##' or a list with elements with (some of) the following named elements.
##' \code{last.date}, \code{event.date}, \code{prog.date}, \code{dth.date} and \code{withdrawn.date}
##' In this case the package will attempt to derive the time column 
##' See the vignette and then eventPrediction:::AddTimeColumn for further details
##' @param site optional column for subject site
##' @param event.type optional column for the event type (e.g. unstable angina) if not included then `Had Event' will be used 
##' @param remove.0.time logical, if TRUE then all subjects with time = NA or 0 are removed from the
##' data set and not included in the object. If FALSE then they are included in the simulation (but not in the model fitting)
##' @param followup A numeric value for the fixed follow up period a subject
##' is followed. This is in days. If there is no fixed followup period then 
##' Inf should be used
##' @return An \code{EventData} object
##' @export
EventData <- function(data,subject,rand.date,has.event,withdrawn,time,site=NULL,event.type=NULL,remove.0.time=FALSE,followup=Inf){
  
  #set time argument
  arg <- if(!is.list(time))  time else NULL
   
  #check columns are in data frame 
  for(x in c(subject,rand.date,has.event,withdrawn,arg)){
    if(!is.null(x) && !x %in% colnames(data)){
      stop(paste("Column name",x,"not found in data frame"))
    }
  }
  
  if(!is.numeric(followup) || length(followup)> 1 || followup <= 0){
    stop("Invalid followup argument")
  }
  
  #validate the input columns, deriving time from other columns if needed
  data[,rand.date] <- if(nrow(data)>0) FixDates(data[,rand.date]) else data[,rand.date]
  time <- if(!is.list(time))  data[,time] else AddTimeColumn(data,rand.date,has.event,withdrawn,subject,time)
  site <- if(is.null(site)) rep(NA,nrow(data)) else data[,site] 
  
  if(is.null(event.type)) 
    event.type <- ifelse(data[,has.event],"Has Event",factor(NA))
  else{
    event.type <- factor(data[,event.type])
    if("" %in% levels(event.type)){
      event.type[event.type==""] <- factor(NA)
      event.type <- droplevels(event.type)
    }
  }
  
  #unvalidated data frame
  subject.data <- data.frame(
    subject=data[,subject],
    rand.date=data[,rand.date],
    time=time,
    has.event=data[,has.event],
    withdrawn=data[,withdrawn],
    site=site,
    event.type=event.type
  )
  
     
  keep <- which(as.character(subject.data$subject) > "" & !is.na(subject.data$rand.date))
  if(length(keep) != nrow(subject.data)){
    warning(paste("Subjects",paste(subject.data[setdiff(1:nrow(subject.data),keep),subject],collapse=", "),
                  "have been removed due to invalid rand.date or missing subject ID."))
  }
  
  subject.data <- subject.data[keep,]
  
  subject.data <- warnNAtime(remove.0.time,subject.data)
  
  subject.data <- IncludeFollowUp(subject.data,followup)
  
  #finally run some validation checks which change the input data and output warnings
  validation.checks <- function(subject.data,idx,warning.message,fn){
    if(length(idx)>0){
      warning(paste("Subjects",paste(subject.data[idx,subject],collapse=", "),warning.message))
      subject.data <- fn(subject.data,idx)
    }
    subject.data
  }
  
  subject.data <- validation.checks(subject.data,
                    idx=which(subject.data$has.event==1 & is.na(subject.data$event.type)),
                    warning.message="have an event and no event type, their event type is set to 'Had Event'",
                    fn=function(subject.data,idx){
                      if(!"Had Event" %in% levels(subject.data$event.type)){
                        levels(subject.data$event.type) <- c(levels(subject.data$event.type),"Had Event")
                      }
                      subject.data$event.type[idx] <- "Had Event"
                      return(subject.data)
                    })
  
  subject.data <- validation.checks(subject.data,
                    idx=which(subject.data$has.event!=1 & !is.na(subject.data$event.type)),
                    warning.message="have an event type but did not have an event",
                    fn=function(subject.data,idx){
                      subject.data$event.type[idx] <- factor(NA)
                      subject.data$event.type <- droplevels(subject.data$event.type)
                      return(subject.data)
                    })
  
  subject.data <- validation.checks(subject.data,
                    idx=which(subject.data$has.event==1& subject.data$withdrawn==1),
                    warning.message="have an event and are withdrawn, assuming they have an event at the given time",
                    fn=function(subject.data,idx){
                      subject.data$withdrawn[idx] <- 0
                      return(subject.data)
                    })

  #further validation occurs in the validity function of the class 
  return(new("EventData", subject.data = subject.data, followup=followup))
}



##' @name show
##' @rdname show-methods
##' @aliases show,EventData-method
##' @export
setMethod("show",
          "EventData",
    function(object) {
         cat("EventData, use object@subject.data$param to access individual columns:\n")
         cat(str(object@subject.data))
     })


##' @name summary
##' @rdname summary-methods
##' @aliases summary,EventData-method
##' @export
setMethod("summary","EventData",
  function(object){
    df <- object@subject.data 
    daysinyear <- standarddaysinyear()
    
    if(nrow(df)>0){
    
      cat(paste("Number of subjects:",nrow(df),"\n"))
      cat(paste("Number of events:",sum(df$has.event==1),"\n"))
    
      if(sum(df$has.event==1) != 0 && length(levels(df$event.type)) != 1){
        tab <- table(df$event.type)
        lapply(names(tab),function(y){cat(paste(" Of type ",y,": ",tab[y],"\n",sep=""))})
      }
    
      cat(paste("Number of withdrawn:", sum(df$withdrawn==1),"\n"))
      cat(paste("First subject randomized:",as.Date(min(df$rand.date),origin="1970-01-01"),"\n"))
      cat(paste("Last subject randomized:",as.Date(max(df$rand.date),origin="1970-01-01"),"\n"))
    
      if(sum(object@subject.data$has.event)>0){    
        cat(paste("First Event:",as.Date(min(LastDate(df[df$has.event==1,]),na.rm=TRUE),origin="1970-01-01"),"\n"))
        cat(paste("Last Event:",as.Date(max(LastDate(df[df$has.event==1,]),na.rm=TRUE),origin="1970-01-01"),"\n"))
      }
      
      av <- average.rec(nrow(df),as.numeric(min(df$rand.date)),as.numeric(max(df$rand.date))) 
      cat(paste("Average recruitment (subjects/day):",roundForceOutputZeros(av,2),"\n"))
    
      if(!is.infinite(object@followup)){
        cat(paste("Subjects followed for ",round(object@followup), " days (",
                round(object@followup/daysinyear,2)," years)\n",sep=""))
      
        cat(paste("Number of subjects censored at end of follow up period:",
                  sum(object@subject.data$censored.at.follow.up),"\n"))
      
      }
    }
    else{
      cat("Empty data frame!\n")
    }
})



##' Fit a survival model to an EventData object
##' @rdname fit-methods
##' @name fit
##' @param object an EventData object
##' @param ... Additional arguments to the function
##' @param dist character, either weibull or loglogistic, which model
##' is to be fit 
##' @return an EventModel object
##' @export
if (!isGeneric("fit")){
  setGeneric("fit", function(object,...) standardGeneric("fit"))
}



##' @name fit
##' @rdname fit-methods
##' @aliases fit,EventData-method
##' @export
setMethod("fit","EventData",function(object,dist="weibull"){
  
  if(!dist %in% c("weibull","loglogistic")){
    stop("dist must be weibull or loglogistic")
  }
  
  if(nrow(object@subject.data)==0)stop("Empty data frame!") 
  if(sum(object@subject.data$has.event)==0){
    stop("Cannot fit a model to a dataset with no events")
  }
  
  #subjects with time = 0 are set to NA for the model fitting
  #so they are ignored
  indat <- object@subject.data
  indat$time <- ifelse(indat$time==0,NA,indat$time)
  
  model<-survreg(Surv(time, has.event) ~ 1, data=indat, dist=dist, y = TRUE)
  
  new("EventModel",model=model,event.data=object,
      simParams=FromDataParam(object=model,type=dist))
})



##' @rdname plot-methods
##' @name plot
##' @aliases plot,EventData,missing-method
##' @export
setMethod( "plot",
  signature( x="EventData", y="missing" ),
  function(x, xlab="log(t)", ylab="log(-log(S(t)))", main="", ...) {
    
    if(nrow(x@subject.data)==0)stop("Empty data frame!")
    if(sum(x@subject.data$has.event)==0){
      stop("Cannot fit a model to a dataset with no events")
    }
    
    model <- survfit(Surv(time, has.event) ~ 1, data=x@subject.data,...)
      
    res <- data.frame(t=model$time, s=model$surv)
    res <- res[res$t>0 & res$s>0 & res$s<1,]
    
    df <- data.frame(x=log(res$t),y=log(-log(res$s)))
    
    p <- ggplot(df, aes_string(x="x", y="y")) + geom_point() +
      stat_smooth(method="lm", se=FALSE, color="red", size=1 ) +
      xlab(xlab) + ylab(ylab) +
      ggtitle(main) + theme_bw() +
      theme(panel.grid.minor = element_line(colour="gray", size=0.5))
    p
  }
)



##' Method to create new \code{EventData} object only including recruitment times
##' @param object The  object from which to create the new \code{EventData} object
##' @rdname OnlyUseRecTimes-methods
##' @return An \code{EventData} object with time=0 and no events or withdrawn subjects
##' @name OnlyUseRecTimes
##' @export
setGeneric("OnlyUseRecTimes",function(object) standardGeneric("OnlyUseRecTimes"))

##' @rdname OnlyUseRecTimes-methods
##' @name OnlyUseRecTimes
##' @aliases OnlyUseRecTimes,EventData-method
##' @export
setMethod("OnlyUseRecTimes","EventData",function(object){
  object@subject.data$withdrawn <- 0
  object@subject.data$has.event <- 0
  object@subject.data$censored.at.follow.up <- 0
  object@subject.data$time <- 0
  object@subject.data$event.type <- factor(NA)
  object@subject.data$event.type <- droplevels(object@subject.data$event.type)
  object
})


##' Cut the Event Data
##' 
##' Create a new EventData from an existing one which shows how the
##' data would have looked on a given date (assuming no lag in reporting of 
##' events and if necessary subjects have been censored at cutoff point).
##' If subject was censored before the cut date then their censored date remains
##' unchanged
##' @param object The EventData object
##' @param date The date to cut the data at
##' @rdname CutData-methods
##' @name CutData
##' @return \code{EventData} object with the data censored at the appropriate point
##' @export
setGeneric("CutData",function(object,date) standardGeneric("CutData"))


##' @rdname CutData-methods
##' @name CutData
##' @aliases CutData,EventData-method
##' @export
setMethod("CutData","EventData",function(object,date){

  cut.date <- FixDates(date)
  
  subject.data <- object@subject.data
  
  #only include subjects who have been randomized by this time
  subject.data <- subject.data[subject.data$rand.date<=cut.date,]
  if(nrow(subject.data)==0){
    stop("Cut date before first subject randomization")
  }
  
  censored.times <- cut.date - subject.data$rand.date + 1
  idx <- which(censored.times < subject.data$time)
  
  if(length(idx)==0&&nrow(subject.data)==nrow(object@subject.data)){
    warning("Cut date is too late (no ongoing subjects at this time) and so has no effect")
  }
  
  
  
  subject.data$time <- pmin(censored.times,subject.data$time)
  subject.data$has.event[idx] <- 0
  subject.data$event.type[idx] <- factor(NA)
  subject.data$event.type <- droplevels(subject.data$event.type)
  subject.data$withdrawn[idx] <- 0
  
  EventData(
    data=subject.data,
    subject="subject",
    rand.date="rand.date",
    has.event="has.event",
    withdrawn="withdrawn",
    time="time",
    site="site",
    event.type="event.type",
    followup=object@followup
  )
  
})

##' This function calculates the event type for pfs data
##' 
##' This function is called before creating the EventData object 
##' @param has.event A vector or whether subjects had event 
##' @param prog.date A vector of dates (or strings in the standard eventPrediction allowed formats)
##' of progression.date or blank if unknown 
##' @param dth.date A vector of dates (or strings in the standard eventPrediction allowed formats)
##' of death date or blank if unknown
##' @return A vector with NA if subject does not have event, 
##' "Death" if dth.date <= prog.date or prog.date blank,
##' "Progression (not death)" if prog.date < dth.date or dth.date blank
##' "Progression (unknown if death)" if both prog.date and dth.date are blank
##' @export
CalculateProgEventTypes <- function(has.event,prog.date,dth.date){
  #called before creating EventDate object so extra validation needed  
  
  if(any(!has.event %in% c(0,1))){
    stop("Invalid has.event argument")
  }
  
  prog.date <- FixDates(prog.date)
  dth.date <- FixDates(dth.date)
  
  if(!(length(prog.date)==length(dth.date) && length(prog.date) ==length(has.event))){
    stop("Invalid lengths of arguments")
  }
  
  unlist(lapply(seq_along(has.event),function(x){
    if(!has.event[x]) return(NA)
    
    dth <- dth.date[x]
    prog <- prog.date[x]
    
    if(!is.na(prog) && !is.na(dth)){
      if(dth<=prog) return("Death")  
      return("Progression (not death)")
    }  
    
    
    #just dth return dth
    if(!is.na(dth)) return("Death")
    
    #just prog return prog
    if(!is.na(prog)) return("Progression (not death)")
    
    #Otherwise are relying on lastDate so do not know which
    return("Progression (unknown if death)")
    
  }))
  
  
}


##' Create an \code{EventData} object with no subjects
##' 
##' @param followup The time subjects are followed until censoring  
##' @export
EmptyEventData <- function(followup=Inf){
  
  if(length(followup) > 1 || !is.numeric(followup) || followup < 0){
    stop("Invalid followup argument")
  }
  
  d <- data.frame(subject=character(0),
                  randDate=numeric(0),
                  has.event=numeric(0),
                  withdrawn=numeric(0),
                  time=numeric(0))
  
  #convert to an empty EventData object
   EventData(data=d,
             subject="subject",
             rand.date="randDate",
             has.event="has.event",
             withdrawn="withdrawn",
             time="time",
             followup=followup)
}

##' Calculate the number of days at risk
##' 
##' For \code{EventData} object this is the number of days
##' at risk of the input data. For \code{SingleSimDetails} object
##' this is the median numberof days at risk in the simulation set at the given time
##' 
##' @param object The object to calculate the days at risk  
##' @param ... Additional arguments passed to the function
##' @param times A vector of dates for calculating cutting the simulated date at in order to
##' calculate the number of days at risk by this point.
##' @name CalculateDaysAtRisk
##' @rdname CalculateDaysAtRisk-methods
##' @export
setGeneric("CalculateDaysAtRisk",function(object,...)standardGeneric("CalculateDaysAtRisk"))


##' @rdname CalculateDaysAtRisk-methods
##' @name CalculateDaysAtRisk
##' @aliases CalculateDaysAtRisk,EventData-method
##' @export
setMethod("CalculateDaysAtRisk","EventData",
          function(object){
            return(sum(object@subject.data$time))            
          }
)

##' Makes a barplot with the number of events over time
##' 
##' @param object The object to create the event versus time plot 
##' @param ... Additional arguments passed to the function
##' @param timeunit The resoloution on the x-axis (month, weeks or quarter)
##' @name EventVsTimePlot
##' @rdname EventVsTimePlot-methods
##' @export
setGeneric("EventVsTimePlot",function(object,timeunit="Months",...)standardGeneric("EventVsTimePlot"))


##' @rdname EventVsTimePlot-methods
##' @name EventVsTimePlot
##' @aliases EventVsTimePlot,EventData-method
##' @export
setMethod("EventVsTimePlot","EventData",
          function( object, timeunit="month" ){
            my.data <- object@subject.data
            my.data$lastdate <- LastDate( my.data )
            my.data$has.event <- as.integer( my.data$has.event )
            my.data$mytimeunit <- as.Date( cut( my.data$lastdate, 
                                                breaks = timeunit ) )
            
            mybreaks <- if( timeunit == "quarter" ) "3 months" else "1 month" 
            
            ggplot( data = my.data, aes_string( x="mytimeunit", y="has.event" )) +
              stat_summary( fun.y = sum, geom = "bar") + 
              scale_x_date(
                date_labels = "%Y-%b",
                date_breaks = mybreaks ) + 
              theme_bw() + 
              theme( axis.text.x = element_text(angle = 45, hjust = 1) ) +
              xlab( "" ) +
              ylab( "Observed events" )
          }
)