From 123d5c43f10bce0055445df1f8f8bb04e96ff59e Mon Sep 17 00:00:00 2001 From: davidsantiagoquevedo Date: Sat, 2 Nov 2024 19:45:27 +0100 Subject: [PATCH 1/6] feat: added function to compute prediction from Cox model --- R/coh_eff_utils.R | 53 +++++++++++++++++++++++++++++++++++++ man/cox_model_prediction.Rd | 24 +++++++++++++++++ 2 files changed, 77 insertions(+) create mode 100644 man/cox_model_prediction.Rd diff --git a/R/coh_eff_utils.R b/R/coh_eff_utils.R index d58cafe..a2c9b54 100644 --- a/R/coh_eff_utils.R +++ b/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) +} diff --git a/man/cox_model_prediction.Rd b/man/cox_model_prediction.Rd new file mode 100644 index 0000000..d5ae091 --- /dev/null +++ b/man/cox_model_prediction.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/coh_eff_utils.R +\name{cox_model_prediction} +\alias{cox_model_prediction} +\title{Internal function to calculate prediction from Cox model} +\usage{ +cox_model_prediction(cox_model, vaccinated_status, unvaccinated_status) +} +\arguments{ +\item{cox_model}{Result from \code{cox_model} function} +} +\value{ +\code{data.frame} containing +\code{time}: time-to-event until \code{at} +\code{logtime}: log of time-to-event +\code{hazard}: estimated hazard from model +\code{surv}: estimated survival probability from model +\code{loglog}: -log(-log(survival probability)) +\code{strata}: vaccinated/unvaccinated status +} +\description{ +Internal function to calculate prediction from Cox model +} +\keyword{internal} From 03f6108a02decea2ff513f29f6c6d6cddeabd6c3 Mon Sep 17 00:00:00 2001 From: davidsantiagoquevedo Date: Sat, 2 Nov 2024 19:46:35 +0100 Subject: [PATCH 2/6] included prediction of cox model in coh_eff_hr --- R/coh_eff_hr.R | 9 ++++++++- man/coh_eff_hr.Rd | 1 + 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/R/coh_eff_hr.R b/R/coh_eff_hr.R index c2ea1c0..fe2a832 100644 --- a/R/coh_eff_hr.R +++ b/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 ) diff --git a/man/coh_eff_hr.Rd b/man/coh_eff_hr.Rd index d3f8c3d..a4f0337 100644 --- a/man/coh_eff_hr.Rd +++ b/man/coh_eff_hr.Rd @@ -18,6 +18,7 @@ coh_eff_hr( \value{ VE (CI95\%), output from \code{cox_model}, +prediction from \code{cox_model}, and output from \code{km_model} } \description{ From 55a68453c653fcef6d2e5eccfd8e19eafd64ee6c Mon Sep 17 00:00:00 2001 From: davidsantiagoquevedo Date: Sat, 2 Nov 2024 19:47:15 +0100 Subject: [PATCH 3/6] update plot loglog to use cox model prediction --- R/coh_eff_plot.R | 23 ++++++++--------------- R/estimate_vaccineff.R | 2 +- man/plot_loglog.Rd | 9 ++++----- 3 files changed, 13 insertions(+), 21 deletions(-) diff --git a/R/coh_eff_plot.R b/R/coh_eff_plot.R index a194f8b..5baafe4 100644 --- a/R/coh_eff_plot.R +++ b/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", diff --git a/R/estimate_vaccineff.R b/R/estimate_vaccineff.R index baac0a1..5e636bb 100644 --- a/R/estimate_vaccineff.R +++ b/R/estimate_vaccineff.R @@ -171,7 +171,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, diff --git a/man/plot_loglog.Rd b/man/plot_loglog.Rd index 2c08197..e7d8a71 100644 --- a/man/plot_loglog.Rd +++ b/man/plot_loglog.Rd @@ -4,17 +4,16 @@ \alias{plot_loglog} \title{Plot Log-Log Test for Proportional Hazards Hypothesis} \usage{ -plot_loglog(km) +plot_loglog(cox_model_prediction) } \arguments{ -\item{km}{Kaplan-Meier estimation created with \code{km_model}.} +\item{cox_model_prediction}{Prediction from Cox model.} } \value{ Log-log plot. } \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). +This function uses the return from the Cox model +to create a log-log plot. } \keyword{internal} From efdbeed5c0ef0ae25ba7dc120a6869a8314f3ecc Mon Sep 17 00:00:00 2001 From: davidsantiagoquevedo Date: Sat, 2 Nov 2024 19:47:27 +0100 Subject: [PATCH 4/6] update plot loglog to use cox model prediction --- tests/testthat/test-estimate_vaccineff.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/testthat/test-estimate_vaccineff.R b/tests/testthat/test-estimate_vaccineff.R index faeb836..5e114f4 100644 --- a/tests/testthat/test-estimate_vaccineff.R +++ b/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") From 2a93e71047ff9bb43aaefda9ca69e68c56955361 Mon Sep 17 00:00:00 2001 From: davidsantiagoquevedo Date: Thu, 14 Nov 2024 21:52:17 +0100 Subject: [PATCH 5/6] increased threshold for snapsho --- tests/testthat/test-match_cohort.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-match_cohort.R b/tests/testthat/test-match_cohort.R index ee1ce10..ef74e7b 100644 --- a/tests/testthat/test-match_cohort.R +++ b/tests/testthat/test-match_cohort.R @@ -72,6 +72,6 @@ test_that("`match_cohort`: test for input validation", { #### Snapshot for summary test_that("`match_cohort`: summary snapshot", { # nolint for (column in c("balance_all", "balance_match", "summary")) { - expect_snapshot_value(matching[column], style = "json2", tolerance = 1e-2) + expect_snapshot_value(matching[column], style = "json2", tolerance = 1e-1) } }) From d5f1b706e0320126686dedc389bfada69ab55bb1 Mon Sep 17 00:00:00 2001 From: davidsantiagoquevedo Date: Thu, 14 Nov 2024 22:11:40 +0100 Subject: [PATCH 6/6] removed sanpshot test for summary --- tests/testthat/_snaps/estimate_vaccineff.md | 24 --------------------- tests/testthat/test-estimate_vaccineff.R | 16 -------------- 2 files changed, 40 deletions(-) delete mode 100644 tests/testthat/_snaps/estimate_vaccineff.md diff --git a/tests/testthat/_snaps/estimate_vaccineff.md b/tests/testthat/_snaps/estimate_vaccineff.md deleted file mode 100644 index 5f18d08..0000000 --- a/tests/testthat/_snaps/estimate_vaccineff.md +++ /dev/null @@ -1,24 +0,0 @@ -# `estimate_vaccineff`: at not null - - Code - summ - Output - [1] "Vaccine Effectiveness at 90 days computed as VE = 1 - HR:" - [2] " VE lower.95 upper.95" - [3] "1 0.9001 0.2199 0.9872" - [4] "" - [5] "Schoenfeld test for Proportional Hazards assumption:" - [6] "p-value = 0.3187" - -# `summary.vaccineff`: basic expectations - - Code - summ - Output - [1] "Vaccine Effectiveness at 60 days computed as VE = 1 - HR:" - [2] " VE lower.95 upper.95" - [3] "1 0.889 0.1241 0.9859" - [4] "" - [5] "Schoenfeld test for Proportional Hazards assumption:" - [6] "p-value = 0.1904" - diff --git a/tests/testthat/test-estimate_vaccineff.R b/tests/testthat/test-estimate_vaccineff.R index 5e114f4..a40741a 100644 --- a/tests/testthat/test-estimate_vaccineff.R +++ b/tests/testthat/test-estimate_vaccineff.R @@ -42,15 +42,6 @@ test_that("`estimate_vaccineff`: basic expectations", { ) }) -# Truncate at `estimate_vaccineff()` -test_that("`estimate_vaccineff`: at not null", { - - # runs without conditions - 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) @@ -69,13 +60,6 @@ test_that("`estimate_vaccineff`: test for input validation", { ) }) -#### Summary -test_that("`summary.vaccineff`: basic expectations", { - # snapshot for summary - summ <- capture.output(summary.vaccineff(ve)) - expect_snapshot(summ) -}) - #### Plot test_that("`plot.vaccineff`: basic expectations", { # test for loglog plot