From ce6890eaab9556780d04d0bb276df4ad3e9f7332 Mon Sep 17 00:00:00 2001 From: "CINO (Christian Haargaard Olsen)" Date: Thu, 19 Sep 2024 15:33:30 +0200 Subject: [PATCH 1/3] Update definition of BreslowDayFunction Call the function from the DescTools package. --- DESCRIPTION | 3 +- R/BreslowDayFunction.R | 91 ++---------------------------------------- 2 files changed, 6 insertions(+), 88 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index b683f6a..deba68d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -47,7 +47,8 @@ Imports: magrittr, checkmate, cli, - usethis + usethis, + DescTools Remotes: hta-pharma/chef, matthew-phelps/testr diff --git a/R/BreslowDayFunction.R b/R/BreslowDayFunction.R index 63e7ab0..e382e49 100644 --- a/R/BreslowDayFunction.R +++ b/R/BreslowDayFunction.R @@ -7,93 +7,10 @@ #' @return A vector with three values statistic - Breslow and Day test #' statistic pval - p value evtl. based on the Tarone test statistic using a #' \eqn{\chi^2(K-1)} distribution +library(DescTools) +library(testthat) breslowdaytest_ <- function(x, odds_ratio = NA, correct = FALSE) { - # Function to perform the Breslow and Day (1980) test including the - # corrected test by Tarone Uses the equations in Lachin (2000), - # Biostatistical Methods, Wiley, p. 124-125. - # - # Programmed by Michael Hoehle - # Code taken originally from a Biostatistical Methods lecture - # held at the Technical University of Munich in 2008. - # - # Params: - # x - a 2x2xK contingency table - # correct - if TRUE Tarones correction is returned - # - # Returns: - # a vector with three values - # statistic - Breslow and Day test statistic - # pval - p value evtl. based on the Tarone test statistic - # using a \chi^2(K-1) distribution - # - - if (is.na(odds_ratio)) { - #Find the common odds_ratio based on Mantel-Haenszel - oddsratio_hat_mh <- stats::mantelhaen.test(x)$estimate - } else { - oddsratio_hat_mh <- odds_ratio - } - - #Number of strata - k <- dim(x)[3] - #Value of the Statistic - x2_hbd <- 0 - #Value of aj, tildeaj and variance_aj - a <- tildea <- variance_a <- numeric(k) - - for (j in 1:k) { - #Find marginals of table j - mj <- apply(x[, , j], MARGIN = 1, sum) - nj <- apply(x[, , j], MARGIN = 2, sum) - - #Solve for tilde(a)_j - coef <- c(-mj[1] * nj[1] * oddsratio_hat_mh, nj[2] - mj[1] - + (oddsratio_hat_mh * (nj[1] + mj[1])), - 1 - oddsratio_hat_mh) - sols <- Re(polyroot(coef)) - #Take the root, which fulfills 0 < tilde(a)_j <= min(n1_j, m1_j) - tildeaj <- sols[(0 < sols) & (sols <= min(nj[1], mj[1]))] - #Observed value - aj <- x[1, 1, j] - - #Determine other expected cell entries - tildebj <- mj[1] - tildeaj - tildecj <- nj[1] - tildeaj - tildedj <- mj[2] - tildecj - - #Compute \hat{\variance}(a_j | \widehat{\odds_ratio}_MH) - variance_aj <- (1 / tildeaj + 1 / tildebj + 1 / tildecj + 1 / tildedj)^(-1) - - #Compute contribution - x2_hbd <- x2_hbd + as.numeric((aj - tildeaj)^2 / variance_aj) - - #Assign found value for later computations - a[j] <- aj - tildea[j] <- tildeaj - variance_a[j] <- variance_aj - } - - # Compute Tarone corrected test - # Add on 2015: The original equation from the 2008 lecture is incorrect - # as pointed out by Jean-Francois Bouzereau. - x2_hbdt <- as.numeric(x2_hbd - (sum(a) - sum(tildea))^2 / sum(variance_a)) - - dname <- deparse(substitute(x)) - - statistic <- if (correct) x2_hbdt else x2_hbd - parameter <- k - 1 - # Compute p-value based on the Tarone corrected test - pval <- 1 - stats::pchisq(statistic, parameter) - method <- - if (correct) - "Breslow-Day Test on Homogeneity of Odds Ratios (with Tarone correction)" - else - "Breslow-Day Test on Homogeneity of Odds Ratios" - names(statistic) <- "X-squared" - names(parameter) <- "df" - structure(list(statistic = statistic, parameter = parameter, - p.value = pval, method = method, data.name = dname), - class = "htest") - + # Call the BreslowDayTest from the DescTools package + DescTools::BreslowDayTest(x, odds_ratio, correct) } From 96a31b2777f50c0e9df14a95d331092edce9bb6a Mon Sep 17 00:00:00 2001 From: "CINO (Christian Haargaard Olsen)" Date: Thu, 19 Sep 2024 15:45:49 +0200 Subject: [PATCH 2/3] Fix issue with having lirbary call in test --- R/BreslowDayFunction.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/BreslowDayFunction.R b/R/BreslowDayFunction.R index e382e49..1ffb68a 100644 --- a/R/BreslowDayFunction.R +++ b/R/BreslowDayFunction.R @@ -7,8 +7,6 @@ #' @return A vector with three values statistic - Breslow and Day test #' statistic pval - p value evtl. based on the Tarone test statistic using a #' \eqn{\chi^2(K-1)} distribution -library(DescTools) -library(testthat) breslowdaytest_ <- function(x, odds_ratio = NA, correct = FALSE) { # Call the BreslowDayTest from the DescTools package From 70b8e681c038f3f20e6b663c32290f25738e8865 Mon Sep 17 00:00:00 2001 From: Matthew Phelps Date: Fri, 20 Sep 2024 09:34:59 +0200 Subject: [PATCH 3/3] Updates due to updates in dependencies --- NAMESPACE | 1 + R/by_strata_by_trt.R | 6 +++--- R/two_by_two_x.R | 6 +++++- tests/testthat/test-two_by_twos.R | 3 +++ 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 7fb6e64..a4a61b9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -27,3 +27,4 @@ export(sd_value) export(use_chefStats) import(data.table) importFrom(Barnard,barnard.test) +importFrom(magrittr,"%>%") diff --git a/R/by_strata_by_trt.R b/R/by_strata_by_trt.R index 74b691e..ed813d7 100644 --- a/R/by_strata_by_trt.R +++ b/R/by_strata_by_trt.R @@ -247,8 +247,8 @@ demographics_continuous <- function(dat, stat <- dat_cell[, .( mean = mean(get(var), na.rm = TRUE), - median = median(get(var), na.rm = TRUE), - sd = sd(get(var), na.rm = TRUE), + median = stats::median(get(var), na.rm = TRUE), + sd = stats::sd(get(var), na.rm = TRUE), min = min(get(var), na.rm = TRUE), max = max(get(var), na.rm = TRUE), n_non_missing = sum(!is.na(get(var))), @@ -580,7 +580,7 @@ sd_value <- function(dat, unique(by = c(subjectid_var)) stat <- dat_cell[[var]] |> - sd() + stats::sd() return( data.table( diff --git a/R/two_by_two_x.R b/R/two_by_two_x.R index f50f75e..7d0eb9d 100644 --- a/R/two_by_two_x.R +++ b/R/two_by_two_x.R @@ -23,8 +23,9 @@ #' @param treatment_var character. The name of the treatment variable in `dat`. #' @param treatment_refval character. The reference value of the treatment variable in `dat`. #' @param subjectid_var character. Name of the subject identifier variable in `dat` (default is "USUBJID"). -#'@return A matrix +#' @return A matrix #' @export +#' @importFrom magrittr %>% #' make_two_by_two_ <- function(dat, @@ -33,6 +34,8 @@ make_two_by_two_ <- treatment_var, treatment_refval, subjectid_var) { + N <- is_cell <- is_event <- INDEX_ <- treatment <- NULL + dat_ <- copy(dat) n_trt_levels <- dat[, unique(dat, by = treatment_var)][[treatment_var]] |> @@ -103,6 +106,7 @@ make_two_by_two_ <- #' @return data.table of 2x2 table in long format #' @noRd ensure_complete_two_by_two <- function(two_by_two_long, treatment_var) { + is_cell <- N <- NULL n_rows_two_by_two_long <- two_by_two_long[(is_cell), .N, by = c("is_event", treatment_var)] |> NROW() diff --git a/tests/testthat/test-two_by_twos.R b/tests/testthat/test-two_by_twos.R index 991dbf8..3cdd9d6 100644 --- a/tests/testthat/test-two_by_twos.R +++ b/tests/testthat/test-two_by_twos.R @@ -37,6 +37,9 @@ expected <- as.matrix(cbind(yes_counts[,.(N)], no_counts[,.(N)])) rownames(expected) <- no_counts$TRT01A colnames(expected) <- c("outcome_YES", "outcome_NO") +# With data.table 1.16, attributes are preserved in a different way, so strip +# these for testing. +attr(dimnames(expected)[[1]], "label") <- NULL expect_identical(actual, expected) })