R/plot.sim.R

In keaven/nphsim: Non proportional hazards sample size and simulation

#' Plot Simulated Trial Outcome Summaries
#'
#' Generate boxplot, display summary of number of death, timing, boundary crossing probability, and user requested simulation-level information 
#'
#' \code{plotsim()} takes the \code{result} data table from the return value of the \code{simtest()} function 
#' (or a table table that has the same structure as \code{result}) and generate the boxplot and display summary statistics.
#'
#' @param d data.table "result" generated by simtest function (or a table table that has the same structure as \code{result})
#' @param y vector of the name(s) of the y varaible(s) for the boxplot. e.g. c("hr")
#' @param dg number of digits for median of y below the table. NULL will disable the display of the y statistics.
#' @param yt y-axis label 
#' @param logs \code{TRUE}=log scale for y-axis
#' @param b whether to show boundary crossing probability at each analysis. This require xeff and xfut exist in the \code{result} data table
#' @param v height of graph area and the summary area below.
#' @examples
#' # Use a gsSurv object as both input for simulation and testing. 
#' # A logrank p-value and HR from cox model is reported.
#' library(gsDesign)
#' gs <- gsSurv (k = 3, test.type = 4, alpha = 0.025, beta = 0.05, timing = c( 0.5,0.75 ), 
#'               sfu = sfHSD , sfupar = c( -4 ), sfl = sfHSD, sflpar = c( -12 ), 
#'               lambdaC = log(2) / 6, hr = 0.65, hr0 = 1, eta = 0.01, 
#'               gamma = c( 2.5,5,7.5,10 ), R = c( 2,2,2,6 ) , S = NULL , T = 15 , minfup = 3 , ratio = 1) 
#' sim1 <- nphsim(nsim=1000,d=gs)
#' test1 <- simtest(x=sim1,anatype='event',method='LR', d=gs)
#' test1$result
#' ## Not useful, demonstration purpose only
#' plotsim(test1$result,y=c("hr","pval"),dg=2,yt="Hazard Ratio",b=1,v=c(2,0.8))  
#' 
#' @export
#' @import ggplot2 grid data.table
plotsim<-function(d
                  ,y
                  ,dg=NULL
                  ,yt="Y Axis Title"
                  ,logs=NULL
                  ,b=NULL
                  ,v=c(2,0.5) 
                  ) 
{
  a=copy(d)
  a[,x:=factor(analysis)]
  ey=substitute(y)

  ### generate data for the boxplot from simulation result
  fig<-melt(a,id.vars=c("sim","analysis"), measure.vars = eval(ey))
  ### generate data for the median (5%,95%) of death and calendar timing from simulation result
  dat<-a[,.(Timing=paste0(round(quantile(t,0.5),1)," (",round(quantile(t,0.05),1),", ",round(quantile(t,0.95),1),")"),
            "No. of Death"=paste0(round(quantile(D,0.5))," (",round(quantile(D,0.05)),", ",round(quantile(D,0.95)),")")),by=analysis]
  dat.m<-melt(dat,id.vars=c("analysis"), measure.vars = c("Timing","No. of Death"))

  ### generate boundary crossing probability controlled by parameter "b"
  if (!is.null(b)) {
    dat1 <-
      a[, .("Upper Bound" = paste0(round(sum(xeff) * 100 / .N, 0), "%"),
            "Lower Bound" = paste0(round(sum(xfut) * 100 / .N, 0), "%")), by = analysis]
    dat1.m <-
      melt(
        dat1,
        id.vars = c("analysis"),
        measure.vars = c("Upper Bound", "Lower Bound")
      )

    dat.m <- rbind(dat.m, dat1.m)
  }

  ### generate median (5%, 95%) of y value based on the simulation result, controlled by "dg"
  if (!is.null(dg)) {
    #dat2 <-a[, lapply(lapply(.SD,median),round,digits=dg), by = analysis, .SDcol=eval(ey)]
    dat2 <-a[, paste0("_",eval(ey),"_"):=lapply(.SD,function(x){paste0(round(median(x,na.rm=TRUE),digits=dg), " (",
                                                                       round(quantile(x,0.05,na.rm=TRUE),digits=dg), ", ",
                                                                       round(quantile(x,0.95,na.rm=TRUE),digits=dg), ")"
    )}),
             by = analysis, .SDcol=eval(ey)]
    dat2.m <-
      melt(
        dat2,
        id.vars = c("analysis"),
        measure.vars = paste0("_",eval(ey),"_")
      )

    dat.m <- rbind(dat.m, dat2.m)
  }
  none <- element_blank()
  g <- ggplot(fig, aes(factor(analysis), value))
  g1 <- g + geom_boxplot(aes(fill = factor(variable)))  +
    labs(y = yt, x = "") +
    theme(panel.border = element_rect(fill = NA), legend.position="none")

  if (isTRUE(logs)) {
    g1<-g1+scale_y_log10()
  }

  if (length(eval(ey))>1) {
    g1<-g1+theme(legend.position = "top",legend.title=none)
  }

  data_table <-
    ggplot(dat.m, aes(
      x = factor(analysis),
      y = factor(
        variable,
        levels = c(rev(paste0("_",eval(ey),"_")),"Lower Bound", "Upper Bound", "Timing", "No. of Death")
      ),
      label = value
    )) +
    geom_text(size = 3) +
    scale_y_discrete() + theme_bw()  +
    theme(
      panel.grid.major = none,
      legend.position = "none",
      panel.border = none,
      axis.text.x = none,
      axis.ticks = none
    ) + theme(plot.margin = unit(c(-0.5,
                                   1, 0, 0), "lines")) + labs(y = "", x = "Analysis Sequence         Median (5%,95%) from Simulation")



  Layout <- grid.layout(nrow = 2, ncol = 1, heights = unit(v, c("null", "null")))
  #grid.show.layout(Layout)
  vplayout <- function(...) {
    grid.newpage()
    pushViewport(viewport(layout = Layout))
  }

  subplot <- function(x, y) viewport(layout.pos.row = x,
                                     layout.pos.col = y)

  mmplot <- function(a, b) {
    vplayout()
    print(a, vp = subplot(1, 1))
    print(b, vp = subplot(2, 1))
  }


  p<-mmplot(g1, data_table)
  return(invisible())
}