R/HazardRatioCoxStat.R

In gpaux/Mediana: Clinical Trial Simulations

######################################################################################################################

# Compute the hazard ratio based on non-missing values in the combined sample

HazardRatioCoxStat = function(sample.list, parameter) {
  
  # Determine the function call, either to generate the statistic or to return description
  call = (parameter[[1]] == "Description")
  
  if (call == FALSE | is.na(call)) {
    
    # Error checks
    if (length(sample.list)!=2)
      stop("Analysis model: Two samples must be specified in the HazardRatioCoxStat statistic.")
    
    # Outcomes in Sample 1
    outcome1 = sample.list[[1]][, "outcome"]
    # Remove the missing values due to dropouts/incomplete observations
    outcome1.complete = outcome1[stats::complete.cases(outcome1)]
    # Observed events in Sample 1 (negation of censoring indicators)
    event1 = !sample.list[[1]][, "patient.censor.indicator"]
    event1.complete = event1[stats::complete.cases(outcome1)]
    # Sample size in Sample 1
    n1 = length(outcome1.complete)
    
    # Outcomes in Sample 2
    outcome2 = sample.list[[2]][, "outcome"]
    # Remove the missing values due to dropouts/incomplete observations
    outcome2.complete = outcome2[stats::complete.cases(outcome2)]
    # Observed events in Sample 2 (negation of censoring indicators)
    event2 = !sample.list[[2]][, "patient.censor.indicator"]
    event2.complete = event2[stats::complete.cases(outcome2)]
    # Sample size in Sample 2
    n2 = length(outcome2.complete)
    
    # Create combined samples of outcomes, censoring indicators (all events are observed) and treatment indicators
    outcome = c(outcome1.complete, outcome2.complete)
    event = c(event1.complete, event2.complete)
    treatment = c(rep(0, n1), rep(1, n2))
    
    # Get the HR from the Cox method
    result = summary(survival::coxph(survival::Surv(outcome, event) ~ treatment))$coef[,"exp(coef)"]
  }
  
  else if (call == TRUE) {
    result = list("Hazard Ratio (Cox)")
  }
  
  return(result)
}
# End of HazardRatioCoxStat