From bbe7658f889825ff885991b4dc321c4d2972bc18 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 7 Dec 2023 16:45:06 +0100 Subject: [PATCH] 92: Add unit tests for correlation functions (#98) Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> --- NAMESPACE | 2 - R/corPFSOS.R | 37 ++++--- R/estimateParams.R | 21 +++- _pkgdown.yaml | 2 - man/PFSOSInteg.Rd | 6 +- man/haz.Rd | 8 ++ man/log_p11.Rd | 8 +- man/p11Integ.Rd | 5 +- man/prepareData.Rd | 7 +- tests/testthat/test-corPFSOS.R | 174 ++++++++++++++++++++++++++++++++- 10 files changed, 236 insertions(+), 34 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 74a677ae..a573ccb1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,7 +22,6 @@ export(ExpQuantOS) export(ExpSurvOS) export(ExpSurvPFS) export(PCWInversionMethod) -export(PFSOSInteg) export(PWCsurvOS) export(PWCsurvPFS) export(PwcOSInt) @@ -63,7 +62,6 @@ export(integrateVector) export(logRankTest) export(log_p11) export(negLogLik) -export(p11Integ) export(passedLogRank) export(piecewise_exponential) export(prepareData) diff --git a/R/corPFSOS.R b/R/corPFSOS.R index ec981441..0337a0a6 100644 --- a/R/corPFSOS.R +++ b/R/corPFSOS.R @@ -1,3 +1,5 @@ +# survPFS ---- + #' PFS Survival Function for Different Transition Models #' #' @param transition (`TransitionParameters`)\cr @@ -53,6 +55,8 @@ survPFS.PWCTransition <- function(transition, t) { ) } +# survOS ---- + #' OS Survival Function for Different Transition Models #' #' @param transition (`TransitionParameters`)\cr @@ -113,6 +117,8 @@ survOS.PWCTransition <- function(transition, t) { ) } +# expval ---- + #' Helper Function for Computing E(PFS^2) #' #' @param x (`numeric`)\cr variable of integration. @@ -145,6 +151,8 @@ expvalOSInteg <- function(x, transition) { x * survOS(transition = transition, t = x) } +# p11 ---- + #' Helper Function for `log_p11()` #' #' @param x (`numeric`)\cr variable of integration. @@ -152,21 +160,17 @@ expvalOSInteg <- function(x, transition) { #' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details. #' #' @return Hazard rate at the specified time for the transition from progression to death. -#' @export -#' -#' @examples -#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) -#' p11Integ(2, transition) +#' @keywords internal p11Integ <- function(x, transition) { haz(transition = transition, t = x, trans = 3) } -#' Probability of Remaining in Progression Between Two Time Points for Different Transition Models. +#' Probability of Remaining in Progression Between Two Time Points for Different Transition Models #' #' @param transition (`TransitionParameters`)\cr #' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details. -#' @param s (`numeric`)\cr lower time point. -#' @param t (`numeric`)\cr higher time point. +#' @param s (`numeric`)\cr lower time points. +#' @param t (`numeric`)\cr higher time points. #' @return This returns the natural logarithm of the probability of remaining in progression (state 1) #' between two time points, conditional on being in state 1 at the lower time point. #' @@ -176,6 +180,11 @@ p11Integ <- function(x, transition) { #' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) #' log_p11(transition, 1, 3) log_p11 <- function(transition, s, t) { + assert_numeric(s, finite = TRUE, any.missing = FALSE, lower = 0) + assert_numeric(t, finite = TRUE, any.missing = FALSE, lower = 0) + assert_true(identical(length(s), length(t))) + assert_true(all(t > s)) + intval <- mapply(function(s, t) { stats::integrate(p11Integ, lower = s, @@ -186,6 +195,8 @@ log_p11 <- function(transition, s, t) { -intval } +# PFSOS ---- + #' Helper Function for `survPFSOS()` #' #' @param u (`numeric`)\cr variable of integration. @@ -193,12 +204,10 @@ log_p11 <- function(transition, s, t) { #' @param transition (`TransitionParameters`)\cr #' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details. #' -#' @return Numeric result of the integrand used to calculate the PFS*OS survival function. -#' @export +#' @note Not all vectors `u` and `t` work here due to assertions in [log_p11()]. #' -#' @examples -#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) -#' PFSOSInteg(1, 2, transition) +#' @return Numeric result of the integrand used to calculate the PFS*OS survival function. +#' @keywords internal PFSOSInteg <- function(u, t, transition) { exp(log_p11(transition, u, t / u) + log(survPFS(transition, u)) + log(haz(transition, u, 1))) } @@ -222,6 +231,8 @@ survPFSOS <- function(t, transition) { }) } +# correlation ---- + #' Correlation of PFS and OS event times for Different Transition Models #' #' @param transition (`TransitionParameters`)\cr diff --git a/R/estimateParams.R b/R/estimateParams.R index 72262579..dc1ca3f2 100644 --- a/R/estimateParams.R +++ b/R/estimateParams.R @@ -11,7 +11,10 @@ #' - id (`integer`): patient id. #' - from (`integer`): start event state. #' - to (`integer`): end event state. -#' - trans (`integer`): transition (1, 2 or 3) identifier. +#' - trans (`integer`): transition (1, 2 or 3) identifier +#' - `1`: Transition from state 0 (stable) to 1 (progression). +#' - `2`: Transition from state 0 (stable) to 2 (death). +#' - `3`: Transition from state 1 (progression) to 2 (death). #' - entry (`numeric`): time at which the patient begins to be at risk for the transition. #' - exit (`numeric`): time at which the patient ends to be at risk for the transition. #' - status (`logical`): event indicator for the transition. @@ -76,11 +79,15 @@ prepareData <- function(data) { #' ) #' negLogLik(transition, prepareData(simData)) negLogLik <- function(transition, data) { - with(data, -sum(log(haz(transition, exit, trans)^status - * survTrans(transition, exit, trans) / survTrans(transition, entry, trans)))) + with( + data, + -sum(log( + haz(transition, exit, trans)^status * survTrans(transition, exit, trans) / + survTrans(transition, entry, trans) + )) + ) } - #' Hazard Function for Different Transition Models #' #' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr @@ -88,6 +95,12 @@ negLogLik <- function(transition, data) { #' @param t (`numeric`)\cr time at which hazard is to be computed. #' @param trans (`integer`)\cr index specifying the transition type. #' +#' @details +#' The transition types are: +#' - `1`: Transition from state 0 (stable) to 1 (progression). +#' - `2`: Transition from state 0 (stable) to 2 (death). +#' - `3`: Transition from state 1 (progression) to 2 (death). +#' #' @return The hazard rate for the specified transition and time. #' @export #' diff --git a/_pkgdown.yaml b/_pkgdown.yaml index c927bdc8..b02f9fda 100644 --- a/_pkgdown.yaml +++ b/_pkgdown.yaml @@ -111,8 +111,6 @@ reference: - expvalPFSInteg - expvalOSInteg - log_p11 - - p11Integ - - PFSOSInteg - survPFS - survPFSOS - survOS diff --git a/man/PFSOSInteg.Rd b/man/PFSOSInteg.Rd index 3c1c1d55..288edcc5 100644 --- a/man/PFSOSInteg.Rd +++ b/man/PFSOSInteg.Rd @@ -20,7 +20,7 @@ Numeric result of the integrand used to calculate the PFS*OS survival function. \description{ Helper Function for \code{survPFSOS()} } -\examples{ -transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) -PFSOSInteg(1, 2, transition) +\note{ +Not all vectors \code{u} and \code{t} work here due to assertions in \code{\link[=log_p11]{log_p11()}}. } +\keyword{internal} diff --git a/man/haz.Rd b/man/haz.Rd index e5aacd2f..ed01e992 100644 --- a/man/haz.Rd +++ b/man/haz.Rd @@ -29,6 +29,14 @@ The hazard rate for the specified transition and time. \description{ Hazard Function for Different Transition Models } +\details{ +The transition types are: +\itemize{ +\item \code{1}: Transition from state 0 (stable) to 1 (progression). +\item \code{2}: Transition from state 0 (stable) to 2 (death). +\item \code{3}: Transition from state 1 (progression) to 2 (death). +} +} \section{Methods (by class)}{ \itemize{ \item \code{haz(ExponentialTransition)}: for an exponential transition model. diff --git a/man/log_p11.Rd b/man/log_p11.Rd index d26a399b..36e6fd4f 100644 --- a/man/log_p11.Rd +++ b/man/log_p11.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/corPFSOS.R \name{log_p11} \alias{log_p11} -\title{Probability of Remaining in Progression Between Two Time Points for Different Transition Models.} +\title{Probability of Remaining in Progression Between Two Time Points for Different Transition Models} \usage{ log_p11(transition, s, t) } @@ -10,16 +10,16 @@ log_p11(transition, s, t) \item{transition}{(\code{TransitionParameters})\cr see \code{\link[=exponential_transition]{exponential_transition()}}, \code{\link[=weibull_transition]{weibull_transition()}} or \code{\link[=piecewise_exponential]{piecewise_exponential()}} for details.} -\item{s}{(\code{numeric})\cr lower time point.} +\item{s}{(\code{numeric})\cr lower time points.} -\item{t}{(\code{numeric})\cr higher time point.} +\item{t}{(\code{numeric})\cr higher time points.} } \value{ This returns the natural logarithm of the probability of remaining in progression (state 1) between two time points, conditional on being in state 1 at the lower time point. } \description{ -Probability of Remaining in Progression Between Two Time Points for Different Transition Models. +Probability of Remaining in Progression Between Two Time Points for Different Transition Models } \examples{ transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) diff --git a/man/p11Integ.Rd b/man/p11Integ.Rd index 0240bf71..002a86a9 100644 --- a/man/p11Integ.Rd +++ b/man/p11Integ.Rd @@ -18,7 +18,4 @@ Hazard rate at the specified time for the transition from progression to death. \description{ Helper Function for \code{log_p11()} } -\examples{ -transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) -p11Integ(2, transition) -} +\keyword{internal} diff --git a/man/prepareData.Rd b/man/prepareData.Rd index 47191111..6cd35015 100644 --- a/man/prepareData.Rd +++ b/man/prepareData.Rd @@ -22,7 +22,12 @@ The output data set contains the following columns: \item id (\code{integer}): patient id. \item from (\code{integer}): start event state. \item to (\code{integer}): end event state. -\item trans (\code{integer}): transition (1, 2 or 3) identifier. +\item trans (\code{integer}): transition (1, 2 or 3) identifier +\itemize{ +\item \code{1}: Transition from state 0 (stable) to 1 (progression). +\item \code{2}: Transition from state 0 (stable) to 2 (death). +\item \code{3}: Transition from state 1 (progression) to 2 (death). +} \item entry (\code{numeric}): time at which the patient begins to be at risk for the transition. \item exit (\code{numeric}): time at which the patient ends to be at risk for the transition. \item status (\code{logical}): event indicator for the transition. diff --git a/tests/testthat/test-corPFSOS.R b/tests/testthat/test-corPFSOS.R index 3439b15f..6d07949e 100644 --- a/tests/testthat/test-corPFSOS.R +++ b/tests/testthat/test-corPFSOS.R @@ -1,4 +1,176 @@ -# corPFSOS ---- +# survPFS ---- + +test_that("survPFS works as expected for exponential transition hazards", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + t <- c(0.4, 1.4) + result <- expect_silent(survPFS(transition, t)) + expected <- ExpSurvPFS(t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02) + expect_identical(result, expected) +}) + +test_that("survPFS works as expected for Weibull transition hazards", { + transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3) + t <- c(0.4, 1.4) + result <- expect_silent(survPFS(transition, t)) + expected <- WeibSurvPFS( + t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02, + p01 = transition$weibull_rates$p01, p02 = transition$weibull_rates$p02 + ) + expect_identical(result, expected) +}) + +test_that("survPFS works as expected for piecewise exponential transition hazards", { + transition <- piecewise_exponential( + h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1), + pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9) + ) + t <- c(0.4, 1.4) + result <- expect_silent(survPFS(transition, t)) + expected <- PWCsurvPFS( + t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02, + pw01 = transition$intervals$pw01, pw02 = transition$intervals$pw02 + ) + expect_identical(result, expected) +}) + +# survOS ---- + +test_that("survOS works as expected for exponential transition hazards", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + t <- c(0.4, 1.4) + result <- expect_silent(survOS(transition, t)) + expected <- ExpSurvOS( + t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02, + h12 = transition$hazards$h12 + ) + expect_identical(result, expected) +}) + +test_that("survOS works as expected for Weibull transition hazards", { + transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3) + t <- c(0.4, 1.4) + result <- expect_silent(survOS(transition, t)) + expected <- WeibSurvOS( + t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02, + h12 = transition$hazards$h12, p01 = transition$weibull_rates$p01, + p02 = transition$weibull_rates$p02, p12 = transition$weibull_rates$p12 + ) + expect_identical(result, expected) +}) + +test_that("survOS works as expected for piecewise exponential transition hazards", { + transition <- piecewise_exponential( + h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1), + pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9) + ) + t <- c(0.4, 1.4) + result <- expect_silent(survOS(transition, t)) + expected <- PWCsurvOS( + t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02, + h12 = transition$hazards$h12, pw01 = transition$intervals$pw01, + pw02 = transition$intervals$pw02, pw12 = transition$intervals$pw12 + ) + expect_identical(result, expected) +}) + +# expval ---- + +# expvalPFSInteg ---- + +test_that("expvalPFSInteg works as expected", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + x <- c(0.4, 10.3) + result <- expect_silent(expvalPFSInteg(x, transition)) + expected <- x * survPFS(transition, x) + expect_identical(result, expected) +}) + +# expvalOSInteg ---- + +test_that("expvalOSInteg works as expected", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + x <- c(0.4, 10.3) + result <- expect_silent(expvalOSInteg(x, transition)) + expected <- x * survOS(transition, x) + expect_identical(result, expected) +}) + +# p11 ---- + +test_that("p11Integ works as expected", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + x <- c(1.45, 0.6) + result <- expect_silent(p11Integ(x, transition)) + expected <- haz(transition = transition, t = x, trans = 3) + expect_identical(result, expected) +}) + +test_that("log_p11 works as expected", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + s <- 1 + t <- 3 + result <- expect_silent(log_p11(transition, s, t)) + expected <- -transition$hazards$h12 * (t - s) + expect_equal(result, expected) +}) + +test_that("log_p11 works as expected for multiple time points", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + s <- c(1, 2.2) + t <- c(4.2, 3.2) + result <- expect_silent(log_p11(transition, s, t)) + expected <- -transition$hazards$h12 * (t - s) + expect_equal(result, expected) +}) + +# PFSOS ---- + +test_that("PFSOSInteg works as expected", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + result <- expect_silent(PFSOSInteg(1, 2, transition)) + expected <- 0.0162823 + expect_equal(result, expected, tolerance = 1e-4) +}) + +test_that("PFSOSInteg works as expected for vectors", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + result <- expect_silent(PFSOSInteg(c(1, 2), c(5, 7), transition)) + expected <- c(0.000134, 0.000492) + expect_equal(result, expected, tolerance = 1e-3) +}) + +test_that("survPFSOS works as expected for exponential transition", { + transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6) + result <- expect_silent(survPFSOS(c(0.2, 0.1, 3.5), transition)) + expected <- c(0.38618, 0.52293, 0.01085) + expect_equal(result, expected, tolerance = 1e-4) +}) + +test_that("survPFSOS works as expected for Weibull transitions", { + transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3) + result <- expect_silent(survPFSOS(c(0.2, 0.1, 2.5), transition)) + expected <- c(0.75197, 0.89561, 0.00056) + expect_equal(result, expected, tolerance = 1e-4) +}) + +test_that("survPFSOS works as expected for piecewise constant hazard transitions", { + transition <- piecewise_exponential( + h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1), + pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9) + ) + result <- expect_silent(survPFSOS(c(0.2, 0.1, 2.5), transition)) + expected <- c(0.43021, 0.5602, 0.03678) + expect_equal(result, expected, tolerance = 1e-4) +}) + +# correlation ---- + +test_that("corTrans works as expected", { + transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3) + result <- expect_silent(corTrans(transition)) + expected <- 0.5993243 + expect_equal(result, expected, tolerance = 1e-4) +}) test_that("corPFSOS returns correct central estimate", { transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)