Create loglog plot using prediction from Cox model

R/coh_eff_hr.R

@@ -12,6 +12,7 @@
 #' @inheritParams estimate_vaccineff
 #' @return VE (CI95%),
 #' output from `cox_model`,
+#' prediction from `cox_model`,
 #' and output from `km_model`
 #' @keywords internal
 
@@ -25,7 +26,7 @@ coh_eff_hr <- function(data_set,
                        start_cohort,
                        end_cohort) {
 
-  # Kaplan-Meier model for loglog curve
+  # Kaplan-Meier model for survival curve
   km <- km_model(data_set = data_set,
     outcome_status_col = outcome_status_col,
     time_to_event_col = time_to_event_col,
@@ -45,6 +46,11 @@ coh_eff_hr <- function(data_set,
     unvaccinated_status = unvaccinated_status
   )
 
+  cx_prediction <- cox_model_prediction(cox_model = cx,
+    vaccinated_status = vaccinated_status,
+    unvaccinated_status = unvaccinated_status
+  )
+
   # Vaccine effectiveness = 1 - HR
   eff <- data.frame(
     VE = 1 - cx$hr,
@@ -57,6 +63,7 @@ coh_eff_hr <- function(data_set,
   ve <- list(
     ve = eff,
     cox_model = cx,
+    cox_model_prediction = cx_prediction,
     kaplan_meier = km
   )
 

R/coh_eff_plot.R

@@ -1,34 +1,27 @@
 #' @title Plot Log-Log Test for Proportional Hazards Hypothesis
 #'
-#' @description This function uses the return from the Kaplan-Meier model
-#' to create a log-log plot. It calculates log(-log(Survival Probability))
-#' and log(Time).
+#' @description This function uses the return from the Cox model
+#' to create a log-log plot.
 #' @importFrom rlang .data
-#' @param km Kaplan-Meier estimation created with `km_model`.
+#' @param cox_model_prediction Prediction from Cox model.
 #' @return Log-log plot.
 #' @keywords internal
 
-plot_loglog <- function(km) {
-  # Calculate log-variables
-  km$loglog <- log(-log(km$surv))
-  km$logtime <- log(km$time)
-
+plot_loglog <- function(cox_model_prediction) {
   # strata levels were defined as order actors to extract like this
-  vaccinated_status <- levels(km$strata)[1]
-  unvaccinated_status <- levels(km$strata)[2]
+  vaccinated_status <- levels(cox_model_prediction$strata)[1]
+  unvaccinated_status <- levels(cox_model_prediction$strata)[2]
   # Plot colors
   vaccinated_color <- "steelblue"
   unvaccinated_color <- "darkred"
-
-  # Plot
-  plt <- ggplot2::ggplot(data = km) +
+  plt <- ggplot2::ggplot(data = cox_model_prediction) +
     ggplot2::geom_step(ggplot2::aes(x = .data$logtime,
                                     y = .data$loglog,
                                     color = .data$strata)
     ) +
     ggplot2::theme_classic() +
     ggplot2::labs(x = "Log[Time to event] (Days)",
-                  y = "Log[-Log[Surv.]]") +
+                  y = "-Log[-Log[Surv.]]") +
     ggplot2::labs(colour = "Vaccine Status") +
     ggplot2::scale_color_manual(
       name = "Vaccine Status",

R/coh_eff_utils.R

@@ -4,6 +4,7 @@
 #' @param model `{survival}` object containing the model
 #' @return `data.frame` with survival data
 #' @keywords internal
+
 extract_surv_model <- function(model, start_cohort, end_cohort) {
   days <- end_cohort - start_cohort
   tte <- seq(0, as.numeric(days) - 1, by = 1)
@@ -23,6 +24,7 @@ extract_surv_model <- function(model, start_cohort, end_cohort) {
 #' "surv", "lower", "upper",
 #' "cumincidence", "cumincidence_lower", "cumincidence_upper"
 #' @keywords internal
+
 km_model <- function(data_set,
                      outcome_status_col,
                      time_to_event_col,
@@ -78,6 +80,7 @@ km_model <- function(data_set,
 #' `{survival}` object with model
 #' `{survival}` object with Schoenfeld test
 #' @keywords internal
+
 cox_model <- function(data_set,
                       outcome_status_col,
                       time_to_event_col,
@@ -125,3 +128,53 @@ cox_model <- function(data_set,
 
   return(cx)
 }
+
+#' @title Internal function to calculate prediction from Cox model
+#'
+#' @inheritParams estimate_vaccineff
+#' @param cox_model Result from `cox_model` function
+#' @return `data.frame` containing
+#' `time`: time-to-event until `at`
+#' `logtime`: log of time-to-event
+#' `hazard`: estimated hazard from model
+#' `surv`: estimated survival probability from model
+#' `loglog`: -log(-log(survival probability))
+#' `strata`: vaccinated/unvaccinated status
+#' @keywords internal
+
+cox_model_prediction <- function(cox_model,
+                                 vaccinated_status,
+                                 unvaccinated_status) {
+  bh <- survival::basehaz(cox_model$model)
+  coef <- stats::coef(cox_model$model)
+  surv_v <- exp(-bh[, 1])^(exp(coef))
+  surv_u <- exp(-bh[, 1])
+
+  loglog_v <- -log(-log(surv_v))
+  loglog_u <- -log(-log(surv_u))
+
+  predicted_u <- data.frame(
+    time = bh[, 2],
+    logtime = log(bh[, 2]),
+    hazard = bh[, 1],
+    surv = surv_u,
+    loglog = loglog_u
+  )
+
+  predicted_v <- data.frame(
+    time = bh[, 2],
+    logtime = log(bh[, 2]),
+    hazard = bh[, 1],
+    surv = surv_v,
+    loglog = loglog_v
+  )
+  predicted_v$strata <- vaccinated_status
+  predicted_u$strata <- unvaccinated_status
+
+  predicted <- rbind(predicted_u, predicted_v)
+  levels(predicted$strata) <- factor(
+    c(vaccinated_status, unvaccinated_status),
+    ordered = TRUE
+  )
+  return(predicted)
+}

R/estimate_vaccineff.R

@@ -178,7 +178,7 @@ plot.vaccineff <- function(x,
   )
 
   if (type == "loglog") {
-    plt <- plot_loglog(x$kaplan_meier)
+    plt <- plot_loglog(x$cox_model_prediction)
   } else if (type == "surv") {
     plt <- plot_survival(x$kaplan_meier,
       percentage = percentage,

tests/testthat/test-estimate_vaccineff.R

@@ -32,7 +32,7 @@ test_that("`estimate_vaccineff`: basic expectations", {
   )
 
   # returns `vaccineff` s3class object
-  ve <-estimate_vaccineff(vaccineff_data, at = 60)
+  ve <- estimate_vaccineff(vaccineff_data, at = 60)
   expect_s3_class(
     ve, "vaccineff"
   )
@@ -46,13 +46,13 @@ test_that("`estimate_vaccineff`: basic expectations", {
 test_that("`estimate_vaccineff`: at not null", {
 
   # runs without conditions
-  ve <-estimate_vaccineff(vaccineff_data, at = 90)
+  ve <- estimate_vaccineff(vaccineff_data, at = 90)
   summ <- capture.output(summary.vaccineff(ve))
   expect_snapshot(summ)
 })
 
 #### Tests for generic methods plot and summary ####
-ve <-estimate_vaccineff(vaccineff_data, at = 60)
+ve <- estimate_vaccineff(vaccineff_data, at = 60)
 
 # test for input validation of plot() and summary() methods
 test_that("`estimate_vaccineff`: test for input validation", {
@@ -80,7 +80,7 @@ test_that("`summary.vaccineff`: basic expectations", {
 test_that("`plot.vaccineff`: basic expectations", {
   # test for loglog plot
   plt <- plot.vaccineff(ve, type = "loglog")
-  expect_identical(plt$labels$y, "Log[-Log[Surv.]]")
+  expect_identical(plt$labels$y, "-Log[-Log[Surv.]]")
   plt <- plot.vaccineff(ve, type = "surv")
   expect_identical(plt$labels$y, "Survival probability")
   expect_s3_class(plt, "ggplot")