diff --git a/DESCRIPTION b/DESCRIPTION index b9afef0..90cac95 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,6 +19,7 @@ Imports: DescTools, dplyr, emmeans, + formatters, lifecycle, lubridate, magrittr, diff --git a/NAMESPACE b/NAMESPACE index de324a1..bf1e170 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,7 +2,9 @@ S3method(print,or_ci) S3method(print,prop_ci) +S3method(print,s_coxph) S3method(print,s_lsmeans) +S3method(print,s_survival) export("%>%") export(Surv) export(derive_bor) @@ -30,6 +32,7 @@ importFrom(dplyr,rowwise) importFrom(dplyr,select) importFrom(dplyr,summarise) importFrom(dplyr,ungroup) +importFrom(formatters,format_value) importFrom(lifecycle,deprecated) importFrom(lubridate,days) importFrom(lubridate,ymd) @@ -45,6 +48,7 @@ importFrom(stats,pbeta) importFrom(stats,quantile) importFrom(stats,rbinom) importFrom(stats,setNames) +importFrom(survival,Surv) importFrom(survival,coxph) importFrom(survival,strata) importFrom(survminer,pairwise_survdiff) diff --git a/R/package.R b/R/package.R index 740262c..e096f7f 100644 --- a/R/package.R +++ b/R/package.R @@ -15,8 +15,13 @@ #' @importFrom survminer pairwise_survdiff #' @importFrom lubridate ymd days #' @importFrom utils combn +#' @importFrom formatters format_value NULL +#' @importFrom survival Surv +#' @export +survival::Surv + utils::globalVariables(c( "ADT", "ADT.x", "ADT.y", "AVAL", "AVALC", "AVALC.x", "AVALC.y", "." )) diff --git a/R/pkg-methods.R b/R/pkg-methods.R index ee4e7ee..aa23358 100644 --- a/R/pkg-methods.R +++ b/R/pkg-methods.R @@ -74,3 +74,203 @@ print.or_ci <- function(x, ...) { invisible(x) } + + +#' @describeIn s_get_survfit prints survival analysis summary from `survfit`. +#' @exportS3Method +#' @keywords internal +print.s_survival <- function(x, ...) { + cat("Surv formula: ", x$params$formula, "\n", sep = "") + cat("Group by: ", paste(unique(x$surv$quantile$group), collapse = ", "), "\n", sep = "") + if (!is.null(x$params$strata)) { + cat("Stratified by: ", paste(x$params$strata, collapse = ", "), "\n", sep = "") + } + cat("Confidence interval type: ", x$params$conf_type, "\n", sep = "") + + cat("\n---\n") + cat("Estimation of Survival Time:\n") + med <- x$surv$median %>% + rowwise() %>% + mutate( + `n(event)` = format_value(c(.data$n, .data$events), format = "xx (xx)"), + Median = format_value(c(.data$median, .data$lower, .data$upper), + format = "xx.xx (xx.xx - xx.xx)" + ) + ) %>% + select("group", "n(event)", "Median") + quant <- x$surv$quantile %>% + rowwise() %>% + mutate( + surv = format_value(c(.data$time, .data$lower, .data$upper), + format = "xx.xx (xx.xx - xx.xx)" + ) + ) %>% + tidyr::pivot_wider( + id_cols = "group", + names_from = "quantile", + names_glue = "Q{quantile}", + values_from = "surv" + ) + rang <- x$surv$range %>% + rowwise() %>% + mutate( + `Range (event)` = format_value(c(.data$event_min, .data$event_max), + format = "(xx.x, xx.x)" + ), + `Range (censor)` = format_value(c(.data$censor_min, .data$censor_max), + format = "(xx.x, xx.x)" + ), + `Range` = format_value(c(.data$min, .data$max), + format = "(xx.x, xx.x)" + ) + ) %>% + select("group", "Range (event)", "Range (censor)", "Range") + list(med, quant, rang) %>% + purrr::reduce(full_join, by = "group") %>% + tidyr::pivot_longer(cols = -1, values_to = "Value", names_to = "Stat") %>% + tidyr::pivot_wider(names_from = "group", values_from = "Value") %>% + tibble::column_to_rownames(var = "Stat") %>% + print() + cat("\n") + + + if (!is.null(x$surv$time_point)) { + cat("---\n") + cat( + "Survival Time at Specified Time Points (", + paste(unique(x$surv$time_point$time), collapse = ","), "):\n", + sep = "" + ) + # print(x$surv$time_point) + x$surv$time_point %>% + rowwise() %>% + mutate( + `Number at risk` = as.character(.data$n.risk), + `Number of event` = as.character(.data$n.event), + `Number of consor` = as.character(.data$n.censor), + `Survival Rate` = format_value(c(.data$surv, .data$lower, .data$upper), + format = "xx.xx (xx.xx - xx.xx)" + ) + ) %>% + select( + "time", "Number at risk", "Number of event", + "Number of consor", "Survival Rate", "group" + ) %>% + ungroup() %>% + dplyr::group_split(.data$time) %>% + purrr::walk(.f = function(dt) { + dt %>% + select(-1) %>% + tidyr::pivot_longer(cols = -c("group"), values_to = "Value", names_to = "Stat") %>% + tidyr::pivot_wider(names_from = "group", values_from = "Value") %>% + tibble::column_to_rownames(var = "Stat") %>% + print() + cat("\n") + }) + } + + if (!is.null(x$surv_diff$rate)) { + cat("---\n") + cat( + "Survival Difference at Specified Time Points (", + paste(unique(x$surv$time_point$time), collapse = ","), "):\n", + sep = "" + ) + x$surv_diff$rate %>% + rowwise() %>% + mutate( + `Diff (Survival Rate)` = format_value(c(.data$surv.diff, .data$lower, .data$upper), + format = "xx.xx (xx.xx - xx.xx)" + ), + `p-value` = format_value(.data$pval, "x.xxxx | (<0.0001)") + ) %>% + select( + "time", "Diff (Survival Rate)", "p-value", "group" + ) %>% + ungroup() %>% + dplyr::group_split(.data$time) %>% + purrr::walk(.f = function(dt) { + dt %>% + select(-1) %>% + tidyr::pivot_longer(cols = -c("group"), values_to = "Value", names_to = "Stat") %>% + tidyr::pivot_wider(names_from = "group", values_from = "Value") %>% + tibble::column_to_rownames(var = "Stat") %>% + print() + cat("\n") + }) + } + + method <- if (x$params$rho == 0) { + "Log-Rank" + } else if (x$params$rho == 1) { + "Peto & Peto" + } + if (!is.null(x$surv_diff$test)) { + cat("---\n") + cat("Hypothesis Testing with ", method, ":\n", sep = "") + x$surv_diff$test %>% + tidyr::pivot_wider(names_from = "comparsion", values_from = "pval") %>% + tibble::column_to_rownames(var = "method") %>% + print() + } + + invisible(x) +} + +#' @describeIn s_get_coxph prints survival analysis summary from `coxph`. +#' @exportS3Method +#' @keywords internal +print.s_coxph <- function(x, ...) { + cat("Surv formula: ", x$params$formula, "\n", sep = "") + cat("Group by: ", paste(x$params$group, collapse = ", "), "\n", sep = "") + if (!is.null(x$params$strata)) { + cat("Stratified by: ", paste(x$params$strata, collapse = ", "), "\n", sep = "") + } + cat("Tie method: ", x$params$ties, "\n", sep = "") + cat("P-value method for HR: ", + switch(x$params$pval_method, + all = paste(c("logtest", "sctest", "waldtest"), collapse = ", "), + log = "logtest", + sc = "sctest", + wald = "waldtest" + ), + "\n", + sep = "" + ) + + cat("\n---\n") + cat("Estimation of Hazard Ratio", + ifelse(!is.null(x$params$strata), " with Stratification", ""), ":\n", + sep = "" + ) + hr <- x$hr %>% + rowwise() %>% + mutate( + `n(event)` = format_value(c(.data$n, .data$events), format = "xx (xx)"), + `Hazard Ratio` = format_value(c(.data$hr, .data$lower, .data$upper), + format = "xx.xx (xx.xx - xx.xx)" + ) + ) %>% + select("comparsion", "n(event)", "Hazard Ratio") + pval <- x$pval %>% + rowwise() %>% + mutate( + pval = format_value(.data$pval, "x.xxxx | (<0.0001)") + ) %>% + select("comparsion", "method", "pval") %>% + tidyr::pivot_wider( + id_cols = "comparsion", + names_from = "method", + names_glue = "p-value ({method})", + values_from = "pval" + ) + + list(hr, pval) %>% + purrr::reduce(full_join, by = "comparsion") %>% + tidyr::pivot_longer(cols = -1, values_to = "Value", names_to = "Stat") %>% + tidyr::pivot_wider(names_from = "comparsion", values_from = "Value") %>% + tibble::column_to_rownames(var = "Stat") %>% + print() + + invisible(x) +} diff --git a/R/summarize-survival.R b/R/summarize-survival.R index d7a0019..2ccc2f0 100644 --- a/R/summarize-survival.R +++ b/R/summarize-survival.R @@ -150,7 +150,8 @@ s_get_survfit <- function(data, # median survival time surv_tb <- if (is.null(km_fit$strata)) { broom::glance(km_fit) %>% - select(n = "nobs", "events", "median", lower = "conf.low", upper = "conf.high") + mutate(group = grps) %>% + select("group", n = "nobs", "events", "median", lower = "conf.low", upper = "conf.high") } else { tibble::as_tibble(summary(km_fit)$table) %>% mutate(group = grps) %>% @@ -218,7 +219,8 @@ s_get_survfit <- function(data, survdiff <- h_pairwise_survdiff( formula = formula, strata = strata, - data = data + data = data, + rho = rho ) tibble::tibble( comparsion = paste(bylist[, 1], bylist[, 2], sep = " vs. "), @@ -320,6 +322,8 @@ s_get_survfit <- function(data, #' ) #' s_get_coxph(data = dat, formula = Surv(LENFOL, FSTAT) ~ AFB) #' +#' s_get_coxph(data = dat, formula = Surv(LENFOL, FSTAT) ~ AFB, pval_method = "sc") +#' #' # specify the stratified log-rank test #' s_get_coxph( #' data = dat, @@ -425,6 +429,7 @@ s_get_coxph <- function(data, pval = pval_tb, params = list( formula = format(formula), + group = paste(group_var, grps, sep = "="), ties = ties, conf_level = conf_level, strata = strata, diff --git a/R/utils.R b/R/utils.R index 0457c17..c835cf9 100644 --- a/R/utils.R +++ b/R/utils.R @@ -13,20 +13,8 @@ #' @return The result of calling `rhs(lhs)`. NULL -#' Create a Survival Object -#' -#' A copy from [survival::Surv] in `survival` package. -#' -#' @inheritDotParams survival::Surv -#' -#' @return An object of class `Surv`. -#' @export -Surv <- function(...) { - survival::Surv(...) -} - #' @describeIn s_get_survfit a modified version of [survminer::pairwise_survdiff()] -#' that can calculate the pairwise comparisons of survival curves regardless of +#' that can calculates the pairwise comparisons of survival curves regardless of #' whether stratification is given. #' #' @export diff --git a/man/Surv.Rd b/man/Surv.Rd deleted file mode 100644 index 871472e..0000000 --- a/man/Surv.Rd +++ /dev/null @@ -1,52 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{Surv} -\alias{Surv} -\title{Create a Survival Object} -\usage{ -Surv(...) -} -\arguments{ -\item{...}{ - Arguments passed on to \code{\link[survival:Surv]{survival::Surv}} - \describe{ - \item{\code{time}}{ - for right censored data, this is the follow up time. For interval - data, the first argument is the starting time for the interval. - } - \item{\code{event}}{ - The status indicator, normally 0=alive, 1=dead. Other choices are - \code{TRUE}/\code{FALSE} (\code{TRUE} = death) or 1/2 (2=death). For - interval censored data, the status indicator is 0=right censored, - 1=event at \code{time}, 2=left censored, 3=interval censored. - For multiple endpoint data the event variable will be a factor, - whose first level is treated as censoring. - Although unusual, the event indicator can be omitted, in which case - all subjects are assumed to have an event. - } - \item{\code{time2}}{ - ending time of the interval for interval censored or counting - process data only. Intervals are assumed to be open on the left and - closed on the right, \code{(start, end]}. For counting process - data, \code{event} indicates whether an event occurred at the end of - the interval. - } - \item{\code{type}}{ - character string specifying the type of censoring. Possible values - are \code{"right"}, \code{"left"}, \code{"counting"}, - \code{"interval"}, \code{"interval2"} or \code{"mstate"}. - } - \item{\code{origin}}{ - for counting process data, the hazard function origin. This option - was intended to be used in conjunction with a model containing - time dependent - strata in order to align the subjects properly when they cross over - from one strata to another, but it has rarely proven useful.} - }} -} -\value{ -An object of class \code{Surv}. -} -\description{ -A copy from \link[survival:Surv]{survival::Surv} in \code{survival} package. -} diff --git a/man/reexports.Rd b/man/reexports.Rd new file mode 100644 index 0000000..cdc12e9 --- /dev/null +++ b/man/reexports.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/package.R +\docType{import} +\name{reexports} +\alias{reexports} +\alias{Surv} +\title{Objects exported from other packages} +\keyword{internal} +\description{ +These objects are imported from other packages. Follow the links +below to see their documentation. + +\describe{ + \item{survival}{\code{\link[survival]{Surv}}} +}} + diff --git a/man/s_get_coxph.Rd b/man/s_get_coxph.Rd index 8ad0ef4..480984a 100644 --- a/man/s_get_coxph.Rd +++ b/man/s_get_coxph.Rd @@ -1,9 +1,12 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/summarize-survival.R -\name{s_get_coxph} +% Please edit documentation in R/pkg-methods.R, R/summarize-survival.R +\name{print.s_coxph} +\alias{print.s_coxph} \alias{s_get_coxph} \title{Summarize Cox Proportional Hazards Regression Model} \usage{ +\method{print}{s_coxph}(x, ...) + s_get_coxph( data, formula, @@ -16,6 +19,8 @@ s_get_coxph( ) } \arguments{ +\item{...}{other arguments to be passed to \code{\link[survival:coxph]{survival::coxph()}}.} + \item{data}{(\code{data.frame})\cr a data frame as input.} \item{formula}{(\code{formula})\cr a formula with survival object.} @@ -33,8 +38,6 @@ default is to present all three methods (Likelihood ratio test, Wald test and Score (logrank) test).} \item{pairwise}{(\code{logical})\cr whether to conduct the pairwise comparison.} - -\item{...}{other arguments to be passed to \code{\link[survival:coxph]{survival::coxph()}}.} } \value{ An object of class \code{s_coxph} is a list contains hazards ratio and p-value tables. @@ -45,6 +48,11 @@ An object of class \code{s_coxph} is a list contains hazards ratio and p-value t This function summarizes \code{coxph} model with or without stratification, and provides hazards ratio, confidence interval and p-value for common clinical survival analysis. } +\section{Functions}{ +\itemize{ +\item \code{print(s_coxph)}: prints survival analysis summary from \code{coxph}. + +}} \examples{ data("whas500") @@ -55,6 +63,8 @@ dat <- whas500 \%>\% ) s_get_coxph(data = dat, formula = Surv(LENFOL, FSTAT) ~ AFB) +s_get_coxph(data = dat, formula = Surv(LENFOL, FSTAT) ~ AFB, pval_method = "sc") + # specify the stratified log-rank test s_get_coxph( data = dat, @@ -75,3 +85,4 @@ dat2 <- whas500 \%>\% ) s_get_coxph(data = dat2, formula = Surv(LENFOL, FSTAT) ~ AFB) } +\keyword{internal} diff --git a/man/s_get_survfit.Rd b/man/s_get_survfit.Rd index 394c619..c105725 100644 --- a/man/s_get_survfit.Rd +++ b/man/s_get_survfit.Rd @@ -1,10 +1,14 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/summarize-survival.R, R/utils.R -\name{s_get_survfit} +% Please edit documentation in R/pkg-methods.R, R/summarize-survival.R, +% R/utils.R +\name{print.s_survival} +\alias{print.s_survival} \alias{s_get_survfit} \alias{h_pairwise_survdiff} \title{Summarize Survival Curves Analysis} \usage{ +\method{print}{s_survival}(x, ...) + s_get_survfit( data, formula, @@ -21,6 +25,8 @@ s_get_survfit( h_pairwise_survdiff(formula, data, strata = NULL, rho = 0) } \arguments{ +\item{...}{other arguments to be passed to \code{\link[survival:survfit]{survival::survfit()}}.} + \item{data}{(\code{data.frame})\cr a data frame as input.} \item{formula}{(\code{formula})\cr a formula of survival model with survival object. @@ -45,8 +51,6 @@ of comparing survival curves. Default is 0 that is log-rank , others options can see \code{\link[survival:survdiff]{survival::survdiff()}}.} \item{pairwise}{(\code{logical})\cr whether to conduct the pairwise comparison.} - -\item{...}{other arguments to be passed to \code{\link[survival:survfit]{survival::survfit()}}.} } \value{ An object of class \code{s_survival} is a list contains several summary tables @@ -68,8 +72,10 @@ time-point survival rates and difference of survival rates. } \section{Functions}{ \itemize{ +\item \code{print(s_survival)}: prints survival analysis summary from \code{survfit}. + \item \code{h_pairwise_survdiff()}: a modified version of \code{\link[survminer:pairwise_survdiff]{survminer::pairwise_survdiff()}} -that can calculate the pairwise comparisons of survival curves regardless of +that can calculates the pairwise comparisons of survival curves regardless of whether stratification is given. }} @@ -131,3 +137,4 @@ s_get_survfit( pairwise = TRUE ) } +\keyword{internal} diff --git a/tests/spelling.R b/tests/spelling.R index 6713838..13f77d9 100644 --- a/tests/spelling.R +++ b/tests/spelling.R @@ -1,3 +1,6 @@ -if(requireNamespace('spelling', quietly = TRUE)) - spelling::spell_check_test(vignettes = TRUE, error = FALSE, - skip_on_cran = TRUE) +if (requireNamespace("spelling", quietly = TRUE)) { + spelling::spell_check_test( + vignettes = TRUE, error = FALSE, + skip_on_cran = TRUE + ) +} diff --git a/tests/testthat/test-pkg-methods.R b/tests/testthat/test-pkg-methods.R index ff69bbe..6310c1b 100644 --- a/tests/testthat/test-pkg-methods.R +++ b/tests/testthat/test-pkg-methods.R @@ -18,3 +18,63 @@ test_that("print.s_lsmeans works as expected for mmrm", { res3 <- capture_output(print(s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = -2, alternative = "greater"))) expect_match(res3, "Null hypothesis is θ inferiority to -2", fixed = TRUE) }) + +test_that("print.s_survival works as expected", { + data("whas500") + dat <- whas500 %>% + mutate( + AFB = factor(AFB, levels = c(1, 0)) + ) + res <- capture_output(print(s_get_survfit(data = dat, formula = Surv(LENFOL, FSTAT) ~ AFB))) + expect_match(res, "Surv formula: Surv(LENFOL, FSTAT) ~ AFB", fixed = TRUE) + expect_match(res, "Group by: AFB=1, AFB=0", fixed = TRUE) + expect_match(res, "Confidence interval type: log-log", fixed = TRUE) + expect_match(res, "Estimation of Survival Time:", fixed = TRUE) + expect_match(res, "Hypothesis Testing with Log-Rank", fixed = TRUE) + + res2 <- capture_output(print(s_get_survfit(data = dat, formula = Surv(LENFOL, FSTAT) ~ 1))) + expect_match(res2, "Group by: Total", fixed = TRUE) + + res3 <- capture_output(print(s_get_survfit( + data = dat, + formula = Surv(LENFOL, FSTAT) ~ AFB, + time_point = c(12, 36, 60) + ))) + expect_match(res3, "Survival Time at Specified Time Points (12,36,60)", fixed = TRUE) + expect_match(res3, "Survival Difference at Specified Time Points (12,36,60)", fixed = TRUE) + + res4 <- capture_output(print(s_get_survfit( + data = dat, + formula = Surv(LENFOL, FSTAT) ~ AFB, + strata = c("AGE", "GENDER") + ))) + expect_match(res4, "Stratified by: AGE, GENDER", fixed = TRUE) +}) + +test_that("print.s_coxph works as expected", { + data("whas500") + dat <- whas500 %>% + mutate( + AFB = factor(AFB, levels = c(1, 0)) + ) + res <- capture_output(print(s_get_coxph(data = dat, formula = Surv(LENFOL, FSTAT) ~ AFB))) + expect_match(res, "Surv formula: Surv(LENFOL, FSTAT) ~ AFB", fixed = TRUE) + expect_match(res, "Group by: AFB=1, AFB=0", fixed = TRUE) + expect_match(res, "Tie method: efro", fixed = TRUE) + expect_match(res, "P-value method for HR: logtest, sctest, waldtest", fixed = TRUE) + expect_match(res, "Estimation of Hazard Ratio", fixed = TRUE) + + res2 <- capture_output(print(s_get_coxph( + data = dat, + formula = Surv(LENFOL, FSTAT) ~ AFB, + strata = c("AGE", "GENDER") + ))) + expect_match(res2, "Stratified by: AGE, GENDER", fixed = TRUE) + + res3 <- capture_output(print(s_get_coxph( + data = dat, + formula = Surv(LENFOL, FSTAT) ~ AFB, + pval_method = "log" + ))) + expect_match(res3, "P-value method for HR: logtest", fixed = TRUE) +}) diff --git a/tests/testthat/test-summarize-survival.R b/tests/testthat/test-summarize-survival.R index 8b7992a..89a1a89 100644 --- a/tests/testthat/test-summarize-survival.R +++ b/tests/testthat/test-summarize-survival.R @@ -12,6 +12,7 @@ test_that("s_get_survfit works as expected with default arguments in single grou expect_equal( res$surv$median, tibble::tibble( + group = "Total", n = 500, events = 215, median = 53.4538,