Merge pull request #20 from hta-pharma/feature/update_breslowday

Update definition of BreslowDayFunction

DESCRIPTION

@@ -47,7 +47,8 @@ Imports:
     magrittr,
     checkmate,
     cli,
-    usethis
+    usethis,
+	DescTools
 Remotes:
     hta-pharma/chef,
     matthew-phelps/testr

NAMESPACE

@@ -27,3 +27,4 @@ export(sd_value)
 export(use_chefStats)
 import(data.table)
 importFrom(Barnard,barnard.test)
+importFrom(magrittr,"%>%")

R/BreslowDayFunction.R

@@ -9,91 +9,6 @@
 #'    \eqn{\chi^2(K-1)} distribution
 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 <http://www.math.su.se/~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)
 }

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(

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()
 

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)
 })