Dev

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: drcHelper
 Title: Collection and verification of helper functions for dose-response analysis
-Version: 0.0.4.9000
+Version: 0.0.5
 Authors@R: c(
     person("Zhenglei", "Gao",  email="zhenglei.gao@bayer.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-4042-310X")),

NAMESPACE

@@ -36,6 +36,8 @@ export(calcTaronesTest)
 export(calculate_noec_rstatix)
 export(cochranArmitageTrendTest)
 export(compare_to_control_fisher)
+export(compare_to_control_welch)
+export(compute_mdd_dunnett)
 export(compute_mdd_williams)
 export(contEndpoint)
 export(convert2Score)
@@ -64,6 +66,7 @@ export(monotonicityTest)
 export(mselect.ED)
 export(mselect.ZG)
 export(mselect.plus)
+export(noec_from_trend_test)
 export(oscillating_response)
 export(pavaMean)
 export(plot.modList)
@@ -74,6 +77,7 @@ export(prelimPlot3)
 export(prelimSummary)
 export(prepDataRSCABS)
 export(rankTransform)
+export(report_dunnett_summary)
 export(reshape_drcData)
 export(runRSCABS)
 export(run_RSCA)
@@ -84,6 +88,7 @@ export(simplifyTreatment)
 export(simulate_dose_response)
 export(stepDownRSCABS)
 export(stepDownTrendTestBinom)
+export(stepDownTrendTest_NOEC)
 export(stepKRSCABS)
 export(step_down_RSCABS)
 export(summaryZG)
@@ -177,6 +182,7 @@ importFrom(stats,sd)
 importFrom(stats,setNames)
 importFrom(stats,summary.lm)
 importFrom(stats,symnum)
+importFrom(stats,t.test)
 importFrom(stats,update)
 importFrom(stats,var)
 importFrom(stats,vcov)

NEWS.md

@@ -1,3 +1,22 @@
-# drcHelper (development version)
+# drcHelper 0.0.5
+
+# drcHelper 0.0.4.9001()
 
 * Initial CRAN submission preparation.
+
+## New Features
+- Added `tsk_auto` for trimmed spearman-karber LD50 calculation.
+- Added `SpearmanKarber_modified` for Spearman-Karber LD50 calculation.
+- Added `compare_to_control_welch` for Welch's t-test with multiple comparison adjustment.
+- Added `noec_from_trend_test` as a wrapper for Stepdown Trend test with a NOEC outcome and corresponding printing method.
+- Added `compute_mdd_williams` and `compute_mdd_dunnett` for MDD calculations, 
+  `report_dunnett_summary` for the reporting of test result.
+- Added many articles to explain concepts and verify implementations in `drcHelper`
+
+## Improvments
+- Added test cases validation reports. 
+  - `/inst/SystemTesting/Consolidated_Dunnett_Report.Rmd`
+
+## Bug Fixes
+- Updated `dunnett_test` to correctly handle non-RCBD replicate IDs.
+

R/MDD.R

@@ -55,3 +55,193 @@ compute_mdd_williams <- function(williams_obj, data, formula) {
 
   tibble::tibble(Dose = doses, MDD_pct = MDD_pct)
 }
