R/functions.R

In LittleBeannie/demo: a R package for demo in Roche inernship

#' Title to predict the final OS based on the interim OS
#'
#' @param formula something like Surv(time, status) ~ trt
#' @param interim.OS
#' @param dataset
#' @param family can be c("exponential","weibull","lognormal")
#'
#' @return
#' @export
#'
#' @import survival
#' @import flexsurv
#' @examples
predict_finalOS <- function(formula, interim.OS = 10, dataset = veteran,family = "exponential"){
  if(family == "exp"){
    fit = flexsurvreg(formula, data = dataset, dist = family)
    rate = exp(fitg$res.t[1])
    final.OS.mean = interim.OS + 1/rate
    final.OS.sd = 1/rate^2
    return( data.frame(mean = final.OS.mean, sd = final.OS.sd) )
  }else if(family == "weibull"){
    fit = flexsurvreg(formula, data = dataset, dist = family)
    shape = exp(fit$coefficients[1])
    scale = exp(fit$coefficients[2])
    k = shape
    lambda = 1/scale^shape
    final.OS.mean =  exp(lambda*interim.OS^k)/lambda^(1/k)*gamma(1/k+1)*(1 - pgamma(lambda*interim.OS^k,shape = 1/k+1,scale = 1))
    final.OS.sd = exp(lambda*interim.OS^k)/lambda^(2/k) *(1- pgamma(lambda*interim.OS^k,shape = 2/k+1,scale = 1))-final.OS.mean^2
    return( data.frame(mean = final.OS.mean, sd = final.OS.sd) )
  }else if(family == "lognormal"){
    fit = flexsurvreg(formula, data = dataset, dist = family)
    mu = exp(fit$coefficients[1])
    sigma = exp(fit$coefficients[2])
    final.OS.mean = exp(2*mu*sigma^2-sigma^2) / (1-pnorm( (log(interim.OS)-mu)/sigma ) ) * (1-pnorm( (log(interim.OS)-sigma^2+mu)/sigma ))
    final.OS.sd = exp(4*mu*sigma^2-4*sigma^2) / (1-pnorm( (log(interim.OS)-mu)/sigma ) ) * (1-pnorm( (log(interim.OS)-sigma^2+mu)/sigma )) -  final.OS.mean^2
    return( data.frame(mean = final.OS.mean, sd = final.OS.sd) )
  }else if(family == "log_logistic"){
  }
}



#'---------------------- May 21- May28 ------------------------

#' genalize linear regression
#' @param x
#' @param y
#' @import glmnet
#' @return the beta coeff
#' @export
#'
genlm <- function(x=matrix(rnorm(20),nrow = 4), y=rnorm(4)){
  modelfit <- glmnet::glmnet(x,y,family = "gaussian")
  return(modelfit$beta)
}


#' A brief introduction to this function.
#' @param x1 explain what is x1
#' @param x2
#' @return explain what is returned by this function.
#' @export
#' @examples foo()
foo <- function(x1 = 2, x2 = 21) {
  y <- x1 * x2
  return(y)
}


#' plot the KM curve
#' @param model a Surv(time, status) object
#' @param covariate it can be 1 if all covariates are included, or it can be a single covariate
#' @param dataset it is the dataset we use
#' @param plot_fancy if it is TRUE, then it uses ggplot2, otherwise it use the base graphics
#'
#' @import survival
#' @import ranger
#' @import ggplot2
#' @import ggfortify
#'
#' @return a figure
#' @export
#'
#' @examples
km_plot <- function(model,covariate,dataset = veteran,plot_fancy=T){
  km_fit <- survfit(model ~ covariate, data=dataset)
  if(plot_fancy ==  T){
    autoplot(km_fit)
  } else {
    plot(km_fit, xlab="Days", main = 'Kaplan Meyer Plot')
  }
}




#' compute log-rank test comparing two or more survival curves.
#'
#' @param model a Surv(time, status) object
#' @param covariate it can be 1 if all covariates are included, or it can be a single covariate
#' @param dataset it is the dataset we use
#'
#' @return a summary of the log rank test
#' @export
#' @import survival
#'
#' @examples
log_rank_test <- function(model,covariate, dataset = veteran){
  test_result <- survdiff(model ~ covariate, data = dataset)
  test_result
}


