R/oxy-plotHazSurv.R

In nphPower: Sample Size Calculation under Non-Proportional Hazards

#' @title Graphic Display of Hazard and Survival Function
#' @description Plot the hazard and survival function of the
#' of control group (from weibull or loglogistic distribution) and treatment group
#' (derived from an arbitrary hazard ratio function)
#'
#'
#' @param bsl_dist a text must be one of (\code{"weibull"},
#' \code{"loglogistic"}) distribution, specified for the control group
#' @param param a vector of length 2, specifying the shape and rate (1/scale)
#' parameter of the \code{bsl_dist} distribution, Default: c(1.2, 0.03)
#' @param fun_list a list of hazard ratio functions comparing treatment group
#' and control group
#' @param end a value specifying the duration of the curve
#' @param tit a vector specifying the titles of each graph, Default: c("Hazard Function", "Survival Function")
#' @param pos a graphic parameter in the form of c(nr,nc). Subsequent
#' figures will be drawn in an nr-by-nc array, Default: c(1, 2)
#' @param hlegend.loc a text indicating the position of legend
#' for the hazard plot. Default: "bottomleft"
#' @return
#' graphics of hazard and survival functions
#' @examples
#' # proportional hazards
#'plotHazSurv(
#'  bsl_dist=c("weibull")
#'  ,param=c(1.2,1/30)
#'  ,fun_list=list(function(x){x^0*0.7})
#'  ,40
#'  ,tit= c("Hazard Function","Survival Function")
#'  ,pos=c(1,2)
#')
#' # crossing hazards
#'plotHazSurv(
#'  bsl_dist=c("weibull")
#'  ,param=c(1.2,1/30)
#'  ,fun_list=list(function(x){1.3*(x<10)+(x>=10)*0.7})
#'  ,40
#'  ,tit= c("Hazard Function","Survival Function")
#'  ,pos=c(1,2)
#')
#' @rdname plotHazSurv
#' @export
#' @importFrom graphics par lines mtext legend
#' @importFrom stats integrate
plotHazSurv <- function(
  bsl_dist=c("weibull","loglogistic")
  ,param=c(1.2,0.03)
  ,fun_list
  ,end
  ,tit= c("Hazard Function","Survival Function")
  ,pos=c(1,2)
){
  t <- seq(0.5,end,length=300)
  fn <- length(tit)

  a <- param[1]
  b <- param[2]
  old.par <- par(no.readonly = TRUE) # all par settings which
  # could be changed.
  on.exit(par(old.par))

 # old <- options(mfrow = pos, mar=c(4,3,3,1))
 # on.exit(options(old), add = TRUE)

  graphics::par(mar=c(4,3,3,1),mfrow=pos)

  for (i in 1:length(fun_list)){
    ##-----------------draw the hazard function------------------##
    if (bsl_dist=="weibull"){
      hr1 <- a*b*(b*t)^(a-1)

    }else if (bsl_dist=="loglogistic"){
      hr1 <- a/b*(t/b)^(a-1)/(1+(t/b)^a)
    }
    Hft <- function(t){
      if (t==0) {1
      }else {Hf(t)}
    }
    hr <- fun_list[[i]](t)*hr1
    plot(t,hr,cex=0.3,ylab="h(t)",cex.main=1,
         ylim = c(0,max(hr,hr1)*1.1),type="l",lty=2,col="blue",xlab="Time")
    graphics::lines(t,hr1, lty = 1,col="red")
    dis <-tit[1]
    graphics::mtext(dis, side=3,line=1,adj=0,cex=0.8)
  graphics::par(xpd=TRUE)
    graphics::legend("bottomleft",legend = c("Treatment Group", "Control Group"), lty=c(1,2), lwd = 1,
           col=c("red","blue"), xpd = TRUE, horiz = FALSE, cex = 0.5)


  }
  for (i in 1:length(fun_list)){
    ##-----------------draw the survival function------------------##
    if (bsl_dist=="weibull"){
      S1 <- exp(-(b*t)^a)
      Hf <- function(t){exp(-1*a*b*stats::integrate( function(x){(x*b)^(a-1)*fun_list[[i]](x)},0,t)$value)}
    }else if (bsl_dist=="loglogistic"){
      S1 <- b^a/(b^a+t^a)
      Hf <- function(t){exp(-1* a/b*stats::integrate( function(x){(x/b)^{a-1}/(1+(x/b)^a)*fun_list[[i]](x)},0,t)$value)}
    }
    Hft <- function(t){
      if (t==0) {1
      }else {Hf(t)}
    }
    S0 <- as.vector(unlist(lapply(t,Hft)))
    plot(t,S0,cex=0.3,cex.main=1,ylab = "S(t)",type="l",
         ylim = c(0,1),lty=2,col="blue",xlab="Time")
    graphics::lines(t,S1,lty=1,col="red")
    dis <-tit[2]
    graphics::mtext(dis, side=3,line=1,adj=0,cex=0.8)
    graphics::legend("topright",legend = c("Treatment Group", "Control Group"), lty=c(1,2), lwd = 1,
           col=c("red","blue"), xpd = TRUE, horiz = FALSE, cex = 0.5)


  }
 # graphics::par(mfrow=c(1,1))
  invisible()

}