R/KM_prob.R

In Lesly1031/AdjKM.CIF: A package for generating covariate-adjusted Kaplan-Meier and cumulative incidence functions

#' Internal functions. Use R code to calculate adjusted survival probability
#'
#' Use input data, time, status,grouping variables, adjusted covariates,
#' whether to use stratified model, and defining reference group as inputs
#'
#' @noRd
#'
#' @return a dataframe
#'
#'

.survival_prob = function(data,coxout,base_res,strata_group,stratified_cox,reference_group){
  #------------ cumulative baseline hazard for each group -----------
  # make the time roughly the same for coxout and base_res for future merge
  coxout$time = as.numeric(as.character(coxout$time))
  coxout$status = as.numeric(as.character(coxout$status))
  coxout$strata = factor(coxout$strata)
  coxout$time = signif(coxout$time,8)
  base_res$time  = signif(base_res$time,8)

  coxout_event = subset(coxout,coxout$status==1)
  unique_eventTime = data.frame(signif(unique(coxout_event$time),8));names(unique_eventTime)="time"
  stratified_cox=stratified_cox
  if(stratified_cox=="No"){
    if(is.null(reference_group)){
      size = sum(coxout$strata==strata_group)
      numtime = nrow(unique_eventTime) # length of unique event times

      coxout_sub = subset(coxout,coxout$strata==strata_group)
      coxout_sub = coxout_sub[order(coxout_sub$time,coxout_sub$status,decreasing=T),]

      # basehaz for each strata
      # base_res$strata =  factor(base_res$strata)
      # baseres_sub = subset(base_res,base_res$strata==strata_group)## change!! strata_group="TAC/MTX"
      baseres_sub=base_res
      baseres_sub = merge(unique_eventTime,baseres_sub,by="time",all.x=TRUE)
      baseres_sub = baseres_sub[order(baseres_sub$time),]
      # t$hazard= na.locf(t$hazard,na.rm=F)
      # t$hazard = na.locf(t$hazard,fromLast = T,na.rm=F)
      baseres_sub= as.numeric(baseres_sub$hazard)



      adjsurv = matrix(0,numtime,1)

      for(i in 1:size){
        expbz = exp(coxout_sub$linear_predictor[i])
        cumuhaz_res = baseres_sub
        surv = exp(-cumuhaz_res)^expbz # estimate of survival function for group 1

        adjsurv = adjsurv + surv
      }
      adjsurv = adjsurv/size
    } else{
      stop("No reference group for non-stratified cox model")
    }
  } else if(stratified_cox=="Yes"){
    if(is.null(reference_group)){
      stop("Please choose  'G&B', or self-define a reference group for stratified cox model")
    } else if(reference_group=="G&B"){ # Reference = "all", Gail and Byar approach

      size = sum(coxout$strata==strata_group)
      numtime = nrow(unique_eventTime) # length of unique event times

      coxout_sub = subset(coxout,coxout$strata==strata_group)
      coxout_sub = coxout_sub[order(coxout_sub$time,coxout_sub$status,decreasing=T),]

      # basehaz for each strata
      base_res$strata =  factor(base_res$strata)
      baseres_sub = subset(base_res,base_res$strata==strata_group)## change!! strata_group="TAC/MTX"

      baseres_sub = merge(unique_eventTime,baseres_sub,by="time",all.x=TRUE)
      baseres_sub = baseres_sub[order(baseres_sub$time),]
      # t$hazard= na.locf(t$hazard,na.rm=F)
      # t$hazard = na.locf(t$hazard,fromLast = T,na.rm=F)
      baseres_sub= as.numeric(baseres_sub$hazard)



      adjsurv = matrix(0,numtime,1)

      for(i in 1:size){
        expbz = exp(coxout_sub$linear_predictor[i])
        cumuhaz_res = baseres_sub
        surv = exp(-cumuhaz_res)^expbz # estimate of survival function for group 1

        adjsurv = adjsurv + surv
      }
      adjsurv = adjsurv/size

    } else{ ## any user defined group
      reference_var = strsplit(reference_group,split=":")[[1]][1]
      reference_level = strsplit(reference_group,split=":")[[1]][2]
      max_num = nrow(subset(data,data[[reference_var]]==reference_level))
      numtime = nrow(unique_eventTime) # length of unique event times

      coxout_sub = subset(coxout,coxout[["strata"]]==reference_level)
      coxout_sub = coxout_sub[order(coxout_sub$time,coxout_sub$status,decreasing=T),]

      # basehaz for each strata
      base_res$strata =  factor(base_res$strata)
      baseres_sub = subset(base_res,base_res$strata==strata_group)## change!! strata_group="TAC/MTX"

      baseres_sub = merge(unique_eventTime,baseres_sub,by="time",all.x=TRUE)
      baseres_sub = baseres_sub[order(baseres_sub$time),]
      baseres_sub= as.numeric(baseres_sub$hazard)



      adjsurv = matrix(0,numtime,1)

      for(i in 1:max_num){
        expbz = exp(coxout_sub$linear_predictor[i])
        cumuhaz_res = baseres_sub
        surv = exp(-cumuhaz_res)^expbz # estimate of survival function for group 1

        adjsurv = adjsurv + surv
      }
      adjsurv = adjsurv/max_num

    }
  }

  adjsurv
}