#' cox regression model
#'
#' @param model a Surv(time, status) object
#' @param covariate it can be 1 if all covariates are included, or it can be a single covariate
#' @param dataset it is the dataset we use
#'
#' @return a summary of the cox model
#' @export
#'
#' @import survival
#' @examples
cox_reg <- function(model, covariate,dataset = veteran){
  cox_model <- coxph(model ~ covariate, data = dataset)
  summary(cox_model)
}


#' fit and predict the survival function
#'
#' @param model a Surv(time, status) object
#' @param covariate it can be 1 if all covariates are included, or it can be a single covariate
#' @param dataset it is the dataset we use
#' @param family the paramatric distribution, assumed distribution for y variable.
#' @param futuretime the time we want to predict
#'
#' @return the predicted value of the survival function
#' @export
#'
#' @import survival
#' @examples
curve_predict <- function(model, covariate, dataset = vetran,family="exp",futuretime=1){
  fit <- survreg(model~covariate,data = dataset, dist = family)
  if(family == "exp"){
    lambda <- exp(-fit$coefficients[1])
    St <- function(t){exp(-lambda*t)}
    predicted_value <- St(futuretime)
    return(predicted_value)

  }else if(family == "Weibull"){
    gamma <- fit$scale
    alpha <- exp(-gamma*fit$coefficients[1])
    St <- function(t){ exp(-alpha*t^gamma) }
    predicted_value <- St(futuretime)
    return(predicted_value)
  }

}




#' Title a fancy KM curve plot with p-value in the plot
#'
#' @param model.formula the survival formula
#' @param pvalue.TF  whether showing the p-value or not
#' @param conf.int.TF whether showing the confidence interval or not
#' @param x.lab the label of x-axis
#' @param x.interval the break of x-axis
#' @param risk.table.style c('absolute', 'percentage','abs_pct') show both absolute number or/and  percentage of individuals at risk
#' @param n.censor.plot.TF to plot the number of censored subjects at time t. As suggested by Marcin Kosinski, This is a good additional feedback to survival curves, so that one could realize: how do survival curves look like, what is the number of risk set AND what is the cause that the risk set become smaller: is it caused by events or by censored events?
#' @param surv.median.line.style horizontal/vertical line at median survival using the argument surv.median.line. Allowed values include one of c(“none”, “hv”, “h”, “v”). v: vertical, h:horizontal.
#' @param legend to change the legend labels.
#'
#' @return a km plot
#' @export
#'
#' @import survival
#' @import survminer
#' @examples
km_plot_fancy <- function(model.formula = survfit(Surv(time, status)~ trt, data=veteran),
                          pvalue.TF = TRUE,
                          conf.int.TF = TRUE,
                          x.lab = "Time in Days",
                          x.interval = 200,
                          risk.table.style = "abs",
                          n.censor.plot.TF= TRUE,
                          surv.median.line.style = "hv",
                          legend = c("group 1", "group 2")){
  ggsurvplot(model.formula,
             pval = pvalue.TF,
             conf.int = conf.int.TF,
             xlab = x.lab,
             break.time.by = x.interval,
             risk.table = risk.table.style,
             risk.table.y.test.col =TRUE,
             risk.table.y.test = FALSE,
             ncensor.plot = n.censor.plot.TF,
             surv.median.line = surv.median.line.style,
             legend.labs = legend,
             palette = c("#E7B800", "#2E9FDF")
             )
}








