Diff of /R/eventModel.R [000000] .. [f2e496]

Switch to side-by-side view

--- a
+++ b/R/eventModel.R
@@ -0,0 +1,116 @@
+# This file contains the public functions associated with
+# the fitted Weibull/loglogistic model for the predict from data simulations  
+
+##' @include eventData.R accrual.R fromDataSimInternal.R simQOutput.R ctrlSpec.R singleSimDetails.R fromDataSimParam.R
+NULL
+
+
+##' An S4 class containing a fitted survival model
+##' of an Eventdata object
+##' @slot model An S3 survreg object
+##' @slot event.data An EventData object used to fit a Weibull survial model
+##' @slot simParams The \code{FromDataSimParam} object which contains the information
+##' needed to generate the survial times
+##' @export
+setClass("EventModel", 
+         slots=list(model="survreg",
+                    event.data="EventData",
+                    simParams="FromDataSimParam")
+         )
+
+
+##' @name show
+##' @rdname show-methods
+##' @aliases show,EventModel-method
+##' @export
+setMethod("show",
+  "EventModel",
+  function(object) {
+    print(object@model)
+    
+})
+
+
+##' @param units Scale for the x-axis. "Days", "Months" or "Years"
+##' @name plot
+##' @rdname plot-methods
+##' @aliases plot,EventModel,missing-method
+##' @export
+setMethod( "plot",
+  signature( x="EventModel", y="missing" ),
+  function(x, units="Days", xlab=paste("Time in study [",units,"]",sep=""),
+           ylab="", main="", ylim=NULL, xlim=NULL, ...) { 
+    
+    xscale <- GetScaleForKM(units,daysinyear)
+    daysinyear <- standarddaysinyear()
+    
+    KM <- survfit(Surv(time, has.event) ~ 1, data=x@event.data@subject.data,...)
+    
+    plot(KM, lwd=c(2,1,1), col=c("red", "black", "black" ), 
+        xlab=xlab, ylab=ylab,
+        xlim=xlim, ylim=ylim,
+        main=main,xscale=xscale)
+    lines(predict(x@model, type="quantile", p=seq(.01,.99,by=.01))[1,]/xscale, 
+         seq(.99,.01,by=-.01), col="brown", type="l", lwd=3)
+    
+    pos <- if(is.null(xlim)) 0.75*max(KM$time/xscale) else xlim[1] + 0.75*(xlim[2]-xlim[1])
+    
+    legend(pos, 1, c( "Data", "Model" ), col=c( "red", "brown" ), lty=c(1,1))
+  
+    }
+)
+
+
+##' Calculate KM risk.table
+##' @rdname risk.table-methods
+##' @name risk.table
+##' @param x An \code{EventModel} object 
+##' @param ... Additional arguments to pass to survfit
+##' @return A data frame contianing the number at risk at the given times
+##' @export
+setGeneric("risk.table",
+           def=function(x,...)
+             standardGeneric("risk.table"))
+
+
+##' @rdname risk.table-methods
+##' @name risk.table
+##' @aliases risk.table,EventModel-method
+##' @param units Scale for the risk table: "Days", "Months" or "Years"
+##' @param risk.times The times for the risk table (in \code{units}). 
+##' If NULL then the 0, 0.25,0.5,0.75,1 quantiles of the KM survival table are used
+##' @export 
+setMethod("risk.table",
+  signature("EventModel"),
+  function(x,units="Days",risk.times=NULL,...){
+    daysinyear <- standarddaysinyear()
+    xscale <- GetScaleForKM(units,daysinyear)
+    KM <- survfit(Surv(time, has.event) ~ 1, data=x@event.data@subject.data,...)
+    if(is.null(risk.times)){
+      risk.times <- quantile(c(0, max(KM$time)))/xscale
+      names(risk.times) <- NULL
+    }  
+    
+    scaled.risk.times <-risk.times*xscale 
+    
+    
+    KM <- summary(KM,time=scaled.risk.times)
+    retVal <- data.frame(time = KM$time/xscale, n.risk = KM$n.risk)
+    
+      
+    n.risk <- unlist(lapply(risk.times,function(x){
+                 
+                  if(x %in% retVal$time){
+                    return(retVal$n.risk[retVal$time==x])
+                  }
+                  return("NA")
+    }))
+    
+    data <- rbind(risk.times,n.risk)
+    rownames(data) <- c("time","n.risk")
+    return(data)
+   
+  }
+)
+
+