R/RRPlot.R

In FRESA.CAD: Feature Selection Algorithms for Computer Aided Diagnosis

RRPlot <-function(riskData=NULL,
                  timetoEvent=NULL,
                  riskTimeInterval=NULL,
                  ExpectedPrevalence=NULL,
                  atRate=c(0.90,0.80),
                  atThr=NULL,
                  plotRR=TRUE,
                  title="",
                  ysurvlim=c(0,1.0))
{

  if (!requireNamespace("corrplot", quietly = TRUE)) {
    install.packages("corrplot", dependencies = TRUE)
  }

  OERatio <- NULL
  OE95ci <- NULL
  OAcum95ci <- NULL
  CumulativeOvs <- NULL
  DCA <- NULL
  OEData <- NULL
  OARatio <- NULL
  netBenefitatP <- NULL
  
  uvalues <- length(unique(riskData[,2]))
  isProbability <- (min(riskData[,2]) >= 0) && (max(riskData[,2]) <= 1.0) && (sd(riskData[,2]) > 0.00001)
  
  
  if (uvalues < 10)
  {
    warning(paste("Not a continous variable;",title))
    return (0)
  }
  ## Removing ties
  deltaR <- 0;
  if (uvalues < nrow(riskData))
  {
    deltaR <- IQR(riskData[,2])/nrow(riskData)
    if (deltaR == 0) 
    {
      deltaR <- sd(riskData[,2])/nrow(riskData)
    }
    dupvalues <- table(riskData[,2])
    dupvalues <- as.numeric(names(dupvalues)[dupvalues>1])
    addNoise <- riskData[,2] %in% dupvalues
    if (isProbability)
    {
        addNoise <- addNoise & (riskData[,2] < 0.999) & (riskData[,2] > 0.001)
    }
    riskData[addNoise,2] <- riskData[addNoise,2] + (runif(sum(addNoise))-0.5)*deltaR*1.0e-6

#    print(sum(addNoise))
  }
  
  isRevesed <- 1.0
  
  if (!isProbability)
  {
    if (mean(riskData[riskData[,1]==1,2]) < mean(riskData[riskData[,1]==0,2]))
    {
      riskData[,2] <- -riskData[,2]
      warning("Reversed Sign")
      isRevesed <- -1.0;
      if (!is.null(atThr))  {atThr <- -atThr;}
    }
  }
  
  riskData <- as.data.frame(riskData)
  if (is.null(rownames(riskData)))
  {
    rownames(riskData) <- as.character(c(1:nrow(riskData)))
  }
  totobs <- nrow(riskData)
  
  ## Getting the threshold values for the specified values
  if (is.null(atThr))
  {
    thr_atP <- quantile(riskData[riskData[,1]==0,2],probs=c(atRate)) - deltaR
    if (atRate[1]<0.5)
    {
      thr_atP <- quantile(riskData[riskData[,1]==1,2],probs=c(atRate)) + deltaR
    }
    if (length(atRate)>1)
    {
      if (atRate[2]>=atRate[1])
      {
        thr_atP <- quantile(riskData[riskData[,1]==1,2],probs=c(1.0-atRate)) + deltaR
      }
      if (atRate[1] < 0.5)
      {
        thr_atP <- quantile(riskData[riskData[,1]==1,2],probs=c(atRate)) + deltaR
      }
      if (thr_atP[1] == thr_atP[2])
      {
        atRate <- atRate[1]
        thr_atP <- thr_atP[1]
      }
    }
  }
  else
  {
    thr_atP <- atThr
    atRate <- atThr
  }
  
  
  
  
        
  numberofEvents <- sum(riskData[,1])
  numberofNoEvents <- sum(riskData[,1]==0)
  
  samplePrevalence <- numberofEvents/totobs;
  ExpectedNoEventsGain <- 1.0
  pre <- samplePrevalence
  
  if (!is.null(ExpectedPrevalence))
  {
    ExpectedNoEvents <- (1.0-ExpectedPrevalence)*totobs*samplePrevalence/ExpectedPrevalence
    ExpectedNoEventsGain <- ExpectedNoEvents/numberofNoEvents
  }
  else
  {
#    print(ExpectedPrevalence)
    if (!is.null(riskTimeInterval) && !is.null(timetoEvent))
    {
      samplePrevalence <- sum(timetoEvent < riskTimeInterval & riskData[,1]==1)/totobs
      pre <- samplePrevalence
    }
    ExpectedPrevalence <- samplePrevalence
 #   print(ExpectedPrevalence)
  }
  
  minRiskAtEvent <- min(riskData[riskData[,1]==1,2])
  thrsWhitinEvents <- riskData[riskData[,2] >= minRiskAtEvent,2]
  names(thrsWhitinEvents) <- rownames(riskData[riskData[,2] >= minRiskAtEvent,])
  thrsWhitinEvents <- thrsWhitinEvents[order(thrsWhitinEvents)]
  nobs <- length(thrsWhitinEvents)
  RR <- numeric(nobs)
  SEN <- numeric(nobs)
  SPE <- numeric(nobs)
  isEvent <- riskData[names(thrsWhitinEvents),1]
  PPV <- numeric(nobs)
  netBenefit <- numeric(nobs)
  URCI <- numeric(nobs)
  LRCI <- numeric(nobs)
  BACC <- numeric(nobs)
  names(RR) <- names(thrsWhitinEvents)
  names(SEN) <- names(thrsWhitinEvents)
  names(SPE) <- names(thrsWhitinEvents)
  names(isEvent) <- names(thrsWhitinEvents)
  names(PPV) <- names(thrsWhitinEvents)
  names(URCI) <- names(thrsWhitinEvents)
  names(LRCI) <- names(thrsWhitinEvents)
  names(BACC) <- names(thrsWhitinEvents)
  names(netBenefit) <- names(thrsWhitinEvents)
  
  mxthr <- max(riskData[,2])
  mithr <- min(riskData[,2])
  
  idx <- 1;
  noMoreLowEventsIdx <- minRiskAtEvent
  for (thr in thrsWhitinEvents)
  {
    
    atLowRisk <- riskData[,2] < thr
    atHighRisk <- !atLowRisk
    LowEvents <- sum(riskData[atLowRisk,1])
    if (LowEvents==0) 
    {
      LowEvents <- 0.1;
      noMoreLowEventsIdx <- thr
    }
    HighEvents <- sum(riskData[atHighRisk,1]);
    SEN[idx] <-  HighEvents/numberofEvents;
    SPE[idx] <- sum(riskData[atLowRisk,1]==0)/numberofNoEvents;
    BACC[idx] <- (SEN[idx] + SPE[idx])/2
    n1 <- sum(atHighRisk)
    n2 <- sum(atLowRisk)
    PPV[idx] <- HighEvents/n1
    LowFraction <- LowEvents/n2;
    HighFraction <- HighEvents/n1;
    x1 <- HighEvents
    x2 <- LowEvents
    if ((n2/numberofNoEvents) > 0.01)
    {
      RR[idx] <- HighFraction/LowFraction
      cidelta <- 1.96*sqrt((n1-x1)/(x1*n1)+(n2-x2)/(x2*n2))
      URCI[idx] <- exp(log(RR[idx])+cidelta)
      LRCI[idx] <- exp(log(RR[idx])-cidelta)
    }
    else
    {
      if (idx>1)
      {
        RR[idx] <- RR[idx-1]
      }
      else
      {
        RR[idx] <- 1.0;
      }
    }
    if (isProbability)
    {
      netBenefit[idx] <- (HighEvents - ExpectedNoEventsGain*sum(riskData[atHighRisk,1]==0)*thrsWhitinEvents[idx]/(1.0001-thrsWhitinEvents[idx]))/totobs;
    }
    idx <- idx + 1;
  }
  colors <- heat.colors(10)
  rgbcolors <- rainbow(10)

  tmop <- par(no.readonly = TRUE)
  

  
  if (isProbability)
  {
    ## Observed vs Cumulative plot
#    par(pty='s',mfrow=c(1,2),cex=0.5)
    par(pty='s')
    xdata <- riskData[order(-riskData[,2]),]
    observed <- numeric(nrow(xdata))
    expected <- numeric(nrow(xdata))
    observed[1] <- xdata[1,1]
    expected[1] <- xdata[1,2]
    pgzero <- xdata[,2]
    
    if (!is.null(timetoEvent) && !is.null(riskTimeInterval))
    {
      atRisktime <- 3.0*mean(timetoEvent[riskData[,1]==1])
      otime <- timetoEvent[order(-riskData[,2])]
      isnoevent <- (otime < atRisktime) & (xdata[,1]==0) # Adjust only if short time to event on censored events
      hazard <- -log(1.00-pgzero[isnoevent])
      pgzero[isnoevent] <- 1.0 - exp(-hazard*otime[isnoevent]/riskTimeInterval)
    }


    for (idx in c(2:nrow(xdata)))
    {
      if (xdata[idx,1]==1)
      {
        expected[idx] <- expected[idx-1] + pgzero[idx];
      }
      else
      {
        expected[idx] <- expected[idx-1] + pgzero[idx]*ExpectedNoEventsGain;
      }
      observed[idx] <- observed[idx-1] + xdata[idx,1];
    }
    totObserved <- sum(xdata[,1])
    tokeep <- (expected > 0.10*totObserved)
    tokeep <- tokeep & (observed > 0.10*totObserved)
    CumulativeOvs <- as.data.frame(cbind(Observed=observed,Cumulative=expected,Included=tokeep))
    
    OEratio <- observed[tokeep]/expected[tokeep]
    OAcum95ci <- c(mean=mean(OEratio),metric95ci(OEratio))
    
    observedCI <- stats::poisson.test(max(observed[tokeep]),max(expected[tokeep]), conf.level = 0.95 )
    OARatio <- observedCI
    OARatio$estimate <- c(OARatio$estimate,OARatio$conf.int,OARatio$p.value)
    names(OARatio$estimate) <- c("O/A","Low","Upper","p.value")
    
    
    rownames(CumulativeOvs) <- rownames(xdata)
    maxobs <- max(c(observed,expected))
    
    if (plotRR)
    {
      plot(expected,observed,
           ylab="Observed",
           xlab="Cumulative Probability",
           xlim=c(0,maxobs),
           ylim=c(0,maxobs),
           main=paste("Cumulative vs. Observed:",title),
           pch=c(5+14*xdata[,1]),
           col=c(1+xdata[,1]),
           cex=c(0.2+xdata[,1]))
           se <- 2*sqrt(observed)
      errbar(expected,observed,observed-se,observed+se,add=TRUE,pch=0,errbar.col="gray",cex=0.25)
      points(expected,observed,
              pch=c(5+14*xdata[,1]),
                 col=c(1+xdata[,1]),
                 cex=c(0.2+xdata[,1]))

      lines(x=c(0,maxobs),y=c(0,maxobs),lty=2)
      legend("bottomright",legend=c("Event","Expected"),
             lty=c(-1,2),
             pch=c(19,-1),
             col=c(2,1),cex=0.75)
      par(tmop)
    }

    
    
    ## Decision curve analysis
    pshape <- 4 + 12*isEvent
    xmax <- min(quantile(thrsWhitinEvents,probs=c(0.95),0.95))
    ymin <- min(quantile(netBenefit,probs=c(0.10)),0)
    DCA <- as.data.frame(cbind(Thrs=thrsWhitinEvents,NetBenefit=netBenefit))
    rownames(DCA) <- names(thrsWhitinEvents)
    

    if (plotRR)
    {
      plot(thrsWhitinEvents,netBenefit,main=paste("Decision Curve Analysis:",title),ylab="Net Benefit",xlab="Threshold",
         ylim=c(ymin,pre),
         xlim=c(0,xmax),
         pch=pshape,col=rgbcolors[1+floor(10*(1.0-SEN))],cex=(0.35 + isEvent))
      fiveper <- as.integer(0.05*length(netBenefit)+0.5)
      range <- c(fiveper:length(netBenefit)-fiveper)
      lfit <-try(loess(netBenefit[range]~thrsWhitinEvents[range],span=0.5));
      if (!inherits(lfit,"try-error"))
      {
        plx <- try(predict(lfit,se=TRUE))
        if (!inherits(plx,"try-error"))
        {
          lines(thrsWhitinEvents[range],plx$fit,lty=1)
          lines(thrsWhitinEvents[range],plx$fit - qt(0.975,plx$df)*plx$se, lty=2)
          lines(thrsWhitinEvents[range],plx$fit + qt(0.975,plx$df)*plx$se, lty=2)
        }
      }
      abline(h=0,col="blue")
      lines(x=c(0,ExpectedPrevalence),y=c(pre,0),col="red")
      legtxt <- sprintf("%3.1f",c(5:0)/5)
      corrplot::colorlegend(rgbcolors, legtxt,xlim=c(0.92*xmax,0.98*xmax),ylim=c(pre*0.45,pre*0.1),cex=0.75)
      text(0.92*xmax,pre*0.5,"SEN")
      
      thrlty <- c(1)
      thrLey <- "H. Thr"
      abline(v=thr_atP[1],col="gray",lty=thrlty)
      if (length(thr_atP)>1) 
      {
        abline(v=thr_atP[2],col="gray",lty=2)
        thrlty <- c(thrlty,2)
        thrLey <- c(thrLey,"L. Thr")
      }
      
      legend("topright",legend=c("No Event","Event","loess fit",thrLey),
             pch=c(4,16,-1,-1,-1),
             col=c(1,1,1,"gray","gray"),
             lty=c(-1,-1,1,thrlty)
             )
      legend("bottomleft",legend=c("Treat All","Treat None"),
             lty=c(1,1),
             col=c("red","blue"),cex=0.5)
      par(tmop)
    }
    
  }
  ## Relative Risk plot
  
  ymax <- quantile(RR,probs=c(0.99))
  pshape <- 2 + 14*isEvent
  RRData <- as.data.frame(cbind(thr=thrsWhitinEvents,
                                Sensitivity=SEN,
                                Specificity=SPE,
                                BACC  = BACC,
                                PPV=PPV,
                                RR=RR,
                                LRR=LRCI,
                                URR=URCI,
                                isEvent=isEvent))
  rownames(RRData) <- names(thrsWhitinEvents)
  lfit <- try(loess(RR~SEN,span=0.5));
#  print(thr_atP)
  ## Get the location of the thresholds
  maxBACC <- max(BACC)
  RRval <- RR[thrsWhitinEvents <= thr_atP[1] & SPE > 0.10]
  thrsval <- thrsWhitinEvents[thrsWhitinEvents <= thr_atP[1] & SPE > 0.10]
  maxRR <- max(RRval)
  thr_values <- c(thr_atP,
                  at_max_BACC=min(thrsWhitinEvents[BACC==maxBACC]),
                  at_max_RR=min(thrsval[RRval==maxRR]),
                  atSPE100=noMoreLowEventsIdx)
  if (isProbability)
  {
    thr_values <- c(thr_values,at_0.5=0.5)
  }
  
  thrLoc <- rep(1,length(thr_values))
#  print(thr_values)
  for (tthr in c(1:length(thr_values)))
  {
    thrLoc[tthr] <- which.min(abs(thrsWhitinEvents-thr_values[tthr]))
  }
#  print(thrLoc)
  thrPoints <- as.data.frame(cbind(isRevesed*thrsWhitinEvents[thrLoc],RR[thrLoc],LRCI[thrLoc],URCI[thrLoc],SEN[thrLoc],SPE[thrLoc],BACC[thrLoc]))
  colnames(thrPoints) <- c("Thr","RR","RR_LCI","RR_UCI","SEN","SPE","BACC")
  if (isProbability)
  {
    thrPoints$NetBenefit <- netBenefit[thrLoc]
    rownames(thrPoints) <- c(paste("@",round(atRate,3),sep=":"),"@MAX_BACC","@MAX_RR","@SPE100","p(0.5)")
  }
  else
  {
    rownames(thrPoints) <- c(paste("@",round(atRate,3),sep=":"),"@MAX_BACC","@MAX_RR","@SPE100")
  }
#  print(thrPoints)
  
  ## Sensitivity, Specificity and RR at top threshold
  lowRisk <- riskData[,2] < thr_atP[1]
  LowEventsFrac <- sum(riskData[lowRisk,1])/sum(lowRisk)
  HighEventsFrac <- sum(riskData[!lowRisk,1])/sum(!lowRisk)
  sensitivity=sum(riskData[!lowRisk,1])/numberofEvents
  who <- names(SEN)[SEN==sensitivity]
  RRAtSen <- c(est=mean(RR[who]),lower=mean(LRCI[who]),upper=mean(URCI[who]))
  
  specificity=sum(riskData[lowRisk,1]==0)/numberofNoEvents
  
  if (plotRR)
  {
    ypmax <- max(c(quantile(URCI,probs=c(0.95),na.rm = TRUE),ymax))

    par(mfrow=c(1,1))
  
    plot(SEN,RR,cex=(0.35 + PPV),
       pch=pshape,
       col=colors[1+floor(10*(1.0-SPE))],
       xlim=c(0,1.15),
#       ylim=c(1.0,ymax+0.5),
       ylim=c(1.0,ypmax*1.2),
#       log="y",
       main=paste("Relative Risk:",title))
    atevent <- isEvent==1
    errbar(SEN[atevent],RR[atevent],URCI[atevent],LRCI[atevent],add=TRUE,pch=0,errbar.col="gray",cex=0.25)
    points(SEN,RR,cex=(0.35 + PPV),
         pch=pshape,
         col=colors[1+floor(10*(1.0-SPE))],
    )  
#    lfit <- try(loess(RR~SEN,span=0.5));
    if (!inherits(lfit,"try-error"))
    {
      plx <- try(predict(lfit,se=TRUE))
      if (!inherits(plx,"try-error"))
      {
        lines(SEN,plx$fit,lty=1)
        lines(SEN,plx$fit - qt(0.975,plx$df)*plx$se, lty=2)
        lines(SEN,plx$fit + qt(0.975,plx$df)*plx$se, lty=2)
      }
    }
  
    legtxt <- sprintf("%3.1f",c(5:0)/5)
    corrplot::colorlegend(heat.colors(10), legtxt,xlim=c(1.05,1.175),ylim=c((ypmax-0.75)*0.45+1.0,ypmax/10+1.0),cex=0.75)
    text(1.075,1,"SPE")
    sizp <- c(5:1)/5.0 + 0.35
    legend("topright",legend=legtxt[1:5],pch=16,cex=sizp)
    legend("bottomleft",legend=c("No","Yes","Loess Fit"),pch=c(2,16,-1),lty=c(0,0,1),cex=0.75)
    text(0.95,ypmax+0.5,"PPV->")


    abline(v=sensitivity,col="blue")
    text(x=sensitivity,y=ypmax*1.2,sprintf("Index(%3.2f)=%4.3f",specificity,isRevesed*thr_atP[1]),pos=4 - 2*(sensitivity>0.5) ,cex=0.7)
    text(x=sensitivity,y=0.9*(ypmax*1.2-1.0)+1.0,
         sprintf("RR(%3.2f)=%4.3f",
                 sensitivity,
                 RRAtSen[1]),
         pos=4 - 2*(sensitivity>0.5),cex=0.7)
    par(tmop)
  }
  
  ## ROC At threshold
  
  if (plotRR)
  {
    par(tmop)
    ROCAnalysis <- predictionStats_binary(riskData,
                                        plotname=paste("ROC:",title),
                                        thr=thr_atP[1],cex=0.85)

    par(tmop)
  }
  else
  {
    ROCAnalysis <- predictionStats_binary(riskData,
                                          thr=thr_atP[1])
  }
  surfit <- NULL
  surdif <- NULL
  LogRankE <- NULL
  cstat <- NULL
  timetoEventData <- NULL
  if (!is.null(timetoEvent))
  {
      timetoEventData <- as.data.frame(cbind(eStatus=riskData[,1],
                           class=1*(riskData[,2]>=thr_atP[1]),
                           eTime=timetoEvent,
                           risk=riskData[,2])
                           )
      if (length(thr_atP)>1)
      {
        timetoEventData$class <- 1*(riskData[,2]>=thr_atP[1]) + 1*(riskData[,2]>=thr_atP[2])
      }
      paletteplot <- c("green", "red")
      if (isRevesed > 0)
      {
        labelsplot <- c("Low Risk",sprintf("High Risk > %5.3f",thr_atP[1]));
        if (length(thr_atP)>1)
        {
          labelsplot <- c(sprintf("Low Risk < %5.3f",thr_atP[2]),sprintf("%5.3f <= Risk < %5.3f",thr_atP[2],thr_atP[1]),
                          sprintf("High Risk >= %5.3f",thr_atP[1]));
          paletteplot <- c("green", "pink","red")
        }
      }
      else
      {
        labelsplot <- c("Low Risk",sprintf("High Risk < %5.3f",-thr_atP[1]));
        if (length(thr_atP)>1)
        {
          labelsplot <- c(sprintf("Low Risk > %5.3f",-thr_atP[2]),sprintf("%5.3f >= Risk > %5.3f",-thr_atP[2],-thr_atP[1]),
                          sprintf("High Risk <= %5.3f",-thr_atP[1]));
          paletteplot <- c("green", "pink","red")
        }
        
      }
      if (isProbability)
      {
        ## Time Plot
        prisk <- riskData[,2]
#        prisk[prisk > 0.999999] <- 0.999999
#        timetoEventData$lammda <- -log(1.0-prisk) # From probability to expected events per time interval
        timetoEventData$lammda <- expectedEventsPerInterval(prisk)
        aliveEvents <- timetoEventData
        aliveEvents <- aliveEvents[order(-aliveEvents$risk),]
        aliveEvents <- aliveEvents[order(aliveEvents$eTime),]
        
        if (!is.null(riskTimeInterval))
        {
          timeInterval <- riskTimeInterval
        }
        else
        {
          timeInterval <- 2*mean(subset(aliveEvents,
                                aliveEvents$eStatus==1)$eTime);
        }
        timetoEventData$expectedTime <- meanTimeToEvent(prisk,timeInterval)
        attr(timetoEventData,"ClassNames") <- labelsplot
        timed <- unique(aliveEvents[aliveEvents$eStatus==1,"eTime"])
        maxtime <- max(timed)
        Observed <- c(1:length(timed))
        Expected <- Observed
        cClass <- Observed
        aliveEvents[aliveEvents$eStatus == 0,"lammda"] <- aliveEvents[aliveEvents$eStatus == 0,"lammda"]*ExpectedNoEventsGain
        passAcum <- 0;
        lastObs <- 0;
        allTimes <- aliveEvents$eTime
        lasttime <- 0

        for (idx in c(1:length(timed)))
        {
          whosum <- (allTimes <= timed[idx]) & (allTimes > lasttime)
          deltatime <- (timed[idx] - lasttime)
          Observed[idx] <- lastObs + sum(aliveEvents[whosum,"eStatus"])
          cClass[idx] <- round(mean(aliveEvents[whosum,"class"]),0)
          lastObs <- Observed[idx]
          eevents <- sum(aliveEvents[allTimes > lasttime,"lammda"])*deltatime/timeInterval
          Expected[idx] <- passAcum + eevents;
          passAcum <- Expected[idx]
          lasttime <- timed[idx]
        }
        
        totObserved <- max(Observed)
        maxevents <- max(c(Observed,Expected))
        tokeep <- (Expected > 0.10*totObserved)
        tokeep <- tokeep & (Observed > 0.10*totObserved)
        OEData <- as.data.frame(cbind(time=timed,Observed=Observed,Expected=Expected,Included=tokeep,class=cClass))
        OEratio <- Observed[tokeep]/Expected[tokeep]
        OE95ci <- c(mean=mean(OEratio),metric95ci(OEratio))

        observedCI <- stats::poisson.test(max(Observed[tokeep]),max(Expected[tokeep]), conf.level = 0.95 )
        OERatio <- observedCI
        OERatio$estimate <- c(OERatio$estimate,OERatio$conf.int,OERatio$p.value)
        ## Estimation of O/E for the threshold values
        names(OERatio$estimate) <- c("O/E","Low","Upper","p.value")
        procat <- timetoEventData$class
        values <- unique(procat)
        values <- values[order(values)]
        names(values) <- c("low",names(thr_atP))
        obs <- numeric()
        expected <- numeric()
        lci <- numeric()
        uci <- numeric()
        lciOE <- numeric()
        uciOE <- numeric()
        pval <- numeric()
        hazards <- timetoEventData$lammda*timetoEventData$eTime/timeInterval
        for (ct in values)
        {
          totobs <- sum(riskData[procat==ct,1])
          obs <-c(obs,totobs)
          expe <- sum(hazards[procat==ct]);
          expected <- c(expected,expe)
          pt <- stats::poisson.test(totobs,1)
          lci <- c(lci,pt$conf.int[1])
          uci <- c(uci,pt$conf.int[2])
          pt2 <- stats::poisson.test(totobs,expe)
          lciOE <- c(lciOE,pt2$conf.int[1])
          uciOE <- c(uciOE,pt2$conf.int[2])
          pval <- c(pval,pt2$p.value)
        }
        totobservedCI <- stats::poisson.test(max(Observed),1, conf.level = 0.95 )
        totalEstimated <- c(max(Observed),totobservedCI$conf.int,max(Expected),max(Observed)/max(Expected),OERatio$conf.int,OERatio$p.value)
        OERatio$atThrEstimates <- as.data.frame(cbind(obs,lci,uci,expected,obs/expected,lciOE,uciOE,pval))
        OERatio$atThrEstimates <- rbind(totalEstimated,OERatio$atThrEstimates)
        colnames(OERatio$atThrEstimates) <- c("Observed","L.CI","H.CI","Expected","O/E","Low","Upper","pvalue")
        rownames(OERatio$atThrEstimates) <- c("Total",names(values))
        ## Now lets plot
        if (plotRR)
        {
          par(mfrow=c(1,1))
          plot(timed,Expected,pch=4,type="b",cex=0.5,
             main=paste("Time vs. Events:",title),
             ylab="Events",
             xlab="Time",
             ylim=c(0,1.05*maxevents),
             xlim=c(0,1.05*maxtime),
             )
          se <- 2*sqrt(Observed)
          errbar(timed,Observed,Observed-se,Observed+se,add=TRUE,pch=0,errbar.col="gray",cex=0.25)
          points(timed,Expected,pch=4,type="b",cex=0.5)
          points(timed,Observed,pch=1,col=paletteplot[1+cClass])
          legend("topleft",legend=c("Expected",labelsplot),pch=c(4,1,1,1),lty=c(1,0,0,0),col=c(1,paletteplot),cex=0.80)
        }
      }
        ## Survival plot
      LogRankE <- EmpiricalSurvDiff(times=timetoEventData$eTime,
                  status=timetoEventData$eStatus,
                  groups=1*(timetoEventData$class > 0),
                  plots=FALSE,main=paste("Kaplan-Meier:",title))
      
      surfit <- survival::survfit(survival::Surv(eTime,eStatus)~class,data = timetoEventData)
      surdif <- survival::survdiff(survival::Surv(eTime,eStatus)~class,data = timetoEventData)
      cstat <- rcorr.cens(-timetoEventData$risk,survival::Surv(timetoEventData$eTime,timetoEventData$eStatus))
      cstatCI <- c(mean=cstat[1],concordance95ci(as.data.frame(cbind(times=timetoEventData$eTime,
                                       status=timetoEventData$eStatus,
                                       preds=-timetoEventData$risk))))
      cstat <- as.list(cstat)
      cstat$cstatCI <- cstatCI
      
      if (plotRR)
      {
        par(mfrow=c(1,1))
        graph <- survminer::ggsurvplot(surfit,
                                       data = timetoEventData, 
                                       conf.int = TRUE, 
                                       legend.labs = labelsplot,
                                       palette = paletteplot,
                                       ylim = ysurvlim,
                                       ggtheme = ggplot2::theme_bw() + 
                                         ggplot2::theme(plot.title = ggplot2::element_text(
                                           hjust = 0.5, face = "bold", size = 16)),
                        title = paste("Kaplan-Meier:",title),
                        risk.table = TRUE,
                        tables.height = 0.2,
                        tables.theme = survminer::theme_cleantable())
        print(graph)
        par(tmop)
        
      }
    

  }
  
  result <- list(CumulativeOvs=CumulativeOvs,
                 OEData=OEData,
                 DCA=DCA,
                 RRData=RRData,
                 timetoEventData=timetoEventData,
                 keyPoints=thrPoints,
                 OERatio=OERatio,
                 OE95ci=OE95ci,
                 OARatio=OARatio,
                 OAcum95ci=OAcum95ci,
                 fit=lfit,
                 ROCAnalysis=ROCAnalysis,
                 prevalence=pre,
                 thr_atP= isRevesed*thr_atP,
                 c.index=cstat,
                 surfit=surfit,
                 surdif=surdif,
                 LogRankE=LogRankE
                 )
  return (result)
}