#' Title
#'
#' @param model1.r1 the Surv(time, status) object in the first round in test 1
#' @param cov1.r1
#' @param dataset1.r1
#'
#' @param model2.r1 the Surv(time, status) object in the first round in test 2
#' @param cov2.r1
#' @param dataset2.r1
#'
#' @param model3.r1 the Surv(time, status) object in the first round in test 3
#' @param cov3.r1
#' @param dataset3.r1
#'
#' @param model1.r2 the Surv(time, status) object in the second round in test 1
#' @param cov1.r2
#' @param dataset1.r2
#'
#' @param model2.r2 the Surv(time, status) object in the second round in test 2
#' @param cov2.r2
#' @param dataset2.r2
#'
#' @param model3.r2 the Surv(time, status) object in the second round in test 3
#' @param cov3.r2
#' @param dataset3.r2
#'
#' @param ratio.compare1
#' @param ratio.compare2
#'
#' @param two.side
#'
#' @return
#' @export
#'
#' @import survival
#' @examples
multi_test_recy_pvalue <- function(model1.r1,cov1.r1,dataset1.r1,
                                   model2.r1,cov2.r1,dataset2.r1,
                                   model3.r1,cov3.r1,dataset3.r1,
                                   model1.r2,cov1.r2,dataset1.r2,
                                   model2.r2,cov2.r2,dataset2.r2,
                                   model3.r2,cov3.r2,dataset3.r2,
                                   ratio.compare1 = c(3,3,19),
                                   ratio.compare2 = c(3,3,19),
                                   two.side=TRUE){
  if(two.side==TRUE){
    pp = 0.025
  } else{
    pp = 0.05
  }

  pp1.r1 = pp*ratio.compare1[1]/sum(ratio.compare1)
  pp2.r1 = pp*ratio.compare1[2]/sum(ratio.compare1)
  pp3.r1 = pp*ratio.compare1[3]/sum(ratio.compare1)

  p.value1.r1 = surv_pvalue(survfit(model1.r1~cov1.r1,data = dataset1.r1))$pval
  p.value2.r1 = surv_pvalue(survfit(model2.r1~cov2.r1,data = dataset2.r1))$pval
  p.value3.r1 = surv_pvalue(survfit(model3.r1~cov3.r1,data = dataset3.r1))$pval

  if( p.value1.r1<pp1.r1 && p.value2.r1<pp1.r2){
    pp3.r1 = pp3.r1 + pp1.r1 + pp2.r1
  }else if (p.value1.r1>pp1.r1 && p.value2.r1>pp2.r1){
    pp3.r1 = pp3.r1
  }else if (p.value1.r1<pp1.r1 && p.value2.r1>pp2.r1){
    pp3.r1 = pp3.r1 + pp1.r1
  }else{
    pp3.r1 = pp3.r1 + pp2.r1
  }


  if(p.value3.r1 > pp3.r1){
    return(data.frame(round = c(rep(1,3),rep(2,3)),
                      group = c("PFS in tGE-WT", "PFS in ITT-WT","OS in ITT-WT",
                                "PFS in tGE-WT", "PFS in ITT-WT","OS in ITT-WT"),
                      p.value = c(p.value1.r1, p.value2.r1,p.value3.r1,
                                  NULL,NULL,NULL),
                      p.value.benchmark = c(pp1.r1,pp2.r1,pp3.r1,
                                            NULL, NULL, NULL)))
  }else{
    pp1.r2 = pp3.r1*ratio.compare2[1]/sum(ratio.compare2)
    pp2.r2 = pp3.r1*ratio.compare2[2]/sum(ratio.compare2)
    pp3.r2 = pp3.r1*ratio.compare2[3]/sum(ratio.compare2)

    p.value1.r2 = surv_pvalue(survfit(model1.r2~cov1.r2,data = dataset1.r2))$pval
    p.value2.r2 = surv_pvalue(survfit(model2.r2~cov2.r2,data = dataset2.r2))$pval
    p.value3.r2 = surv_pvalue(survfit(model3.r2~cov3.r2,data = dataset3.r2))$pval

    if( p.value1.r2<pp1.r2 && p.value2.r2 <pp2.r2){
      pp3.r2 = pp3.r2 + pp1.r2 + pp2.r2
    }else if (p.value1.r2>pp1.r2 && p.value2.r2>pp2.r2){
      pp3.r2 <- pp3.r2
    }else if (p.value1.r2<pp1.r2 && p.value2.r2>pp2.r2){
      pp3.r2 <-pp3.r2+pp1.r2
    }else{
      pp3.r2 <- pp3.r2+pp2.r2
    }

    return(data.frame(round = c(rep(1,3), rep(2,3)),
                      group = c("PFS in tGE-WT", "PFS in ITT-WT","OS in ITT-WT",
                                "pfs in tGE-WT", "PFS in ITT-WT","OS in ITT-WT"),
                      p.value = c(p.value1.r1, p.value2.r1,p.value3.r1,
                                  p.value1.r2, p.value2.r2,p.value3.r2),
                      p.value.benchmark = c(pp1.r1,pp2.r1,pp3.r1,
                                            pp1.r2,pp2.r2,pp3.r2) ) )
  }

}





#' Title Drawing Forest Plot for Cox proportional hazards model
#'
#' @param model a Surv(time, status) object
#' @param covariate covariate it can be 1 if all covariates are included, or it can be a single covariate
#' @param dataset dataset it is the dataset we use
#'
#' @return a forest plot
#' @export
#'
#' @import survival
#' @import survminer
#' @examples
forest_plot <- function(modelformula,dataset = veteran){
  fit.coxph <- coxph(modelformula, data = dataset)
  ggforest(fit.coxph, data = dataset)
}