+
+
+
+
+#' Calculate MDD% for a Dunnett Test Result
+#'
+#' @param dunnett_obj The result object from your dunnett_test function.
+#' @param data The original dataframe used for the test.
+#' @param formula The formula used for the test, e.g., Response ~ Dose.
+#' @return A tibble with Dose, MDD, and MDD_pct.
+#' Calculate MDD% for a Dunnett Test Result (Robust Version 2)
+#'
+#' @param dunnett_obj The result object from your dunnett_test function.
+#' @param alternative The alternative hypothesis ("less", "greater", "two.sided").
+#' @param data The original dataframe used for the test.
+#' @param formula The formula used for the test.
+#' @return A tibble with Dose and MDD_pct.
+#' @export
+compute_mdd_dunnett <- function(dunnett_obj, alternative, data, formula) {
+
+  # 1. Validate input
+  # --- START OF FIX ---
+  # Determine the actual results table from the input object
+  results_table <- NULL
+  if (is.data.frame(dunnett_obj)) {
+    # This handles the direct tibble from broom_dunnett
+    results_table <- dunnett_obj
+  } else if (is.list(dunnett_obj) && "results_table" %in% names(dunnett_obj)) {
+    # This handles the output from the original dunnett_test function
+    results_table <- dunnett_obj$results_table
+  } else {
+    stop("Input 'dunnett_obj' is not a recognized result format.")
+  }
+  # --- END OF FIX ---
+  if (nrow(results_table) == 0) return(tibble::tibble(Dose = numeric(), MDD_pct = numeric()))
+
+  # 2. Extract variable names and ensure dose column is factored
+  resp_name <- all.vars(formula)[1]
+  dose_name <- all.vars(formula)[2]
+  local_data <- data
+  local_data$dose_factor_col <- factor(local_data[[dose_name]], levels = sort(unique(local_data[[dose_name]])))
+  factor_col_name <- "dose_factor_col"
+
+  # 3. Derive the adjusted critical value (Tcrit) using the explicit 'alternative'
+  first_row <- results_table[1, ]
+  Tcrit <- NA
+
+  if (alternative == "less") {
+    Tcrit <- (first_row$conf.high - first_row$estimate) / first_row$std.error
+  } else if (alternative == "greater") {
+    Tcrit <- (first_row$estimate - first_row$conf.low) / first_row$std.error
+  } else if (alternative == "two.sided") {
+    Tcrit <- (first_row$conf.high - first_row$estimate) / first_row$std.error
+  }
+
+  if (is.na(Tcrit) || is.infinite(Tcrit) || Tcrit < 0) {
+    stop("Could not derive a valid Tcrit from confidence intervals for the given alternative.")
+  }
+
+  # 4. Get control mean
+  control_level <- levels(local_data[[factor_col_name]])[1]
+  mu_c <- mean(local_data[local_data[[factor_col_name]] == control_level, resp_name], na.rm = TRUE)
+
+  # 5. Calculate MDD and MDD%
+  SE_diff <- first_row$std.error
+  MDD <- Tcrit * SE_diff
+  MDD_pct <- 100 * MDD / abs(mu_c)
+
+  # 6. Create the final tibble
+  doses <- as.numeric(gsub(" - 0.*", "", results_table$comparison))
+  tibble::tibble(Dose = doses, MDD_pct = rep(MDD_pct, length(doses)))
+}
+
+
+#' Generate a Comprehensive Dunnett Test Summary Report
+#'
+#' This high-level wrapper function performs a Dunnett's test to compare multiple
+#' treatment groups against a single control group. It is designed to be flexible,
+#' allowing for different underlying test functions and handling both standard
+#' fixed-effects (ANOVA) and mixed-effects models. The function combines the
+#' statistical test results with a calculated Minimum Detectable Difference (MDD%)
+#' into a single, formatted output table.
+#'
+#' @param formula A model formula, e.g., `Response ~ Dose` for a fixed-effects model,
+#'   or `Response ~ Dose + (1|Batch)` for a mixed-effects model.
+#' @param data A data frame containing the data specified in the formula.
+#' @param dunnett_test_func A function (passed as a character string or as a function object)
+#'   that performs the Dunnett test and is expected to return an object of class `glht`.
+#'   See `multcomp::glht` for details.
+#' @param alternative A character string specifying the alternative hypothesis. Must be one
+#'   of "two.sided", "less" (default), or "greater".
+#' @param include_random_effect A logical value. Set to `TRUE` if the `formula` includes
+#'   a random effect term (e.g., `(1|Batch)`), which requires a mixed-effects model.
+#'   Defaults to `FALSE`.
+#' @param ... Additional arguments to be passed to the function specified in
+#'   `dunnett_test_func`. A common argument is `alpha` to set the significance
+#'   level (e.g., `alpha = 0.05`).
+#'
+#' @return A data frame containing the comprehensive results for each comparison against
+#'   the control, with columns for the comparison name, estimate difference, test
+#'   statistic, critical value, standard error of the difference, MDD%, and p-value.
+#'
+#' @export
+#' @seealso [multcomp::glht()]
+#'
+#' @examples
+#' \dontrun{
+#' # Generate sample data
+#' set.seed(42)
+#' my_data <- data.frame(
+#'   Dose = factor(rep(c(0, 5, 20, 50), each = 6)),
+#'   Response = c(rnorm(6, 100, 10), rnorm(6, 90, 12),
+#'                rnorm(6, 85, 10), rnorm(6, 80, 11))
+#' )
+#'
+#' # Run the summary function for a one-sided test (less than control)
+#' summary_results <- report_dunnett_summary(
+#'   formula = Response ~ Dose,
+#'   data = my_data,
+#'   dunnett_test_func = "dunnett_test", # Your custom test function
+#'   alternative = "less",
+#'   alpha = 0.05
+#' )
+#'
+#' print(summary_results)
+#' }
+report_dunnett_summary <- function(formula, data, dunnett_test_func, alternative = "less",
+                                   include_random_effect = FALSE,...) {
+
+  resp_name <- all.vars(formula)[1]
+  dose_name <- all.vars(formula)[2]
+
+  # 1. Preliminary Summary
+  prelim_stats <- prelimSummary(data, dose_col = dose_name, response_col = resp_name)
+
+  # 2. Run the Dunnett Test
+  if(dunnett_test_func == "dunnett_test") dunnett_results_obj <- dunnett_test(
+    data = data,
+    response_var = resp_name,
+    dose_var = dose_name,
+    alternative = alternative, # Pass alternative to the test
+    include_random_effect = include_random_effect,
+    ...
+  )
+  if(dunnett_test_func == "broom_dunnett") dunnett_results_obj <- broom_dunnett(formula,
+    data = data,
+    alternative = alternative, # Pass alternative to the test
+    ...
+  )
+
+  # 3. MDD% Calculation
+  mdd_results <- compute_mdd_dunnett(
+    dunnett_obj = dunnett_results_obj,
+    alternative = alternative, # Pass the same alternative to the MDD function
+    data = data,
+    formula = formula
+  )
+  mdd_results$Dose <- factor(mdd_results$Dose)
+  names(mdd_results)[1] <- dose_name
+  # --- START OF FIX ---
+  # Determine the actual results table from the input object
+  results_table <- NULL
+  if (is.data.frame(dunnett_results_obj)) {
+    # This handles the direct tibble from broom_dunnett
+    results_table <- dunnett_results_obj
+  } else if (is.list(dunnett_results_obj) && "results_table" %in% names(dunnett_results_obj)) {
+    # This handles the output from the original dunnett_test function
+    results_table <- dunnett_results_obj$results_table
+  } else {
+    stop("Input 'dunnett_obj' is not a recognized result format.")
+  }
+  # --- END OF FIX ---
+  # 4. Tidy and Join Results
+  test_output <- results_table %>%
+    dplyr::mutate(Dose = as.numeric(gsub(" - 0.*", "", comparison))) %>%
+    dplyr::select(Dose, statistic, p.value, significant)
+  test_output$Dose <- factor(test_output$Dose)
+  names(test_output)[1] <- dose_name
+  final_report <- prelim_stats %>%
+    dplyr::left_join(test_output, by = dose_name) %>%
+    dplyr::left_join(mdd_results, by = dose_name) %>%
+    dplyr::rename(MDD_Dunnett_Percent = MDD_pct, t_statistic = statistic)
+
+  # 5. Reorder and return
+  final_report %>%
+    dplyr::select(
+      Dose, Mean, SD, CV, `% Inhibition`,
+      t_statistic, p.value, significant, MDD_Dunnett_Percent
+    )
+}