From fcd69c89cf38e9e1f29cb90ff41656890d26219b Mon Sep 17 00:00:00 2001 From: Zhenglei Gao Date: Thu, 2 Oct 2025 18:01:52 +0200 Subject: [PATCH 1/5] updated MDD --- R/MDD.R | 146 +++++ R/broom.R | 114 ++-- R/data_description.R | 4 +- R/quantal.R | 3 +- inst/DevNotes/function_dev/dev_dunnett.R | 519 ------------------ man/DixonQ.Rd | 4 +- man/Tarone.test.Rd | 3 +- man/broom_dunnett.Rd | 2 + man/compute_mdd_dunnett.Rd | 26 + man/report_dunnett_summary.Rd | 32 ++ ...t_for_Data_with_Hierarchical_Structure.Rmd | 18 + .../articles/MDD-in-Regulatory-Context.Rmd | 77 +-- .../Verification-of-Williams-Test.Rmd | 11 +- vignettes/knitr-setup.R | 4 +- 14 files changed, 334 insertions(+), 629 deletions(-) delete mode 100644 inst/DevNotes/function_dev/dev_dunnett.R create mode 100644 man/compute_mdd_dunnett.Rd create mode 100644 man/report_dunnett_summary.Rd diff --git a/R/MDD.R b/R/MDD.R index 150467f..ab31ba0 100644 --- a/R/MDD.R +++ b/R/MDD.R @@ -55,3 +55,149 @@ 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. + +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 (Robust Version) +#' +#' @param formula The model formula, e.g., Response ~ Dose. +#' @param data The dataframe containing the data. +#' @param dunnett_test_func The dunnett_test function you use. +#' @param alternative The alternative hypothesis to be used for the test. +#' @param ... Additional arguments passed to your dunnett_test_func. +#' @return A single tibble combining descriptive stats, test results, and MDD%. +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 + ) +} diff --git a/R/broom.R b/R/broom.R index c105df0..27ea59e 100644 --- a/R/broom.R +++ b/R/broom.R @@ -125,14 +125,21 @@ broom_williams <- function(x, method = c("Williams_PMCMRplus", "Williams_JG"), . direction <- "increasing" } - - wt <- drcHelper::williamsTest_JG(df=model_data,resp=resp_name, trt = factor_name, direction = direction) + ## Browse[1]> factor_name + ## [1] "factor(Dose)" + ## This is an edge case where the factor_name is not really included in the model_data frame. + # --- PATCH --- + # Clean the factor name to handle cases like "factor(Dose)" + # This removes "factor(" and the closing ")" to get the clean variable name. + clean_factor_name <- gsub("factor\\(|\\)", "", factor_name) + # --- END PATCH --- + wt <- drcHelper::williamsTest_JG(df=model_data,resp=resp_name, trt = clean_factor_name, direction = direction) if(is.null(list(...)$direction) || list(...)$direction == "decreasing"){ - comparisons <- paste(wt$dose,"- 0 >= 0") + comparisons <- paste( wt[,clean_factor_name],"- 0 >= 0") }else{ - comparisons <- paste(wt$dose,"- 0 <= 0") + comparisons <- paste( wt[,clean_factor_name],"- 0 <= 0") } wt <- wt[rev(seq_len(nrow(wt))), , drop = FALSE] @@ -196,56 +203,56 @@ broom_williams <- function(x, method = c("Williams_PMCMRplus", "Williams_JG"), . #' @importFrom PMCMRplus dunnettTest #' @importFrom DescTools DunnettTest #' @export -broom_dunnett <- function(x, method = c("Dunnett_multcomp", "Dunnett_DescTools", "Dunnett_PMCMRplus"), ...) { +broom_dunnett <- function(x, method = c("Dunnett_multcomp", "Dunnett_DescTools", "Dunnett_PMCMRplus"), + alternative = c("two.sided", "less", "greater"), alpha=0.05,...) { method <- match.arg(method) + alternative <- match.arg(alternative) + # ... (Dunnett_PMCMRplus and Dunnett_DescTools blocks are unchanged) ... if (method == "Dunnett_PMCMRplus") { result <- try({ - dt <- PMCMRplus::dunnettTest(x, ...) - - # Extract summary information - dt_summary <- summaryZG_dunnett(dt) - comparisons <- rownames(dt_summary) - - # Create standardized output + dt <- PMCMRplus::dunnettTest(x, alternative = alternative, ...) + control_name <- colnames(dt$p.value)[1] + treatment_names <- rownames(dt$p.value) + if (alternative == "less") { + comparisons <- paste(treatment_names, "-", control_name, ">=", "0") + } else if (alternative == "greater") { + comparisons <- paste(treatment_names, "-", control_name, "<=", "0") + } else { + comparisons <- paste(treatment_names, "-", control_name, "==", "0") + } tibble::tibble( comparison = comparisons, - estimate = NA, - p.value = dt_summary$`Pr(>|t|)`, - conf.low = NA, - conf.high = NA, + estimate = NA_real_, + statistic = as.vector(dt$statistic), + p.value = as.vector(dt$p.value), + std.error = NA_real_, + conf.low = NA_real_, + conf.high = NA_real_, method = "Dunnett_PMCMRplus" ) }, silent = TRUE) } else if (method == "Dunnett_DescTools") { result <- try({ - dt <- DescTools::DunnettTest(x, ...) - - # Extract the first element if it's a list (DescTools returns a list for multiple controls) - if (is.list(dt) && !is.data.frame(dt)) { - dt <- dt[[1]] - } - - # Convert to tibble with standardized column names + if(alternative != "two.sided") stop("DescTools package does not accept one-sided test, try other methods instead. We recommend multcomp as the standard approach. ") + dt_list <- DescTools::DunnettTest(x, alternative = alternative, ...) + dt <- dt_list[[1]] tibble::tibble( comparison = rownames(dt), estimate = dt[, "diff"], + statistic = NA_real_, p.value = dt[, "pval"], - std.error = NA, + std.error = NA_real_, conf.low = dt[, "lwr.ci"], conf.high = dt[, "upr.ci"], method = "Dunnett_DescTools" ) }, silent = TRUE) } else if (method == "Dunnett_multcomp") { - ##browser() result <- try({ - # For multcomp, we need to create an aov model first if (inherits(x, "formula")) { model_data <- list(...)$data - if (is.null(model_data)) { - model_data <- environment(x) - } + if (is.null(model_data)) { model_data <- environment(x) } aov_model <- aov(x, data = model_data) } else if (inherits(x, c("aov", "lm"))) { aov_model <- x @@ -253,44 +260,31 @@ broom_dunnett <- function(x, method = c("Dunnett_multcomp", "Dunnett_DescTools", stop("x must be a formula, aov, or lm object for Dunnett_multcomp method") } - # Extract the factor name from the formula - if (inherits(x, "formula")) { - terms <- terms(x) - factor_name <- attr(terms, "term.labels")[1] - } else { - factor_name <- names(aov_model$xlevels)[1] - } - - # Get control level - control <- list(...)$control - if (!is.null(control)) { - ## control <- 1 # Default to first level - stop("Please make sure control is the first level of the dose group") - } - - # Create contrast matrix for Dunnett - contrasts_call <- call("mcp") - ##contrasts_call[[2]] <- as.name(factor_name) - contrasts_call[[2]] <- "Dunnett" - names(contrasts_call)[2] <- factor_name - ##contrasts_call$base <- control + factor_name <- names(aov_model$xlevels)[1] + # --- THIS IS THE CORRECTED PART --- + # Create a named list where the name is the factor and the value is "Dunnett" + contrast_arg <- list("Dunnett") + names(contrast_arg) <- factor_name - dunnett_contrasts <- eval(contrasts_call) + # Now create the multiple comparison object + dunnett_mcp <- do.call(multcomp::mcp, contrast_arg) + # This correctly creates the equivalent of `mcp(Dose = "Dunnett")` + # ------------------------------------ - # Run the test - glht_result <- multcomp::glht(aov_model, linfct = dunnett_contrasts) + glht_result <- multcomp::glht(aov_model, linfct = dunnett_mcp, alternative = alternative) summary_result <- summary(glht_result) conf_int <- confint(glht_result) - # Create standardized output tibble::tibble( comparison = names(summary_result$test$coefficients), estimate = summary_result$test$coefficients, + statistic = summary_result$test$tstat, p.value = summary_result$test$pvalues, std.error = summary_result$test$sigma, conf.low = conf_int$confint[, "lwr"], conf.high = conf_int$confint[, "upr"], + significant = summary_result$test$pvalues < alpha, # Basic significance flag method = "Dunnett_multcomp" ) }, silent = TRUE) @@ -298,15 +292,7 @@ broom_dunnett <- function(x, method = c("Dunnett_multcomp", "Dunnett_DescTools", if (inherits(result, "try-error")) { warning("Dunnett test failed. Error: ", attr(result, "condition")$message) - return(tibble::tibble( - comparison = character(), - estimate = numeric(), - p.value = numeric(), - std.error = numeric(), - conf.low = numeric(), - conf.high = numeric(), - method = character() - )) + return(tibble::tibble()) } return(result) diff --git a/R/data_description.R b/R/data_description.R index b12b6f5..990a9ab 100644 --- a/R/data_description.R +++ b/R/data_description.R @@ -101,8 +101,8 @@ NULL #' Dixon's outlier test critical Q table #' #' @author Zhenglei Gao -#' @references DIXON, W. J. (1950) Analysis of extreme values. Ann. Math. Stat. 21, 488–506. -#' DEAN, R. B., DIXON, W. J. (1951) Simplified statistics for small numbers of observation. Anal. Chem. 23, 636–638. +#' @references DIXON, W. J. (1950) Analysis of extreme values. Ann. Math. Stat. 21, 488-506. +#' DEAN, R. B., DIXON, W. J. (1951) Simplified statistics for small numbers of observation. Anal. Chem. 23, 636-638. "DixonQ" diff --git a/R/quantal.R b/R/quantal.R index 8b613bb..30a865e 100644 --- a/R/quantal.R +++ b/R/quantal.R @@ -285,7 +285,8 @@ compare_to_control_fisher <- function(data, factor_col, success_col, failure_col #' #' @author \href{https://stats.stackexchange.com/users/173082/reinstate-monica}{Ben O'Neill} #' @references \url{https://stats.stackexchange.com/a/410376/6378} and -#' R. E. TARONE, Testing the goodness of fit of the binomial distribution, Biometrika, Volume 66, Issue 3, December 1979, Pages 585–590, \url{https://doi.org/10.1093/biomet/66.3.585} +#' R. E. TARONE, Testing the goodness of fit of the binomial distribution, Biometrika, V +#' olume 66, Issue 3, December 1979, Pages 585-590, \url{https://doi.org/10.1093/biomet/66.3.585} #' @importFrom stats pnorm #' @export #' @examples diff --git a/inst/DevNotes/function_dev/dev_dunnett.R b/inst/DevNotes/function_dev/dev_dunnett.R deleted file mode 100644 index 6fa5f85..0000000 --- a/inst/DevNotes/function_dev/dev_dunnett.R +++ /dev/null @@ -1,519 +0,0 @@ -#' Conduct Dunnett Test with Various Model Specifications -#' -#' This function performs Dunnett's test for comparing multiple treatment levels to a control -#' using various model specifications, including options for random effects and variance structures. -#' -#' @param data A data frame containing the dose-response data -#' @param response_var Name of the response variable column -#' @param dose_var Name of the dose/treatment variable column -#' @param block_var Name of the blocking/tank variable column (optional) -#' @param control_level The level of dose_var to use as control (default is minimum dose) -#' @param include_random_effect Logical, whether to include random effects for blocks/tanks -#' @param variance_structure Character, specifying the variance structure: -#' "homoscedastic" (default) or "heteroscedastic" -#' @param alpha Significance level for determining NOEC (default = 0.05) -#' @param conf_level Confidence level for intervals (default = 0.95) -#' @param return_model Logical, whether to return the fitted model object (default = FALSE) -#' -#' @return A list containing the Dunnett test results, NOEC value, and optionally the model object -#' @export -#' -#' @importFrom multcomp glht mcp -#' @importFrom lme4 lmer -#' @importFrom nlme gls lme varIdent -#' @importFrom stats as.formula -dunnett_test <- function(data, - response_var = "Response", - dose_var = "Dose", - block_var = "Tank", - control_level = NULL, - include_random_effect = TRUE, - variance_structure = c("homoscedastic", "heteroscedastic"), - alpha = 0.05, - conf_level = 0.95, - return_model = FALSE) { - - # Input validation - if (!is.data.frame(data)) { - stop("Data must be a data frame") - } - - if (!response_var %in% names(data)) { - stop(paste("Response variable", response_var, "not found in data")) - } - - if (!dose_var %in% names(data)) { - stop(paste("Dose/treatment variable", dose_var, "not found in data")) - } - - # Ensure dose variable is a factor - if (!is.factor(data[[dose_var]])) { - data[[dose_var]] <- factor(data[[dose_var]]) - } - - # Set control level if not specified - if (is.null(control_level)) { - # Use the minimum dose level as control - control_level <- levels(data[[dose_var]])[1] - } else { - # Ensure control_level is in the levels - if (!as.character(control_level) %in% levels(data[[dose_var]])) { - stop("Control level not found in dose variable levels") - } - } - - # Match variance structure argument - variance_structure <- match.arg(variance_structure) - - # Check if block variable exists when random effects are requested - if (include_random_effect && !block_var %in% names(data)) { - stop(paste("Block/tank variable", block_var, "not found in data")) - } - - # Create formula strings - fixed_formula_str <- paste(response_var, "~", dose_var) - fixed_formula <- as.formula(fixed_formula_str) - - # Fit the appropriate model based on specifications - if (include_random_effect) { - if (variance_structure == "homoscedastic") { - # Mixed model with homoscedastic errors - message("Fitting mixed model with homoscedastic errors") - model <- lme4::lmer( - as.formula(paste(fixed_formula_str, "+ (1|", block_var, ")")), - data = data - ) - } else { - # Mixed model with heteroscedastic errors by dose level - message("Fitting mixed model with heteroscedastic errors") - model <- nlme::lme( - fixed = fixed_formula, - random = as.formula(paste("~ 1 |", block_var)), - weights = nlme::varIdent(form = as.formula(paste("~ 1 |", dose_var))), - data = data - ) - } - } else { - if (variance_structure == "homoscedastic") { - # Linear model with homoscedastic errors - message("Fitting linear model with homoscedastic errors") - model <- stats::lm(fixed_formula, data = data) - } else { - # GLS model with heteroscedastic errors by dose level - message("Fitting GLS model with heteroscedastic errors") - ## browser() checked correct - model <- nlme::gls( - fixed_formula, - weights = nlme::varIdent(form = as.formula(paste("~ 1 |", dose_var))), - data = data - ) - } - } - - # Create a list for the mcp call - mcp_list <- list("Dunnett") - names(mcp_list) <- dose_var - - # Set control level if not the first level - if (control_level != levels(data[[dose_var]])[1]) { - control_idx <- which(levels(data[[dose_var]]) == as.character(control_level)) - mcp_list$base <- control_idx - } - - # Perform Dunnett test - ##browser() - dunnett_result <- multcomp::glht(model, linfct = do.call(multcomp::mcp, mcp_list)) - dunnett_summary <- summary(dunnett_result, test = multcomp::adjusted("single-step")) - dunnett_confint <- confint(dunnett_result, level = conf_level) - - # Extract p-values and format comparison results - p_values <- dunnett_summary$test$pvalues - comparisons <- rownames(as.data.frame(dunnett_result$linfct)) - - # Create a data frame with results - results_df <- data.frame( - comparison = comparisons, - estimate = dunnett_summary$test$coefficients, - std.error = dunnett_summary$test$sigma, - statistic = dunnett_summary$test$tstat, - p.value = p_values, - conf.low = dunnett_confint$confint[, "lwr"], - conf.high = dunnett_confint$confint[, "upr"], - significant = p_values < alpha - ) - - # Determine NOEC (No Observed Effect Concentration) - # Extract dose levels from comparison strings and convert to numeric - dose_levels <- sapply(strsplit(comparisons, " - "), function(x) x[1]) - - # Convert to numeric if possible - numeric_doses <- suppressWarnings(as.numeric(dose_levels)) - if (all(!is.na(numeric_doses))) { - dose_levels <- numeric_doses - } - - # Find the highest dose with non-significant effect - significant_effects <- p_values < alpha - if (all(significant_effects)) { - noec <- min(dose_levels) # All doses show effects, NOEC is below lowest dose - noec_message <- "All tested doses show significant effects. NOEC is below the lowest tested dose." - } else if (!any(significant_effects)) { - noec <- max(dose_levels) # No doses show effects, NOEC is at or above highest dose - noec_message <- "No significant effects detected at any dose. NOEC is at or above the highest tested dose." - } else { - # Find the highest non-significant dose - non_sig_doses <- dose_levels[!significant_effects] - sig_doses <- dose_levels[significant_effects] - - # Ensure we're working with proper numeric values for comparison - if (is.numeric(non_sig_doses) && is.numeric(sig_doses)) { - noec <- max(non_sig_doses[non_sig_doses < max(sig_doses)]) - } else { - # If doses aren't numeric, just return the highest non-significant level - noec <- non_sig_doses[length(non_sig_doses)] - } - noec_message <- paste("NOEC determined as", noec) - } - - # Prepare return object - result <- list( - dunnett_test = dunnett_summary, - results_table = results_df, - noec = noec, - noec_message = noec_message, - model_type = paste0( - ifelse(include_random_effect, "Mixed", "Fixed"), - " model with ", - variance_structure, - " errors" - ), - control_level = control_level, - alpha = alpha - ) - - if (return_model) { - result$model <- model - } - - class(result) <- "dunnett_test_result" - - return(result) -} - - - - - - - - - - - - - - - - - - - - -#' Conduct Dunnett Test with Various Model Specifications -#' -#' This function performs Dunnett's test for comparing multiple treatment levels to a control -#' using various model specifications, including options for random effects and variance structures. -#' -#' @param data A data frame containing the dose-response data -#' Conduct Dunnett Test with Various Model Specifications -#' -#' This function performs Dunnett's test for comparing multiple treatment levels to a control -#' using various model specifications, including options for random effects and variance structures. -#' -#' @param data A data frame containing the dose-response data -#' @param response_var Name of the response variable column -#' @param dose_var Name of the dose/treatment variable column -#' @param block_var Name of the blocking/tank variable column (optional) -#' @param control_level The level of dose_var to use as control (default is minimum dose) -#' @param include_random_effect Logical, whether to include random effects for blocks/tanks -#' @param variance_structure Character, specifying the variance structure: -#' "homoscedastic" (default) or "heteroscedastic" -#' @param alpha Significance level for determining NOEC (default = 0.05) -#' @param conf_level Confidence level for intervals (default = 0.95) -#' @param return_model Logical, whether to return the fitted model object (default = FALSE) -#' -#' @return A list containing the Dunnett test results, NOEC value, and optionally the model object -#' @export -#' -#' @importFrom multcomp glht mcp -#' @importFrom lme4 lmer -#' @importFrom nlme gls lme varIdent -#' @importFrom stats as.formula -dunnett_test <- function(data, - response_var = "Response", - dose_var = "Dose", - block_var = "Tank", - control_level = NULL, - include_random_effect = TRUE, - variance_structure = c("homoscedastic", "heteroscedastic"), - alpha = 0.05, - conf_level = 0.95, - return_model = FALSE) { - - # Input validation - if (!is.data.frame(data)) { - stop("Data must be a data frame") - } - - if (!response_var %in% names(data)) { - stop(paste("Response variable", response_var, "not found in data")) - } - - if (!dose_var %in% names(data)) { - stop(paste("Dose/treatment variable", dose_var, "not found in data")) - } - - # Ensure dose variable is a factor - if (!is.factor(data[[dose_var]])) { - data[[dose_var]] <- factor(data[[dose_var]]) - } - - # Set control level if not specified - if (is.null(control_level)) { - # Use the minimum dose level as control - control_level <- levels(data[[dose_var]])[1] - } else { - # Ensure control_level is in the levels - if (!as.character(control_level) %in% levels(data[[dose_var]])) { - stop("Control level not found in dose variable levels") - } - } - - # Match variance structure argument - variance_structure <- match.arg(variance_structure) - - # Check if block variable exists when random effects are requested - if (include_random_effect && !block_var %in% names(data)) { - stop(paste("Block/tank variable", block_var, "not found in data")) - } - - # Create formula strings - fixed_formula_str <- paste(response_var, "~", dose_var) - fixed_formula <- as.formula(fixed_formula_str) - - # Fit the appropriate model based on specifications - if (include_random_effect) { - if (variance_structure == "homoscedastic") { - # Mixed model with homoscedastic errors - message("Fitting mixed model with homoscedastic errors") - model <- lme4::lmer( - as.formula(paste(fixed_formula_str, "+ (1|", block_var, ")")), - data = data - ) - } else { - # Mixed model with heteroscedastic errors by dose level - message("Fitting mixed model with heteroscedastic errors") - model <- nlme::lme( - fixed = fixed_formula, - random = as.formula(paste("~ 1 |", block_var)), - weights = nlme::varIdent(form = as.formula(paste("~ 1 |", dose_var))), - data = data - ) - } - } else { - if (variance_structure == "homoscedastic") { - # Linear model with homoscedastic errors - message("Fitting linear model with homoscedastic errors") - model <- stats::lm(fixed_formula, data = data) - } else { - # GLS model with heteroscedastic errors by dose level - message("Fitting GLS model with heteroscedastic errors") - model <- nlme::gls( - fixed_formula, - weights = nlme::varIdent(form = as.formula(paste("~ 1 |", dose_var))), - data = data - ) - } - } - - # Create contrast for Dunnett test - # This is the corrected part that properly handles variable names - linfct <- NULL - - if (inherits(model, "lmerMod") || inherits(model, "lm") || inherits(model, "lme")) { - # For lmer,lme and lm models - dunnett_args <- list(model) - mc_formula <- paste(dose_var, "= 'Dunnett'") - mc_call <- call("mcp") - mc_call[[dose_var]] <- "Dunnett" - - # Set control level if not the first level - if (control_level != levels(data[[dose_var]])[1]) { - mc_call$base <- which(levels(data[[dose_var]]) == as.character(control_level)) - } - - dunnett_args$linfct <- mc_call - dunnett_result <- do.call(multcomp::glht, dunnett_args) - - } else if (inherits(model, "gls")) { - # For nlme models (lme, gls) - # Create a contrast matrix manually - browser() - n_levels <- nlevels(data[[dose_var]]) - control_idx <- which(levels(data[[dose_var]]) == as.character(control_level)) - - # Create Dunnett contrast matrix - K <- matrix(0, n_levels - 1, n_levels) - row_idx <- 1 - for (i in 1:n_levels) { - if (i != control_idx) { - K[row_idx, i] <- 1 # Treatment level - #K[row_idx, control_idx] <- -1 # Control level - row_idx <- row_idx + 1 - } - } - - # Create row names for the contrast matrix - level_names <- levels(data[[dose_var]]) - row_names <- character(n_levels - 1) - row_idx <- 1 - for (i in 1:n_levels) { - if (i != control_idx) { - row_names[row_idx] <- paste(level_names[i], "-", level_names[control_idx]) - row_idx <- row_idx + 1 - } - } - rownames(K) <- row_names - - # Create the contrast - linfct <- multcomp::glht(model, linfct = K) - dunnett_result <- linfct - } - - # Get test results - dunnett_summary <- summary(dunnett_result, test = multcomp::adjusted("single-step")) - dunnett_confint <- confint(dunnett_result, level = conf_level) - - # Extract p-values and format comparison results - p_values <- dunnett_summary$test$pvalues - ##browser() - comparisons <- rownames(as.data.frame(dunnett_result$linfct)) - - # Create a data frame with results - results_df <- data.frame( - comparison = comparisons, - estimate = dunnett_summary$test$coefficients, - std.error = dunnett_summary$test$sigma, - statistic = dunnett_summary$test$tstat, - p.value = p_values, - conf.low = dunnett_confint$confint[, "lwr"], - conf.high = dunnett_confint$confint[, "upr"], - significant = p_values < alpha - ) - - # Determine NOEC (No Observed Effect Concentration) - # Extract dose levels from comparison strings and convert to numeric - dose_levels <- sapply(strsplit(comparisons, " - "), function(x) x[1]) - - # Convert to numeric if possible - numeric_doses <- suppressWarnings(as.numeric(dose_levels)) - if (all(!is.na(numeric_doses))) { - dose_levels <- numeric_doses - } - - # Find the highest dose with non-significant effect - significant_effects <- p_values < alpha - if (all(significant_effects)) { - noec <- min(dose_levels) # All doses show effects, NOEC is below lowest dose - noec_message <- "All tested doses show significant effects. NOEC is below the lowest tested dose." - } else if (!any(significant_effects)) { - noec <- max(dose_levels) # No doses show effects, NOEC is at or above highest dose - noec_message <- "No significant effects detected at any dose. NOEC is at or above the highest tested dose." - } else { - # Find the highest non-significant dose - non_sig_doses <- dose_levels[!significant_effects] - sig_doses <- dose_levels[significant_effects] - - # Ensure we're working with proper numeric values for comparison - if (is.numeric(non_sig_doses) && is.numeric(sig_doses)) { - noec <- max(non_sig_doses[non_sig_doses < max(sig_doses)]) - } else { - # If doses aren't numeric, just return the highest non-significant level - noec <- non_sig_doses[length(non_sig_doses)] - } - noec_message <- paste("NOEC determined as", noec) - } - - # Prepare return object - result <- list( - dunnett_test = dunnett_summary, - results_table = results_df, - noec = noec, - noec_message = noec_message, - model_type = paste0( - ifelse(include_random_effect, "Mixed", "Fixed"), - " model with ", - variance_structure, - " errors" - ), - control_level = control_level, - alpha = alpha - ) - - if (return_model) { - result$model <- model - } - - class(result) <- "dunnett_test_result" - - return(result) -} - -#' Print method for dunnett_test_result objects -#' -#' @param x A dunnett_test_result object -#' @param ... Additional arguments passed to print methods -#' -#' @export -print.dunnett_test_result <- function(x, ...) { - cat("Dunnett Test Results\n") - cat("-------------------\n") - cat("Model type:", x$model_type, "\n") - cat("Control level:", x$control_level, "\n") - cat("Alpha level:", x$alpha, "\n\n") - - cat("Results Table:\n") - print(x$results_table, row.names = FALSE) - - cat("\nNOEC Determination:\n") - cat(x$noec_message, "\n") -} - -#' Plot method for dunnett_test_result objects -#' -#' @param x A dunnett_test_result object -#' @param ... Additional arguments passed to plot methods -#' -#' @importFrom ggplot2 ggplot aes geom_point geom_errorbar theme_minimal labs geom_hline -#' @export -plot.dunnett_test_result <- function(x, ...) { - # Extract data for plotting - plot_data <- x$results_table - plot_data$comparison <- factor(plot_data$comparison,levels=plot_data$comparison) - # Create the plot - p <- ggplot2::ggplot(plot_data, ggplot2::aes(x = comparison, y = estimate, color = significant)) + - ggplot2::geom_point(size = 3) + - ggplot2::geom_errorbar(ggplot2::aes(ymin = conf.low, ymax = conf.high), width = 0.2) + - ggplot2::geom_hline(yintercept = 0, linetype = "dashed") + - ggplot2::theme_minimal() + - ggplot2::labs( - title = paste("Dunnett Test Results:", x$model_type), - subtitle = paste("NOEC =", x$noec), - x = "Comparison", - y = "Difference from Control", - color = "Significant" - ) + - ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) - - return(p) -} diff --git a/man/DixonQ.Rd b/man/DixonQ.Rd index 5d74f8a..654f913 100644 --- a/man/DixonQ.Rd +++ b/man/DixonQ.Rd @@ -14,8 +14,8 @@ DixonQ Dixon's outlier test critical Q table } \references{ -DIXON, W. J. (1950) Analysis of extreme values. Ann. Math. Stat. 21, 488–506. -DEAN, R. B., DIXON, W. J. (1951) Simplified statistics for small numbers of observation. Anal. Chem. 23, 636–638. +DIXON, W. J. (1950) Analysis of extreme values. Ann. Math. Stat. 21, 488-506. +DEAN, R. B., DIXON, W. J. (1951) Simplified statistics for small numbers of observation. Anal. Chem. 23, 636-638. } \author{ Zhenglei Gao diff --git a/man/Tarone.test.Rd b/man/Tarone.test.Rd index bd93c22..04a6684 100644 --- a/man/Tarone.test.Rd +++ b/man/Tarone.test.Rd @@ -25,7 +25,8 @@ Tarone.test(N, M) } \references{ \url{https://stats.stackexchange.com/a/410376/6378} and -R. E. TARONE, Testing the goodness of fit of the binomial distribution, Biometrika, Volume 66, Issue 3, December 1979, Pages 585–590, \url{https://doi.org/10.1093/biomet/66.3.585} +R. E. TARONE, Testing the goodness of fit of the binomial distribution, Biometrika, V +olume 66, Issue 3, December 1979, Pages 585-590, \url{https://doi.org/10.1093/biomet/66.3.585} } \author{ \href{https://stats.stackexchange.com/users/173082/reinstate-monica}{Ben O'Neill} diff --git a/man/broom_dunnett.Rd b/man/broom_dunnett.Rd index 516cb19..b691a2a 100644 --- a/man/broom_dunnett.Rd +++ b/man/broom_dunnett.Rd @@ -7,6 +7,8 @@ broom_dunnett( x, method = c("Dunnett_multcomp", "Dunnett_DescTools", "Dunnett_PMCMRplus"), + alternative = c("two.sided", "less", "greater"), + alpha = 0.05, ... ) } diff --git a/man/compute_mdd_dunnett.Rd b/man/compute_mdd_dunnett.Rd new file mode 100644 index 0000000..c5bb530 --- /dev/null +++ b/man/compute_mdd_dunnett.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MDD.R +\name{compute_mdd_dunnett} +\alias{compute_mdd_dunnett} +\title{Calculate MDD\% for a Dunnett Test Result} +\usage{ +compute_mdd_dunnett(dunnett_obj, alternative, data, formula) +} +\arguments{ +\item{dunnett_obj}{The result object from your dunnett_test function.} + +\item{alternative}{The alternative hypothesis ("less", "greater", "two.sided").} + +\item{data}{The original dataframe used for the test.} + +\item{formula}{The formula used for the test.} +} +\value{ +A tibble with Dose, MDD, and MDD_pct. +Calculate MDD\% for a Dunnett Test Result (Robust Version 2) + +A tibble with Dose and MDD_pct. +} +\description{ +Calculate MDD\% for a Dunnett Test Result +} diff --git a/man/report_dunnett_summary.Rd b/man/report_dunnett_summary.Rd new file mode 100644 index 0000000..5f74c4b --- /dev/null +++ b/man/report_dunnett_summary.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MDD.R +\name{report_dunnett_summary} +\alias{report_dunnett_summary} +\title{Generate a Comprehensive Dunnett Test Summary Report (Robust Version)} +\usage{ +report_dunnett_summary( + formula, + data, + dunnett_test_func, + alternative = "less", + include_random_effect = FALSE, + ... +) +} +\arguments{ +\item{formula}{The model formula, e.g., Response ~ Dose.} + +\item{data}{The dataframe containing the data.} + +\item{dunnett_test_func}{The dunnett_test function you use.} + +\item{alternative}{The alternative hypothesis to be used for the test.} + +\item{...}{Additional arguments passed to your dunnett_test_func.} +} +\value{ +A single tibble combining descriptive stats, test results, and MDD\%. +} +\description{ +Generate a Comprehensive Dunnett Test Summary Report (Robust Version) +} diff --git a/vignettes/Dunnetts_Test_for_Data_with_Hierarchical_Structure.Rmd b/vignettes/Dunnetts_Test_for_Data_with_Hierarchical_Structure.Rmd index 390ff82..d59c388 100644 --- a/vignettes/Dunnetts_Test_for_Data_with_Hierarchical_Structure.Rmd +++ b/vignettes/Dunnetts_Test_for_Data_with_Hierarchical_Structure.Rmd @@ -216,3 +216,21 @@ result4$results_table$estimate -s4$test$coefficients result4$results_table$p.value - s4$test$pvalues < 1e-05 ``` +## Dunnett's Test with Different R Packages + +We recommend to use the `multcomp` package to perform the routine Dunnett's test. + +```{r} +# Define Methods and Run All +methods_to_test <- c("Dunnett_PMCMRplus", "Dunnett_DescTools", "Dunnett_multcomp") + +all_results <- purrr::map_dfr(methods_to_test, ~ broom_dunnett( + x = Response ~ factor(Dose), + data = dat_medium, + method = .x, + alternative = "two.sided" +)) + +all_results %>% knitr::kable(.,digits = 3) +``` + diff --git a/vignettes/articles/MDD-in-Regulatory-Context.Rmd b/vignettes/articles/MDD-in-Regulatory-Context.Rmd index ce9b4ae..a7896fa 100644 --- a/vignettes/articles/MDD-in-Regulatory-Context.Rmd +++ b/vignettes/articles/MDD-in-Regulatory-Context.Rmd @@ -17,7 +17,7 @@ source("../knitr-setup.R") ``` -In the regulatory context, MDD% stands for the minimum percent change (relative to the control mean) that would be deemed statistically significant given the test's critical value, the residual variability, and the replication. There are two common ways to define it; the often used way in ecotoxicology is ofent the significance-threshold version (no power term). +In the regulatory context, MDD% stands for the minimum percent change (relative to the control mean) that would be deemed statistically significant given the test's critical value, the residual variability, and the replication. This **IS NOT** a power based MDD as we usually would have defined in statistical power analysis. There are two common ways to define it; the often used way in ecotoxicology is ofent the significance-threshold version (no power term). ## Key definitions @@ -346,61 +346,70 @@ Reusing existing test objects and building a modular system around them is far m ```{r} library(drcHelper) Results <- broom_williams(Response ~ factor(Dose), data = dat_medium, - test = "Williams_JG", alpha = 0.05, alternative = "less") + method = "Williams_PMCMRplus", alpha = 0.05, alternative = "less") compute_mdd_williams(Results,data=dat_medium,formula = Response ~ factor(Dose)) +Results <- dunnett_test(data = dat_medium,response_var = "Response",tank_var = NULL,include_random_effect = FALSE,dose_var = "Dose", alpha = 0.05, alternative = "less") ``` +```{r} +Results <- broom_dunnett(Response ~ factor(Dose), data = dat_medium, + method = "Dunnett_multcomp", alpha = 0.05, alternative = "two.sided") +Results +Results <- dunnett_test(data = dat_medium,response_var = "Response",tank_var = NULL,include_random_effect = FALSE,dose_var = "Dose", alpha = 0.05, alternative = "two.sided") +Results -# Implementation in the Validation for Williams -Here's a self-contained helper to compute per-dose MDD% for Williams and return a tidy table. It uses your existing helpers and broom_williams output. +Results <- broom_dunnett(Response ~ factor(Dose), data = dat_medium, + method = "Dunnett_multcomp", alpha = 0.05, alternative = "less") +Results +Results <- dunnett_test(data = dat_medium,response_var = "Response",tank_var = NULL,include_random_effect = FALSE,dose_var = "Dose", alpha = 0.05, alternative = "less") +compute_mdd_dunnett(Results,data = dat_medium,formula = Response~factor(Dose),alternative = "less") -# How to integrate into the validation workflow +Results <- dunnett_test(data = dat_medium,response_var = "Response",tank_var = NULL,include_random_effect = FALSE,dose_var = "Dose", alpha = 0.05, alternative = "two.sided") +compute_mdd_dunnett(Results,data = dat_medium,formula = Response~factor(Dose),alternative = "two.sided") +Results -- In the Williams run_actual handler, compute the prelim summary as before (for Mean, SD, %Inhibition, CV), and also compute MDD%: +``` -```r -prelim <- compute_prelim_summary(endpoint_data) -mdd_tbl <- compute_mdd_williams(endpoint_data, alternative = alternative) +### Overall reporting -# Return both; the generic runner can join MDD% when expected rows exist -list(actual_df = actual_df, group_means = prelim, mdd = mdd_tbl) -``` +- both functions `dunnett_test` and `broom_dunnet` can be used now. +- note `dunnett_test` has more options to perform test including random effect and inhomogeneous variance, while `broom_dunnet` is the simple version that is similar to Williams' test. -- In your expected parsing (parse_expected_metrics for Williams), you already capture "MDD%": +```{r} +# This call will now work correctly for any alternative. +final_dunnett_table <- report_dunnett_summary( + formula = Response ~ Dose, + data = dat_medium, + dunnett_test_func = "dunnett_test", + alternative = "less", + alpha = 0.05 +) -```r -out$mdd <- pick("\\bMDD%\\b"); -names(out$mdd)[2] <- "Expected_MDDpct" -``` +print(final_dunnett_table) -- In the generic runner's join stage for Williams, add: +# This call will now work correctly for any alternative. +final_dunnett_table <- report_dunnett_summary( + formula = Response ~ Dose, + data = dat_medium, + dunnett_test_func = "broom_dunnett", + alternative = "less", + alpha = 0.05 +) -```{r eval=FALSE} -if (!is.null(exp_splits$mdd) && nrow(exp_splits$mdd)) { - mdd_join <- exp_splits$mdd %>% - dplyr::left_join(actual$mdd[, c("Dose", "Actual_MDDpct")], by = "Dose") %>% - dplyr::mutate( - Endpoint = endpoint, - Diff = abs(Actual_MDDpct - Expected_MDDpct), - Status = dplyr::case_when( - is.na(Expected_MDDpct) | is.na(Actual_MDDpct) ~ "MISSING", - Diff <= tolerance ~ "PASS", TRUE ~ "FAIL") - ) %>% - dplyr::transmute(Endpoint, Dose, metric = "MDD%", Actual = Actual_MDDpct, Expected = Expected_MDDpct, Status) - pieces <- c(pieces, list(mdd_join)) -} +print(final_dunnett_table) ``` -### Handling "greater" vs "smaller" -- The formula for MDD% doesn't change with direction; it's the minimum absolute difference required for significance. You typically compare magnitudes (positive values). If you want to encode a "directional" MDD% (signed), you can attach sign = ifelse(alternative=="greater", +1, -1) and report MDD%*sign, but most reports use a positive threshold. +### Handling "greater" vs "smaller" in Williams' test + +- The formula for MDD% doesn't change with direction for Williams' test; it's the minimum absolute difference required for significance. You typically compare magnitudes (positive values). If you want to encode a "directional" MDD% (signed), you can attach sign = ifelse(alternative=="greater", +1, -1) and report MDD%*sign, but most reports use a positive threshold. ### Unequal sample sizes and heteroscedasticity diff --git a/vignettes/articles/Verification-of-Williams-Test.Rmd b/vignettes/articles/Verification-of-Williams-Test.Rmd index aa4bff1..eab03fb 100644 --- a/vignettes/articles/Verification-of-Williams-Test.Rmd +++ b/vignettes/articles/Verification-of-Williams-Test.Rmd @@ -15,17 +15,17 @@ library(drcHelper) -Short answer: Yes for the PMCMRplus method. broom_williams doesn’t have a formal alternative argument, but it forwards any extra arguments via ... to PMCMRplus::williamsTest, so you can pass alternative = "less" or alternative = "greater". For the Williams_JG method, you should use direction = "decreasing" or direction = "increasing" instead of alternative. +## Which Package to Use -How to verify with your dat_medium +Short answer: Use `drcHelper::broom_williams` with the PMCMRplus method. It handles formula gracefully. `broom_williams` doesn't have a formal alternative argument, but it forwards any extra arguments via ... to PMCMRplus::williamsTest, so you can pass alternative = "less" or alternative = "greater". For the Williams_JG method, you should use direction = "decreasing" or direction = "increasing" instead of alternative. + +### How to use the function with dat_medium - Ensure the control is the first level of the dose factor (lowest dose first). - Call broom_williams with method = "Williams_PMCMRplus" and pass alternative via ... - Observe the "comparison" column: for alternative = "greater" you’ll see "<= 0"; for alternative = "less" you’ll see ">= 0". - - -# Prepare the data: control first level +### Prepare the data: control first level ```{r} dm <- dat_medium dm$Dose_factor <- factor(dm$Dose, levels = sort(unique(dm$Dose))) @@ -67,6 +67,7 @@ bw_inc <- drcHelper::broom_williams(Response ~ Dose_factor, print(bw_inc) ``` +**Note that the t'-stat output from StatCharms can be of incorrect signs at it returns the absolute value of the Williams' t-statistic, which should match the sign of the estimated differences** ## What to look for diff --git a/vignettes/knitr-setup.R b/vignettes/knitr-setup.R index 66ed9f2..ff4b8c8 100644 --- a/vignettes/knitr-setup.R +++ b/vignettes/knitr-setup.R @@ -7,9 +7,11 @@ knitr::opts_chunk$set( # e.g., for 'My-Vignette.Rmd', figures will be named 'My-Vignette-chunk-label-1.png' # and placed in the 'articles' directory during the pkgdown build. ## fig.path = paste0(tools::file_path_sans_ext(basename(knitr::current_input())), "-"), - + # Provide a default, non-descriptive alt text. # This is better than no alt text, but it's highly recommended to # provide specific, descriptive alt text for each plot using fig.alt in the chunk options. fig.alt = "A plot generated from an R code chunk." ) + +options(knitr.kable.NA = "") From 808614f91bf5bbe1127738cfd9f3ae3c55fc7577 Mon Sep 17 00:00:00 2001 From: Zhenglei Gao Date: Thu, 2 Oct 2025 20:42:08 +0200 Subject: [PATCH 2/5] fledge: Bump version to 0.0.4.9001 --- DESCRIPTION | 2 +- NEWS.md | 97 +++++++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 96 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5a9da13..0931312 100644 --- a/DESCRIPTION +++ b/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.4.9001 Authors@R: c( person("Zhenglei", "Gao", email="zhenglei.gao@bayer.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-4042-310X")), diff --git a/NEWS.md b/NEWS.md index 61287bd..6905f01 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,96 @@ -# drcHelper (development version) + + +# drcHelper 0.0.4.9001 + +## Features + +- Track knitr-setup.R and update .gitignore. + +## Uncategorized + +- Merge branch 'dev' of https://github.com/Bayer-Group/drcHelper into dev. + +- Introduced USER_QUESTIONS_ANSWERED.md for detailed answers to user queries regarding validation results and data quality issues. + +- Implemented comprehensive_multi_study_analysis.R to analyze multiple studies with Dunnett tests, summarizing results and endpoints. + +- Created detailed_individual_analysis.R for in-depth validation of individual function groups, including actual vs expected comparisons. + +- Developed detailed_validation_functions.R with enhanced validation functions providing detailed row-by-row results for Dunnett tests. + +- Generated final_validation_summary.R to summarize the overall validation results, including success rates and study breakdowns. + +- Added working_detailed_analysis.R for a practical demonstration of validation processes and handling of data quality issues. + +- Created definitive_multi_study_multi_endpoint_demo.R to showcase the package's capability to handle multiple studies and endpoints effectively. + +- Implemented `run_dunnett_validation` function to support multi-endpoint validation. + +- Created scripts to check study IDs, data columns, and expected results for function groups. + +- Added debugging scripts to investigate issues with expected results and data linking for FG00225. + +- Developed tests for individual function groups and multi-endpoint validation. + +- Enhanced error handling and output reporting in validation functions. + +- Included checks for European decimal notation in dose conversion. + +- Established a framework for validating T-statistics, P-values, and means against expected results. + +- Created a comprehensive validation summary for the Dunnett Test, detailing data quality issues and actions taken. + +- Implemented a test script (`test_dunnett_call.R`) to evaluate the functionality of the `dunnett_test` function with specific study data. + +- Developed a validation script (`test_fixed_patterns.R`) to compare actual results against expected values using corrected patterns. + +- Added a script (`test_updated_tolerances.R`) to test the Dunnett Test with updated tolerances for t-values and p-values. + +- Included data cleaning steps to handle invalid placeholders and ensure proper statistical comparisons. + +- Enhanced error handling in test scripts to capture and report issues during validation. + +- Implemented a new script to analyze which endpoints have count vs continuous data. + +- Added detailed output for each endpoint, including data presence checks for Total, Alive, and Dead columns. + +- Developed a corrected validation function that checks count data per endpoint instead of per study. + +- Enhanced the function to handle both count and continuous data appropriately, including control level determination. + +- Tested the fixed validation logic across multiple function groups with comprehensive output. + +- Implemented a comprehensive test script (test_comprehensive.R) to verify the correct matching logic for data, focusing on the Myriophyllum study (MOCK0065) and other studies. + +- Created a focused test script (test_matching_issue.R) to demonstrate specific matching problems and the correct logic for handling measurement variables. + +- Developed a test script (test_matching_logic.R) to verify the structure of test data and results, ensuring the matching logic is correctly applied. + +- Added a validation test script (test_validation.R) to simulate actual testing functions and confirm the matching logic works as intended. + +- Enhanced documentation and comments throughout the scripts to clarify the matching rules and logic applied. + +- Introduced a comprehensive guide for the Statistical Test Validation Framework in `Statistical_Test_Framework_Guide.Rmd`. + +- Created configuration file `test_framework_config.R` to define available statistical tests, their properties, and tolerance settings. + +- Implemented report generation script `generate_test_reports.R` for creating validation reports for statistical tests. + +- Developed a master template `statistical_test_template.Rmd` for generating detailed validation reports with dynamic content. + +- Established a modular design for easy extension and integration of new statistical tests within the framework. + +- Merge pull request #20 from Bayer-Group/dev. + + Dev + +- Merge pull request #17 from Bayer-Group/dev. + +- Merge pull request #13 from Bayer-Group/copilot/fix-7be30736-01c8-4e5e-9ad6-508c41c5d437. + +- Merge pull request #11 from Bayer-Group/dev. + + added quantal example data so that page deployment can be run + +- Initial CRAN submission preparation. -* Initial CRAN submission preparation. From 7a21afb4a028b356624deb078089d527e23311ca Mon Sep 17 00:00:00 2001 From: Zhenglei Gao Date: Thu, 2 Oct 2025 21:23:42 +0200 Subject: [PATCH 3/5] Added MDD related functions, added compare_to_control_welch, bumped version to 0.0.4.9001 --- NAMESPACE | 4 + NEWS.md | 104 +---- R/MDD.R | 62 ++- R/continuous_tests.R | 151 ++++++ R/{quantal.R => quantal_tests.R} | 0 man/Tarone.test.Rd | 2 +- man/compare_to_control_fisher.Rd | 2 +- man/compare_to_control_welch.Rd | 68 +++ man/create_contingency_table.Rd | 2 +- man/many_to_one_fisher_test.Rd | 2 +- man/report_dunnett_summary.Rd | 61 ++- vignettes/articles/A-Notes-on-Power.Rmd | 142 ++++++ .../articles/MDD-in-Regulatory-Context.Rmd | 441 +++++++----------- vignettes/articles/Quantal-Data.Rmd | 90 ++++ 14 files changed, 743 insertions(+), 388 deletions(-) create mode 100644 R/continuous_tests.R rename R/{quantal.R => quantal_tests.R} (100%) create mode 100644 man/compare_to_control_welch.Rd create mode 100644 vignettes/articles/A-Notes-on-Power.Rmd diff --git a/NAMESPACE b/NAMESPACE index 9ce011d..0c3c4ee 100644 --- a/NAMESPACE +++ b/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) @@ -74,6 +76,7 @@ export(prelimPlot3) export(prelimSummary) export(prepDataRSCABS) export(rankTransform) +export(report_dunnett_summary) export(reshape_drcData) export(runRSCABS) export(run_RSCA) @@ -177,6 +180,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) diff --git a/NEWS.md b/NEWS.md index 6905f01..09b441f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,96 +1,18 @@ - +# drcHelper 0.0.4.9001() -# drcHelper 0.0.4.9001 +* Initial CRAN submission preparation. -## Features +## 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 +- 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` -- Track knitr-setup.R and update .gitignore. +## Improvments +- Added test cases validation reports. + - `/inst/SystemTesting/Consolidated_Dunnett_Report.Rmd` -## Uncategorized - -- Merge branch 'dev' of https://github.com/Bayer-Group/drcHelper into dev. - -- Introduced USER_QUESTIONS_ANSWERED.md for detailed answers to user queries regarding validation results and data quality issues. - -- Implemented comprehensive_multi_study_analysis.R to analyze multiple studies with Dunnett tests, summarizing results and endpoints. - -- Created detailed_individual_analysis.R for in-depth validation of individual function groups, including actual vs expected comparisons. - -- Developed detailed_validation_functions.R with enhanced validation functions providing detailed row-by-row results for Dunnett tests. - -- Generated final_validation_summary.R to summarize the overall validation results, including success rates and study breakdowns. - -- Added working_detailed_analysis.R for a practical demonstration of validation processes and handling of data quality issues. - -- Created definitive_multi_study_multi_endpoint_demo.R to showcase the package's capability to handle multiple studies and endpoints effectively. - -- Implemented `run_dunnett_validation` function to support multi-endpoint validation. - -- Created scripts to check study IDs, data columns, and expected results for function groups. - -- Added debugging scripts to investigate issues with expected results and data linking for FG00225. - -- Developed tests for individual function groups and multi-endpoint validation. - -- Enhanced error handling and output reporting in validation functions. - -- Included checks for European decimal notation in dose conversion. - -- Established a framework for validating T-statistics, P-values, and means against expected results. - -- Created a comprehensive validation summary for the Dunnett Test, detailing data quality issues and actions taken. - -- Implemented a test script (`test_dunnett_call.R`) to evaluate the functionality of the `dunnett_test` function with specific study data. - -- Developed a validation script (`test_fixed_patterns.R`) to compare actual results against expected values using corrected patterns. - -- Added a script (`test_updated_tolerances.R`) to test the Dunnett Test with updated tolerances for t-values and p-values. - -- Included data cleaning steps to handle invalid placeholders and ensure proper statistical comparisons. - -- Enhanced error handling in test scripts to capture and report issues during validation. - -- Implemented a new script to analyze which endpoints have count vs continuous data. - -- Added detailed output for each endpoint, including data presence checks for Total, Alive, and Dead columns. - -- Developed a corrected validation function that checks count data per endpoint instead of per study. - -- Enhanced the function to handle both count and continuous data appropriately, including control level determination. - -- Tested the fixed validation logic across multiple function groups with comprehensive output. - -- Implemented a comprehensive test script (test_comprehensive.R) to verify the correct matching logic for data, focusing on the Myriophyllum study (MOCK0065) and other studies. - -- Created a focused test script (test_matching_issue.R) to demonstrate specific matching problems and the correct logic for handling measurement variables. - -- Developed a test script (test_matching_logic.R) to verify the structure of test data and results, ensuring the matching logic is correctly applied. - -- Added a validation test script (test_validation.R) to simulate actual testing functions and confirm the matching logic works as intended. - -- Enhanced documentation and comments throughout the scripts to clarify the matching rules and logic applied. - -- Introduced a comprehensive guide for the Statistical Test Validation Framework in `Statistical_Test_Framework_Guide.Rmd`. - -- Created configuration file `test_framework_config.R` to define available statistical tests, their properties, and tolerance settings. - -- Implemented report generation script `generate_test_reports.R` for creating validation reports for statistical tests. - -- Developed a master template `statistical_test_template.Rmd` for generating detailed validation reports with dynamic content. - -- Established a modular design for easy extension and integration of new statistical tests within the framework. - -- Merge pull request #20 from Bayer-Group/dev. - - Dev - -- Merge pull request #17 from Bayer-Group/dev. - -- Merge pull request #13 from Bayer-Group/copilot/fix-7be30736-01c8-4e5e-9ad6-508c41c5d437. - -- Merge pull request #11 from Bayer-Group/dev. - - added quantal example data so that page deployment can be run - -- Initial CRAN submission preparation. +## Bug Fixes +- Updated `dunnett_test` to correctly handle non-RCBD replicate IDs. diff --git a/R/MDD.R b/R/MDD.R index ab31ba0..8fe11c4 100644 --- a/R/MDD.R +++ b/R/MDD.R @@ -72,7 +72,7 @@ compute_mdd_williams <- function(williams_obj, data, formula) { #' @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 @@ -129,14 +129,58 @@ compute_mdd_dunnett <- function(dunnett_obj, alternative, data, formula) { } -#' Generate a Comprehensive Dunnett Test Summary Report (Robust Version) +#' 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 +#' ) #' -#' @param formula The model formula, e.g., Response ~ Dose. -#' @param data The dataframe containing the data. -#' @param dunnett_test_func The dunnett_test function you use. -#' @param alternative The alternative hypothesis to be used for the test. -#' @param ... Additional arguments passed to your dunnett_test_func. -#' @return A single tibble combining descriptive stats, test results, and MDD%. +#' print(summary_results) +#' } report_dunnett_summary <- function(formula, data, dunnett_test_func, alternative = "less", include_random_effect = FALSE,...) { @@ -147,7 +191,7 @@ report_dunnett_summary <- function(formula, data, dunnett_test_func, alternative 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( + if(dunnett_test_func == "dunnett_test") dunnett_results_obj <- dunnett_test( data = data, response_var = resp_name, dose_var = dose_name, diff --git a/R/continuous_tests.R b/R/continuous_tests.R new file mode 100644 index 0000000..959ee58 --- /dev/null +++ b/R/continuous_tests.R @@ -0,0 +1,151 @@ +#' Perform Many-to-One Welch's t-Tests +#' +#' Compares each level of a factor to a specified control level using Welch's +#' t-test, which does not assume equal variances. It calculates various statistics, +#' including the Minimum Detectable Difference as a percentage of the control mean. +#' +#' @param data A data frame containing the data. +#' @param factor_col String name of the column containing the grouping factor (e.g., "dose"). +#' @param response_col String name of the column containing the continuous response variable. +#' @param control_level The level of the factor to be used as the control group. +#' If NULL (default), the first factor level is used. +#' @param p.adjust.method Method for adjusting p-values for multiple comparisons. +#' See `?p.adjust` for options (default: "holm"). +#' @param alternative A string specifying the alternative hypothesis. Must be one of +#' "two.sided" (default), "greater", or "less". +#' @param alpha The significance level used to calculate the critical value and MDD%. +#' Default is 0.05. +#' +#' @return A data frame with results for each comparison, including estimates, +#' t-statistic, standard error, critical value, MDD%, raw p-value, +#' and adjusted p-value. +#' @seealso [drcHelper::compare_to_control_fisher()] +#' @export +#' @importFrom stats aggregate p.adjust qt t.test +#' @concept NOEC +#' +#' @examples +#' # Generate example data +#' set.seed(42) +#' my_data <- data.frame( +#' treatment = rep(c("Control", "Dose1", "Dose2", "Dose3"), each = 5), +#' response = c(rnorm(5, 10, 1.5), rnorm(5, 9, 2), +#' rnorm(5, 11, 1.8), rnorm(5, 13, 2.2)) +#' ) +#' +#' # Run Welch's t-test comparing each dose to "Control" +#' welch_results <- compare_to_control_welch( +#' data = my_data, +#' factor_col = "treatment", +#' response_col = "response", +#' control_level = "Control" +#' ) +#' +#' print(welch_results) + +compare_to_control_welch <- function(data, factor_col, response_col, + control_level = NULL, + p.adjust.method = "holm", + alternative = "two.sided", + alpha = 0.05) { + + # --- 1. Validate Inputs --- + if (!all(c(factor_col, response_col) %in% colnames(data))) { + stop("One or more specified columns not found in the data frame.") + } + if (!is.numeric(data[[response_col]])) { + stop(paste("Response column '", response_col, "' must be numeric.")) + } + + factor_data <- data[[factor_col]] + + if (is.null(control_level)) { + control_level <- if(is.factor(factor_data)) levels(factor_data)[1] else sort(unique(factor_data))[1] + message(paste("`control_level` not specified. Using '", control_level, "' as the control.")) + } + + if (!control_level %in% factor_data) { + stop(paste("Control level '", control_level, "' not found in the factor column.")) + } + + # --- 2. Pre-calculate Group Statistics --- + # *** THIS IS THE CORRECTED LINE *** + # The formula object is now the first argument, and it is not named. + group_stats <- aggregate( + as.formula(paste(response_col, "~", factor_col)), + data = data, + FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)) + ) + + group_summary <- do.call(data.frame, group_stats) + colnames(group_summary) <- c(factor_col, "mean", "sd", "n") + + # --- 3. Prepare for Iteration --- + control_data <- data[[response_col]][factor_data == control_level] + control_stats <- group_summary[group_summary[[factor_col]] == control_level, ] + + if(control_stats$mean == 0){ + warning("Control group mean is zero. MDD% will be calculated as Inf or NaN.") + } + + test_levels <- unique(as.character(factor_data)) + test_levels <- test_levels[test_levels != control_level] + + results_list <- list() + + # --- 4. Loop Through Each Test Group and Perform Test --- + for (level in test_levels) { + test_data <- data[[response_col]][factor_data == level] + test_stats <- group_summary[group_summary[[factor_col]] == level, ] + + test_result <- t.test( + x = test_data, + y = control_data, + alternative = alternative, + var.equal = FALSE + ) + + se_diff <- sqrt(test_stats$sd^2 / test_stats$n + control_stats$sd^2 / control_stats$n) + df <- test_result$parameter + t_critical <- if (alternative == "two.sided") { + qt(1 - alpha / 2, df) + } else { + qt(1 - alpha, df) + } + + mdd_absolute <- t_critical * se_diff + mdd_percent <- (mdd_absolute / abs(control_stats$mean)) * 100 + + results_list[[level]] <- data.frame( + p_value = test_result$p.value, + estimate_diff = test_result$estimate[1] - test_result$estimate[2], + statistic = test_result$statistic, + critical_value = t_critical, + se_diff = se_diff, + `MDD%` = mdd_percent, + check.names = FALSE + ) + } + + # --- 5. Combine, Adjust p-values, and Format Output --- + if (length(results_list) == 0) { + message("No test groups to compare against the control.") + return(data.frame()) + } + + final_results <- do.call(rbind, results_list) + final_results[[factor_col]] <- rownames(final_results) + rownames(final_results) <- NULL + + final_results$p_adjusted <- p.adjust(final_results$p_value, method = p.adjust.method) + + final_results <- final_results[, c( + factor_col, "estimate_diff", "statistic", "critical_value", "se_diff", + "MDD%", "p_value", "p_adjusted" + )] + + numeric_cols <- sapply(final_results, is.numeric) + final_results[numeric_cols] <- lapply(final_results[numeric_cols], round, 4) + + return(final_results) +} diff --git a/R/quantal.R b/R/quantal_tests.R similarity index 100% rename from R/quantal.R rename to R/quantal_tests.R diff --git a/man/Tarone.test.Rd b/man/Tarone.test.Rd index 04a6684..397973a 100644 --- a/man/Tarone.test.Rd +++ b/man/Tarone.test.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/quantal.R +% Please edit documentation in R/quantal_tests.R \name{Tarone.test} \alias{Tarone.test} \title{Tarone's Z Test} diff --git a/man/compare_to_control_fisher.Rd b/man/compare_to_control_fisher.Rd index 2316ad6..f765658 100644 --- a/man/compare_to_control_fisher.Rd +++ b/man/compare_to_control_fisher.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/quantal.R +% Please edit documentation in R/quantal_tests.R \name{compare_to_control_fisher} \alias{compare_to_control_fisher} \title{Perform Fisher's Exact Test Comparing Each Level to Control} diff --git a/man/compare_to_control_welch.Rd b/man/compare_to_control_welch.Rd new file mode 100644 index 0000000..9cab9dd --- /dev/null +++ b/man/compare_to_control_welch.Rd @@ -0,0 +1,68 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/continuous_tests.R +\name{compare_to_control_welch} +\alias{compare_to_control_welch} +\title{Perform Many-to-One Welch's t-Tests} +\usage{ +compare_to_control_welch( + data, + factor_col, + response_col, + control_level = NULL, + p.adjust.method = "holm", + alternative = "two.sided", + alpha = 0.05 +) +} +\arguments{ +\item{data}{A data frame containing the data.} + +\item{factor_col}{String name of the column containing the grouping factor (e.g., "dose").} + +\item{response_col}{String name of the column containing the continuous response variable.} + +\item{control_level}{The level of the factor to be used as the control group. +If NULL (default), the first factor level is used.} + +\item{p.adjust.method}{Method for adjusting p-values for multiple comparisons. +See \code{?p.adjust} for options (default: "holm").} + +\item{alternative}{A string specifying the alternative hypothesis. Must be one of +"two.sided" (default), "greater", or "less".} + +\item{alpha}{The significance level used to calculate the critical value and MDD\%. +Default is 0.05.} +} +\value{ +A data frame with results for each comparison, including estimates, +t-statistic, standard error, critical value, MDD\%, raw p-value, +and adjusted p-value. +} +\description{ +Compares each level of a factor to a specified control level using Welch's +t-test, which does not assume equal variances. It calculates various statistics, +including the Minimum Detectable Difference as a percentage of the control mean. +} +\examples{ +# Generate example data +set.seed(42) +my_data <- data.frame( + treatment = rep(c("Control", "Dose1", "Dose2", "Dose3"), each = 5), + response = c(rnorm(5, 10, 1.5), rnorm(5, 9, 2), + rnorm(5, 11, 1.8), rnorm(5, 13, 2.2)) +) + +# Run Welch's t-test comparing each dose to "Control" +welch_results <- compare_to_control_welch( + data = my_data, + factor_col = "treatment", + response_col = "response", + control_level = "Control" +) + +print(welch_results) +} +\seealso{ +\code{\link[=compare_to_control_fisher]{compare_to_control_fisher()}} +} +\concept{NOEC} diff --git a/man/create_contingency_table.Rd b/man/create_contingency_table.Rd index c79342c..8f28440 100644 --- a/man/create_contingency_table.Rd +++ b/man/create_contingency_table.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/quantal.R +% Please edit documentation in R/quantal_tests.R \name{create_contingency_table} \alias{create_contingency_table} \title{Create Contingency Table from Count Data} diff --git a/man/many_to_one_fisher_test.Rd b/man/many_to_one_fisher_test.Rd index 76d59ed..ac5f74d 100644 --- a/man/many_to_one_fisher_test.Rd +++ b/man/many_to_one_fisher_test.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/quantal.R +% Please edit documentation in R/quantal_tests.R \name{many_to_one_fisher_test} \alias{many_to_one_fisher_test} \title{Many-to-One Pairwise Fisher's Exact Test} diff --git a/man/report_dunnett_summary.Rd b/man/report_dunnett_summary.Rd index 5f74c4b..c5141ca 100644 --- a/man/report_dunnett_summary.Rd +++ b/man/report_dunnett_summary.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/MDD.R \name{report_dunnett_summary} \alias{report_dunnett_summary} -\title{Generate a Comprehensive Dunnett Test Summary Report (Robust Version)} +\title{Generate a Comprehensive Dunnett Test Summary Report} \usage{ report_dunnett_summary( formula, @@ -14,19 +14,64 @@ report_dunnett_summary( ) } \arguments{ -\item{formula}{The model formula, e.g., Response ~ Dose.} +\item{formula}{A model formula, e.g., \code{Response ~ Dose} for a fixed-effects model, +or \code{Response ~ Dose + (1|Batch)} for a mixed-effects model.} -\item{data}{The dataframe containing the data.} +\item{data}{A data frame containing the data specified in the formula.} -\item{dunnett_test_func}{The dunnett_test function you use.} +\item{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 \code{glht}. +See \code{multcomp::glht} for details.} -\item{alternative}{The alternative hypothesis to be used for the test.} +\item{alternative}{A character string specifying the alternative hypothesis. Must be one +of "two.sided", "less" (default), or "greater".} -\item{...}{Additional arguments passed to your dunnett_test_func.} +\item{include_random_effect}{A logical value. Set to \code{TRUE} if the \code{formula} includes +a random effect term (e.g., \code{(1|Batch)}), which requires a mixed-effects model. +Defaults to \code{FALSE}.} + +\item{...}{Additional arguments to be passed to the function specified in +\code{dunnett_test_func}. A common argument is \code{alpha} to set the significance +level (e.g., \code{alpha = 0.05}).} } \value{ -A single tibble combining descriptive stats, test results, and MDD\%. +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. } \description{ -Generate a Comprehensive Dunnett Test Summary Report (Robust Version) +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. +} +\examples{ +\dontrun{ +# Assuming `dunnett_test` and helper functions are defined in your environment +# and the 'multcomp' package is loaded. + +# 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) +} +} +\seealso{ +\code{\link[multcomp:glht]{multcomp::glht()}} } diff --git a/vignettes/articles/A-Notes-on-Power.Rmd b/vignettes/articles/A-Notes-on-Power.Rmd new file mode 100644 index 0000000..6cf131b --- /dev/null +++ b/vignettes/articles/A-Notes-on-Power.Rmd @@ -0,0 +1,142 @@ +--- +title: "A Notes on Power" +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(drcHelper) +``` + + +## Concepts + +In statistical analysis, tests can be categorized as liberal or conservative based on their propensity to reject the null hypothesis. Here are examples of both types of tests: + +### Liberal Tests + +1. **Least Significant Difference (LSD) Test**: + - Used for pairwise comparisons following ANOVA. + - Does not adjust for multiple comparisons, making it more prone to Type I errors. + +2. **Unadjusted t-tests for Multiple Comparisons**: + - Performing multiple t-tests without any correction for multiple testing increases the chance of Type I errors. + +3. **Chi-Square Test Without Yates' Continuity Correction**: + - In small sample sizes, not applying Yates' continuity correction can make the test more liberal. + +### Conservative Tests + +1. **Bonferroni Correction**: + - Adjusts the significance level for multiple comparisons by dividing the alpha level by the number of tests. + - Reduces the likelihood of Type I errors but can be overly conservative, increasing the risk of Type II errors. + +2. **Tukey's Honest Significant Difference (HSD) Test**: + - Used for multiple comparisons following ANOVA. + - Controls the family-wise error rate, making it more conservative than LSD. + +3. **Holm-Bonferroni Method**: + - A stepwise approach to control the family-wise error rate. + - Less conservative than Bonferroni but still provides strong control over Type I errors. + +4. **Yates' Continuity Correction for Chi-Square Tests**: + - Applied to 2x2 contingency tables to make the test more conservative, especially with small sample sizes. + +### Implications of Choosing Liberal vs. Conservative Tests + +- **Liberal Tests**: These are useful when the primary concern is maximizing the power to detect differences, and the cost of Type I errors is low. However, they should be used with caution when multiple comparisons are involved. + +- **Conservative Tests**: These are preferred when controlling for Type I errors is crucial, such as in confirmatory studies where false positives could lead to incorrect scientific conclusions. However, they may increase the risk of Type II errors, potentially missing true effects. + +Choosing between liberal and conservative tests depends on the research context, study design, and the balance between Type I and Type II error risks. Researchers should carefully consider these factors when selecting statistical methods for their analyses. + + +## A Powerful Test + +No, a more powerful test is not necessarily the same as a more liberal test. These terms refer to different aspects of statistical testing: + +### Power of a Test + +- **Definition**: The power of a statistical test is the probability that the test correctly rejects a false null hypothesis (i.e., it detects an effect when there is one). Higher power means a greater ability to detect true effects. +- **Factors Influencing Power**: Larger sample sizes, larger effect sizes, and lower variability increase the power of a test. The significance level (alpha) also affects power; a higher alpha can increase power but also increases the risk of Type I errors. + +### Liberal Test + +- **Definition**: A liberal test is one that has a higher probability of rejecting the null hypothesis, which can lead to more Type I errors (false positives). It is less stringent in terms of controlling for these errors. +- **Characteristics**: Typically, liberal tests have a higher Type I error rate than the nominal level, making them more prone to finding statistically significant results even when there is no true effect. + +### Differences Between Power and Liberalism + +1. **Power**: + - Focuses on the test's ability to detect true effects (avoiding Type II errors, or false negatives). + - A test can be powerful without being liberal if it maintains the correct Type I error rate while maximizing the ability to detect true effects. + +2. **Liberalism**: + - Focuses on the test's tendency to reject the null hypothesis, potentially at the expense of increasing Type I errors. + - A liberal test may appear more "powerful" in terms of rejecting the null hypothesis, but this can be misleading because it might also reject the null when it is true. + +### Conclusion + +While both power and liberalism relate to the outcomes of hypothesis testing, they address different types of errors and have different implications for statistical decision-making. Ideally, a test should be both powerful and correctly control the Type I error rate, achieving a balance between detecting true effects and avoiding false positives. + + +## Nonparametric Tests + +### Many-to-one Wilcoxon v.s. Dunn + +The many-to-one version of Dunn's test is more powerful compaired to multiple wilcoxon rank sum test because it maintains the relationship between all groups throughout the analysis by using ranks from the complete dataset to compute test statistics. This integrated approach provides better statistical efficiency compared to the separate pairwise comparisons in Wilcoxon tests. + +1. uses complete ranking information rather than pairwise rankings. +2. more efficiently uses the data structure by considering all groups simultaneously +3. reduced multiple comparison penalty in many-to-one setup +4. Dunn's test handles tied ranks more efficiently across the entire dataset, while Wilcoxon rank sum test processes ties separately for each comparison. + +This reminds me of how ANOVA is generally more powerful than multiple t-tests - it's a similar principle of using all available information simultaneously rather than breaking it into pieces. + + +```{r} +library(PMCMRplus) +## Data set PlantGrowth +## Global test + +prelimPlot1(PlantGrowth %>% mutate(Dose=group,Response=weight)) +kruskalTest(weight ~ group, data = PlantGrowth) + +## Conover's many-one comparison test +## single-step means p-value from multivariate t distribution +ans <- kwManyOneConoverTest(weight ~ group, data = PlantGrowth, + p.adjust.method = "single-step") +summary(ans) + +## Conover's many-one comparison test +ans <- kwManyOneConoverTest(weight ~ group, data = PlantGrowth, + p.adjust.method = "holm") +summary(ans) + +## Dunn's many-one comparison test +ans <- kwManyOneDunnTest(weight ~ group, data = PlantGrowth, + p.adjust.method = "holm") +summary(ans) + +## Nemenyi's many-one comparison test +ans <- kwManyOneNdwTest(weight ~ group, data = PlantGrowth, + p.adjust.method = "holm") +summary(ans) + +## Many one U test +ans <- manyOneUTest(weight ~ group, data = PlantGrowth, + p.adjust.method = "holm") +summary(ans) + +## Chen Test +ans <- chenTest(weight ~ group, data = PlantGrowth, + p.adjust.method = "holm") +summary(ans) +``` + + diff --git a/vignettes/articles/MDD-in-Regulatory-Context.Rmd b/vignettes/articles/MDD-in-Regulatory-Context.Rmd index a7896fa..76cf028 100644 --- a/vignettes/articles/MDD-in-Regulatory-Context.Rmd +++ b/vignettes/articles/MDD-in-Regulatory-Context.Rmd @@ -58,281 +58,10 @@ Notes - Computes MDD = Tcrit × SE_diff and MDD% = 100 × MDD / control_mean. - Returns a tidy tibble with dose, n_c, n_t, df (when applicable), SE_diff, Tcrit, MDD, MDD_pct, method, alpha, alternative, and control_mean. -## Code: compute_mdd_generic +## Usage examples -```{r} -compute_mdd_generic <- function(x, data, - test = c("Williams_PMCMRplus", "Williams_JG", "Dunnett", "t", "Welch"), - alpha = 0.05, - alternative = c("smaller", "greater", "two-sided")) { - test <- match.arg(test) - alternative <- match.arg(alternative) - - # Parse formula - if (!inherits(x, "formula")) stop("x must be a formula, e.g., Response ~ Dose") - terms_obj <- terms(x) - resp_name <- rownames(attr(terms_obj, "factors"))[1] - factor_name <- attr(terms_obj, "term.labels")[1] - - df <- data - if (!is.data.frame(df)) stop("data must be a data.frame") - - # Prepare dose factor with control first (lowest dose) - dose_raw <- df[[factor_name]] - dose_num <- suppressWarnings(as.numeric(gsub(",", ".", as.character(dose_raw)))) - # If numeric conversion fails everywhere, treat as factor and try to sort levels by numeric content - if (all(is.na(dose_num))) { - # Attempt to parse numeric part from strings - dose_num <- suppressWarnings(as.numeric(gsub(",", ".", gsub("[^0-9.,-]", "", as.character(dose_raw))))) - } - # Build factor levels: sort unique numeric doses; fallback to existing factor levels - if (!all(is.na(dose_num))) { - levels_sorted <- sort(unique(dose_num)) - df$Dose_factor <- factor(dose_num, levels = levels_sorted) - } else { - df$Dose_factor <- if (is.factor(dose_raw)) factor(dose_raw) else factor(dose_raw) - warning("Dose values not numeric; using existing factor levels. Ensure control is the first level.") - } - - # Basic stats per group - resp <- df[[resp_name]] - ctrl_level <- levels(df$Dose_factor)[1] - ctrl_idx <- df$Dose_factor == ctrl_level - mu_c <- mean(resp[ctrl_idx], na.rm = TRUE) - n_c <- sum(ctrl_idx) - - # One-way ANOVA for pooled MSE (used by Williams/Dunnett/t pooled) - aov_fit <- stats::aov(resp ~ df$Dose_factor) - mse <- tryCatch(summary(aov_fit)[[1]]["Residuals", "Mean Sq"], error = function(e) NA_real_) - df_resid <- tryCatch(summary(aov_fit)[[1]]["Residuals", "Df"], error = function(e) NA_real_) - - # Per-dose counts and SDs - group_n <- aggregate(resp ~ df$Dose_factor, FUN = function(z) sum(!is.na(z))) - group_sd <- aggregate(resp ~ df$Dose_factor, FUN = stats::sd) - names(group_n) <- c("Dose_factor", "n_t") - names(group_sd) <- c("Dose_factor", "sd_t") - - # Helper to compute one-/two-sided quantile tail - qtail <- function(df, alt, alpha) { - if (alt == "two-sided") stats::qt(1 - alpha / 2, df = df) else stats::qt(1 - alpha, df = df) - } - - # Initialize result container - res_list <- list() - - if (test %in% c("Williams_PMCMRplus", "Williams_JG")) { - # Direction for Williams - direction <- if (alternative == "smaller") "decreasing" else "increasing" - bw <- try(drcHelper::broom_williams(x, data = df, method = if (test == "Williams_PMCMRplus") "Williams_PMCMRplus" else "Williams_JG", direction = direction), silent = TRUE) - if (inherits(bw, "try-error") || nrow(bw) == 0) stop("broom_williams failed to produce results.") - - bw <- as.data.frame(bw) - # Map Dose from comparison ("a - 0 <= 0" or "a - 0 >= 0") - comp_clean <- gsub("\\s*(<=|>=)\\s*0\\s*$", "", as.character(bw$comparison)) - dose_contrast <- suppressWarnings(as.numeric(gsub(",", ".", sub(" - 0$", "", comp_clean)))) - if (all(is.na(dose_contrast)) && test == "Williams_JG") { - # Fallback: assign ascending treatment doses excluding control - trt_levels <- levels(df$Dose_factor)[-1] - dose_contrast <- suppressWarnings(as.numeric(trt_levels)) - dose_contrast <- dose_contrast[seq_len(min(length(dose_contrast), nrow(bw)))] - # If lengths mismatch, pad with NA - if (length(dose_contrast) < nrow(bw)) { - dose_contrast <- c(dose_contrast, rep(NA_real_, nrow(bw) - length(dose_contrast))) - } - } - names(bw)[names(bw) == "`t'-crit"] <- "Tcrit" - # Merge per-dose counts - trt_n <- merge(data.frame(Dose_factor = factor(dose_contrast, levels = levels(df$Dose_factor))), - group_n, by = "Dose_factor", all.x = TRUE) - n_t <- trt_n$n_t - # SE_diff via pooled MSE - if (is.na(mse)) stop("Cannot compute pooled MSE from ANOVA; MDD cannot be computed.") - SE_diff <- sqrt(mse * (1 / n_c + 1 / n_t)) # $$SE_{diff} = \sqrt{MSE \cdot (1/n_c + 1/n_t)}$$ - MDD <- bw$Tcrit * SE_diff # $$MDD = T_{crit} \cdot SE_{diff}$$ - MDD_pct <- if (!is.na(mu_c) && mu_c != 0) 100 * MDD / mu_c else NA_real_ # $$MDD\% = 100 \cdot MDD / \mu_c$$ - res_list <- list( - dose = dose_contrast, - n_c = rep(n_c, length(n_t)), - n_t = n_t, - df = rep(df_resid, length(n_t)), - SE_diff = SE_diff, - Tcrit = bw$Tcrit, - MDD = MDD, - MDD_pct = MDD_pct, - method = test - ) - } else if (test == "Dunnett") { - # Need multcomp - if (!requireNamespace("multcomp", quietly = TRUE)) stop("Package 'multcomp' is required for Dunnett MDD.") - # Build Dunnett contrasts and glht - contrasts_call <- as.call(list(as.name("mcp"))) - contrasts_call[[2]] <- "Dunnett" - names(contrasts_call)[2] <- factor_name - dunnett_contrasts <- eval(contrasts_call) - # Construct aov from formula x - aov_model <- stats::aov(x, data = df) - gl <- multcomp::glht(aov_model, linfct = dunnett_contrasts) - sm <- summary(gl) - # Dose from rownames of coefficients like "0.0448 - 0" - comp <- names(sm$test$coefficients) - dose_contrast <- suppressWarnings(as.numeric(gsub(",", ".", sub(" - 0$", "", comp)))) - # Standard errors per contrast come from sigma in summary (vector) - SE_diff <- as.numeric(sm$test$sigma) - # Critical value from qfunction - qfun <- sm$test$qfunction - Tcrit <- tryCatch(as.numeric(qfun(conf.level = 1 - alpha, adjusted = TRUE)), error = function(e) NA_real_) - if (is.na(Tcrit)) { - warning("Could not obtain Dunnett adjusted critical value; using unadjusted qt with residual df.") - Tcrit <- qtail(df_resid, alternative, alpha) - } - MDD <- Tcrit * SE_diff - MDD_pct <- if (!is.na(mu_c) && mu_c != 0) 100 * MDD / mu_c else NA_real_ - # n_t per contrast - trt_n <- merge(data.frame(Dose_factor = factor(dose_contrast, levels = levels(df$Dose_factor))), - group_n, by = "Dose_factor", all.x = TRUE) - n_t <- trt_n$n_t - res_list <- list( - dose = dose_contrast, - n_c = rep(n_c, length(n_t)), - n_t = n_t, - df = rep(df_resid, length(SE_diff)), - SE_diff = SE_diff, - Tcrit = rep(Tcrit, length(SE_diff)), - MDD = MDD, - MDD_pct = MDD_pct, - method = "Dunnett_multcomp" - ) - } else if (test == "t") { - # Pooled-variance per contrast (control vs each treatment) - # Compute control SD - sd_c <- stats::sd(resp[ctrl_idx]) - trt_levels <- levels(df$Dose_factor)[-1] - out_rows <- lapply(trt_levels, function(lv) { - trt_idx <- df$Dose_factor == lv - n_t <- sum(trt_idx) - sd_t <- stats::sd(resp[trt_idx]) - # pooled variance - sp2 <- ((n_c - 1) * sd_c^2 + (n_t - 1) * sd_t^2) / (n_c + n_t - 2) # $$s_p^2 = ((n_c-1)s_c^2 + (n_t-1)s_t^2) / (n_c+n_t-2)$$ - SE_diff <- sqrt(sp2 * (1 / n_c + 1 / n_t)) # $$SE_{diff} = \sqrt{s_p^2 (1/n_c + 1/n_t)}$$ - df_ct <- n_c + n_t - 2 - Tcrit <- qtail(df_ct, alternative, alpha) - MDD <- Tcrit * SE_diff - MDD_pct <- if (!is.na(mu_c) && mu_c != 0) 100 * MDD / mu_c else NA_real_ - data.frame(dose = suppressWarnings(as.numeric(lv)), n_c = n_c, n_t = n_t, - df = df_ct, SE_diff = SE_diff, Tcrit = Tcrit, MDD = MDD, MDD_pct = MDD_pct, - method = "Student_t") - }) - res_df <- do.call(rbind, out_rows) - return(res_df) - } else if (test == "Welch") { - # Welch SE and df per contrast - sd_c <- stats::sd(resp[ctrl_idx]) - trt_levels <- levels(df$Dose_factor)[-1] - out_rows <- lapply(trt_levels, function(lv) { - trt_idx <- df$Dose_factor == lv - n_t <- sum(trt_idx) - sd_t <- stats::sd(resp[trt_idx]) - SE_diff <- sqrt(sd_c^2 / n_c + sd_t^2 / n_t) # $$SE_{diff} = \sqrt{s_c^2/n_c + s_t^2/n_t}$$ - df_w <- (sd_c^2 / n_c + sd_t^2 / n_t)^2 / - ((sd_c^2 / n_c)^2 / (n_c - 1) + (sd_t^2 / n_t)^2 / (n_t - 1)) # $$df_{Welch} = \frac{(s_c^2/n_c + s_t^2/n_t)^2}{(s_c^2/n_c)^2/(n_c-1) + (s_t^2/n_t)^2/(n_t-1)}$$ - Tcrit <- qtail(df_w, alternative, alpha) - MDD <- Tcrit * SE_diff - MDD_pct <- if (!is.na(mu_c) && mu_c != 0) 100 * MDD / mu_c else NA_real_ - data.frame(dose = suppressWarnings(as.numeric(lv)), n_c = n_c, n_t = n_t, - df = df_w, SE_diff = SE_diff, Tcrit = Tcrit, MDD = MDD, MDD_pct = MDD_pct, - method = "Welch_t") - }) - res_df <- do.call(rbind, out_rows) - return(res_df) - } - - # Assemble and return tibble for Williams/Dunnett branches - res_df <- data.frame( - Dose = res_list$dose, - n_c = res_list$n_c, - n_t = res_list$n_t, - df = res_list$df, - SE_diff = res_list$SE_diff, - Tcrit = res_list$Tcrit, - MDD = res_list$MDD, - MDD_pct = res_list$MDD_pct, - method = res_list$method, - alpha = alpha, - alternative = alternative, - control_mean = mu_c - ) - # Order by Dose ascending - res_df <- res_df[order(res_df$Dose), ] - res_df -} - -``` - -Usage examples - -- Williams PMCMRplus (your MOCK0065 Growth Rate): - -```{r eval=FALSE} -# Helper to convert dose strings to numeric, handling various formats -convert_dose <- function(x) { - if (length(x) == 0) return(numeric(0)) - xc <- as.character(x) - xc <- trimws(xc) - xc[xc %in% c("", "n/a", "NA")] <- NA_character_ - # Normalize decimal commas and keep scientific notation (e.g., "4,48E-2" -> "4.48E-2") - xc <- gsub(",", ".", xc, fixed = TRUE) - out <- suppressWarnings(as.numeric(xc)) - return(out) -} - -convert_numeric <- function(x) { - if (length(x) == 0) return(numeric(0)) - xc <- as.character(x) - xc <- trimws(xc) - xc[xc %in% c("", "n/a", "NA")] <- NA_character_ - xc <- gsub(",", ".", xc, fixed = TRUE) - suppressWarnings(as.numeric(xc)) -} -dose_from_comparison <- function(comp_vec) { - if (length(comp_vec) == 0) return(numeric(0)) - vapply(comp_vec, function(s) { - if (is.na(s)) return(NA_real_) - parts <- strsplit(s, " - ", fixed = TRUE)[[1]] - convert_dose(parts[1]) - }, FUN.VALUE = numeric(1)) -} -w_ep <- test_cases_data %>% dplyr::filter(`Study ID` == "MOCK0065", Endpoint == "Growth Rate") -w_ep$Dose_numeric <- convert_dose(w_ep$Dose) -w_ep$Dose_factor <- factor(w_ep$Dose_numeric, levels = sort(unique(w_ep$Dose_numeric))) -res_william <- broom_williams(Response ~ Dose_factor, data = w_ep, - test = "Williams_PMCMRplus", alpha = 0.05, alternative = "less") -mdd_w_pm <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, - test = "Williams_PMCMRplus", alpha = 0.05, alternative = "smaller") -print(mdd_w_pm) -``` - - -- Dunnett: -```{r eval=FALSE} -mdd_dunn <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, - test = "Dunnett", alpha = 0.05, alternative = "smaller") -print(mdd_dunn) -``` -```{r eval=FALSE} -- Student’s t (pooled, control vs each dose): - -mdd_t <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, - test = "t", alpha = 0.05, alternative = "smaller") -print(mdd_t) -``` -- Welch: -```{r eval=FALSE} -mdd_welch <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, - test = "Welch", alpha = 0.05, alternative = "two-sided") -print(mdd_welch) -``` -Notes and options +### Notes and options - Alternative handling: MDD uses a critical value magnitude; the choice of one-sided vs two-sided changes the quantile. If your expected MDD% was computed with two-sided thresholds, set alternative = "two-sided". - For Dunnett, the adjusted Tcrit is taken from `summary(glht)$test$qfunction`. If unavailable, we fallback to unadjusted qt with residual df. - For `Williams_JG`, comparison strings might not carry doses; we assign doses in ascending order excluding control. If you need exact mapping, we can inspect the internal `williamsTest_JG` output via drcHelper to align rows to doses. @@ -343,6 +72,9 @@ Notes and options Reusing existing test objects and building a modular system around them is far more efficient and robust. + +## Williams' Test MDD + ```{r} library(drcHelper) Results <- broom_williams(Response ~ factor(Dose), data = dat_medium, @@ -352,6 +84,7 @@ compute_mdd_williams(Results,data=dat_medium,formula = Response ~ factor(Dose)) Results <- dunnett_test(data = dat_medium,response_var = "Response",tank_var = NULL,include_random_effect = FALSE,dose_var = "Dose", alpha = 0.05, alternative = "less") ``` +## Dunnett's Test MDD ```{r} Results <- broom_dunnett(Response ~ factor(Dose), data = dat_medium, @@ -405,7 +138,98 @@ final_dunnett_table <- report_dunnett_summary( print(final_dunnett_table) ``` +## Welch's t-test (Multiple Comparison) + +Before we explain the MDD calculation for Welch's t-test implemented in `drcHelper`, we need to explain the difference between test sensitivity (power) and error rate control. + +### 1. What MDD Represents: The Power of a *Single* Test + +The Minimum Detectable Difference (MDD) is a measure of statistical power. It answers the question: + +> "For a **single pairwise comparison**, given its specific sample sizes and data variability, what is the smallest true difference between means that I would have a high probability of detecting as statistically significant at a given alpha level (e.g., $\alpha = 0.05$)?" + +The calculation reflects this focus on a single test: +$MDD = t_{critical} \times SE_{diff}$ + +* **$SE_{diff} $(Standard Error of the Difference):** This depends only on the standard deviations and sample sizes of the two groups being compared. It's a property of the data. +* **$t_{critical} $(Critical Value):** This is derived from the t-distribution using two parameters: + 1. The degrees of freedom (determined by the sample sizes and variances). + 2. The chosen significance level for that **single test** (`alpha`, usually 0.05). + +Notice that nowhere in this calculation is there a parameter for the *total number of tests being performed*. The MDD is an intrinsic property of the individual test's design and sensitivity. + +### 2. What P-Value Adjustment Represents: Controlling Error Across *Multiple* Tests +P-value adjustment methods (like Holm or Bonferroni) are designed to solve a different problem: controlling the **Family-Wise Error Rate (FWER)**. They answer the question: + +> "Now that I have performed **many tests at once**, how can I adjust my significance threshold to ensure the probability of making even one false positive (Type I error) across the *entire family* of tests remains low (e.g., at 0.05)?" + +These methods work by making it "harder" for any individual test to be declared significant. They do this by either: +* Decreasing the p-value threshold for significance (e.g., Bonferroni's $\alpha_{adj} = \alpha / k$). +* Increasing the raw p-values to create `p_adjusted` values, which are then compared to the original `alpha` of 0.05. + +The key takeaway is that p-value adjustment is a **post-hoc correction** applied to the *results* of the tests, not a change to the *design or inherent power* of the individual tests themselves. + +### Could You Calculate an "Adjusted MDD"? + +Theoretically, yes. You could calculate a more stringent MDD by using an adjusted critical value based on a Bonferroni-corrected alpha (e.g., $\alpha_{adj} = 0.05 / k$). This would result in a larger $t_{critical} $and therefore a larger MDD. + +However, this is **not standard practice** and can be confusing. The value of the MDD is in understanding the sensitivity of your experimental design for a given pairwise comparison. It tells you the inherent power of your test. Mixing in the multi-comparison adjustment muddies this clear interpretation. + +**Conclusion:** The function `compare_to_control_welch` implements the standard and most useful approach: + +1. It calculates the **MDD** based on the properties of each **individual t-test** (using the single-test `alpha` of 0.05). This tells you about the sensitivity of your experiment. +2. It then **separately** calculates the **adjusted p-values**. This provides a corrected measure for making final conclusions about significance while controlling the overall error rate. + +```{r} +# Generate example data +set.seed(42) +my_data <- data.frame( + treatment = rep(c("Control", "Dose1", "Dose2", "Dose3"), each = 5), + response = c(rnorm(5, 10, 1.5), rnorm(5, 9, 2), + rnorm(5, 11, 1.8), rnorm(5, 13, 2.2)) +) + +# Run Welch's t-test comparing each dose to "Control" +welch_results <- compare_to_control_welch( + data = my_data, + factor_col = "treatment", + response_col = "response", + control_level = "Control" +) + +print(welch_results) +``` + +### Run two-sample Welch's t-test + +```{r} +## Two sample test! +welch_results <- compare_to_control_welch( + data = my_data%>%filter(treatment %in% c("Control","Dose1")), + factor_col = "treatment", + response_col = "response", + control_level = "Control" +) +print(welch_results) +``` + + +```{r} +# Run Welch's t-test comparing each dose to "Control" +welch_results <- compare_to_control_welch( + data = dat_medium, + factor_col = "Dose", + response_col = "Response", + control_level = "0" +) + +# print(welch_results) +welch_results %>% knitr::kable(.,digits = 2) +``` + + +## Notes ### Handling "greater" vs "smaller" in Williams' test @@ -414,7 +238,7 @@ print(final_dunnett_table) ### Unequal sample sizes and heteroscedasticity - The SE_diff above already handles unequal n via 1/n_c + 1/n_t. -- If variances are heterogeneous and a Welch-type approach is used, you would estimate SE_diff via group-specific variances. Williams is based on a pooled-variance ANOVA, so the pooled MSE is appropriate here. +- If variances are heterogeneous and a Welch-type approach is used, you would estimate SE_diff via group-specific variances. Williams is based on a pooled-variance ANOVA, so the pooled MSE is appropriate here.We probably need to adjust if the expected MDD% uses a different SE (e.g., contrast-specific SE in Williams instead of the simple (1/n_c + 1/n_t) form). ## Optional: power-based MDD% @@ -423,6 +247,71 @@ print(final_dunnett_table) - Then $$MDD\%_{power} = 100 \cdot MDD_{power} / \mu_c$$ - df can be taken from the residual df of the ANOVA model. You can parameterize β (e.g., 0.2 for 80% power) and compute the quantiles using stats::qt. -## Notes -We probably need to adjust if the expected MDD% uses a different SE (e.g., contrast-specific SE in Williams instead of the simple (1/n_c + 1/n_t) form). + + +# In test cases validation + +- Williams PMCMRplus (your MOCK0065 Growth Rate): + +```{r eval=FALSE} +# Helper to convert dose strings to numeric, handling various formats +convert_dose <- function(x) { + if (length(x) == 0) return(numeric(0)) + xc <- as.character(x) + xc <- trimws(xc) + xc[xc %in% c("", "n/a", "NA")] <- NA_character_ + # Normalize decimal commas and keep scientific notation (e.g., "4,48E-2" -> "4.48E-2") + xc <- gsub(",", ".", xc, fixed = TRUE) + out <- suppressWarnings(as.numeric(xc)) + return(out) +} + +convert_numeric <- function(x) { + if (length(x) == 0) return(numeric(0)) + xc <- as.character(x) + xc <- trimws(xc) + xc[xc %in% c("", "n/a", "NA")] <- NA_character_ + xc <- gsub(",", ".", xc, fixed = TRUE) + suppressWarnings(as.numeric(xc)) +} + +dose_from_comparison <- function(comp_vec) { + if (length(comp_vec) == 0) return(numeric(0)) + vapply(comp_vec, function(s) { + if (is.na(s)) return(NA_real_) + parts <- strsplit(s, " - ", fixed = TRUE)[[1]] + convert_dose(parts[1]) + }, FUN.VALUE = numeric(1)) +} +w_ep <- test_cases_data %>% dplyr::filter(`Study ID` == "MOCK0065", Endpoint == "Growth Rate") +w_ep$Dose_numeric <- convert_dose(w_ep$Dose) +w_ep$Dose_factor <- factor(w_ep$Dose_numeric, levels = sort(unique(w_ep$Dose_numeric))) +res_william <- broom_williams(Response ~ Dose_factor, data = w_ep, + test = "Williams_PMCMRplus", alpha = 0.05, alternative = "less") +mdd_w_pm <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, + test = "Williams_PMCMRplus", alpha = 0.05, alternative = "smaller") +print(mdd_w_pm) +``` + + +- Dunnett: +```{r eval=FALSE} +mdd_dunn <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, + test = "Dunnett", alpha = 0.05, alternative = "smaller") +print(mdd_dunn) +``` +```{r eval=FALSE} +- Student’s t (pooled, control vs each dose): + +mdd_t <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, + test = "t", alpha = 0.05, alternative = "smaller") +print(mdd_t) +``` +- Welch: +```{r eval=FALSE} +mdd_welch <- compute_mdd_generic(Response ~ Dose_factor, data = w_ep, + test = "Welch", alpha = 0.05, alternative = "two-sided") +print(mdd_welch) +``` + diff --git a/vignettes/articles/Quantal-Data.Rmd b/vignettes/articles/Quantal-Data.Rmd index 1e42e7d..658ff2d 100644 --- a/vignettes/articles/Quantal-Data.Rmd +++ b/vignettes/articles/Quantal-Data.Rmd @@ -513,3 +513,93 @@ result <- cochranArmitageTrendTest(successes = as.numeric(second_arcsin_data_2$I result ``` + +## A Note on Implementation + +Of course! I can help you understand the differences between these two R functions. + +There are two wrapper function for performing many-to-one fisher's test: `many_to_one_fisher_test` and `compare_to_control_fisher`. Both functions are designed to solve a similar problem: performing multiple Fisher's exact tests to compare several treatment groups against a single control or reference group. However, they are written with different philosophies, accept different data formats, and use different underlying packages. + +Here is a detailed breakdown of their key differences: + +### 1. Input Data Format + +This is the most significant practical difference for a user. + +* `many_to_one_fisher_test`: **Requires a pre-summarized contingency table** (a `matrix` or `data.frame`). The rows should represent your groups (e.g., `dose_0`, `dose_10`, `dose_50`), and the two columns should represent the counts of the two outcomes (e.g., `alive`, `dead`). You must aggregate your raw data into this table format *before* calling the function. + +* `compare_to_control_fisher`: **Works directly with a "long" format data frame**. This data frame should have a column for your groups (e.g., `dose`), a column for success counts (e.g., `alive`), and a column for failure counts (e.g., `dead`). The function performs the necessary aggregation (summing counts for each group) internally. + +### 2. Dependencies & "Ecosystem" + +The functions are built upon different sets of tools. + +* `many_to_one_fisher_test`: Is built as an extension of the **`rstatix` package** and the wider **Tidyverse**. It uses functions from `rstatix`, `dplyr`, and `purrr`. It is designed to fit seamlessly into a workflow that already uses these packages. + +* `compare_to_control_fisher`: Is more of a **standalone, "base R" style function**. Although it uses a `%>%` pipe for clarity, its core statistical logic relies only on the built-in `stats` package (`fisher.test`, `p.adjust`). It has fewer external dependencies for its core functionality. + +### 3. Core Implementation + +The internal logic reflects their different dependencies. + +* `many_to_one_fisher_test`: + * It loops through the non-reference groups. + * For each group, it creates a 2x2 table and passes it to **`rstatix::fisher_test()`**. + * The `rstatix` function is a wrapper around the base `stats::fisher.test()`, but it formats the output into a tidy tibble. + * It uses `rstatix::adjust_pvalue()` for p-value correction. + +* `compare_to_control_fisher`: + * It manually calculates the sum of successes and failures for the control group. + * It then loops through the test groups, calculating their sums. + * Inside the loop, it manually constructs a 2x2 `matrix`. + * It passes this matrix directly to the base **`stats::fisher.test()`**. + * It collects the results in a new data frame and then uses the base `stats::p.adjust()` for p-value correction. + +### 4. Output Format and Content + +The functions return differently structured results. + +* `many_to_one_fisher_test`: Returns a **`tibble` with the special `rstatix_test` class**. This makes it compatible with other `rstatix` functions. The columns are standardized (`group1`, `group2`, `n`, `p`, `p.adj`). It does not return the odds ratio or confidence interval in the default output. + +* `compare_to_control_fisher`: Returns a **standard `data.frame`**. The output is more detailed by default, including: + * The group level being tested. + * The raw p-value (`p_value`). + * The adjusted p-value (`p_adjusted`). + * The calculated **odds ratio** (`odds_ratio`). + * The **confidence interval** for the odds ratio (`ci_lower`, `ci_upper`). + +### 5. Function Arguments and Flexibility + +* `many_to_one_fisher_test`: + * The reference group **must** be specified via `ref.group`. + * It can pass additional arguments (like `alternative` or `conf.level`) to `fisher.test()` via the `...` parameter. + * The `detailed = TRUE` argument can be used to get more information from the underlying `rstatix::fisher_test` call. + +* `compare_to_control_fisher`: + * Has a smart default for the `control_level` (it chooses the first level if not specified). + * Has dedicated arguments for `alternative` and `conf.level`, making them more explicit. + + +### Summary Table + +| Feature | `many_to_one_fisher_test` | `compare_to_control_fisher` | +| :--- | :--- | :--- | +| **Input Data** | **Contingency Table** (summarized counts) | **Long Data Frame** (raw counts per replicate) | +| **Dependencies** | Tidyverse: `rstatix`, `dplyr`, `purrr` | Primarily `stats` (Base R) | +| **Core Function Call** | `rstatix::fisher_test()` | `stats::fisher.test()` | +| **Output Object** | `rstatix_test` tibble | Standard `data.frame` | +| **Output Content** | `p`, `p.adj` | `p_value`, `p_adjusted`, **odds ratio**, **CI** | +| **Reference Group** | Must be specified | Defaults to first level if not specified | + +### Which One Should You Use? + +* Choose **`many_to_one_fisher_test`** if: + * You are already heavily using the `rstatix` package and Tidyverse workflow. + * Your data is already summarized into a contingency table. + * You prefer the standardized output format of `rstatix`. + +* Choose **`compare_to_control_fisher`** if: + * You have your data in a long format with columns for successes and failures. + * You prefer a self-contained function with fewer package dependencies. + * You want the odds ratio and its confidence interval directly in your results table. + From dae374a7f1f6e877dad1adaa8ff397995c01395b Mon Sep 17 00:00:00 2001 From: Zhenglei Gao Date: Fri, 3 Oct 2025 16:36:27 +0200 Subject: [PATCH 4/5] added noec_from_trend_test function, report_dunnett_summary function, update some vignettes, added jdata for Jockheere-Terpstra test --- NAMESPACE | 2 + NEWS.md | 6 +- R/data_description.R | 23 ++ R/stepDownTrendTest_wrapper.R | 187 ++++++++++ _pkgdown.yml | 1 + data/jdata.rda | Bin 0 -> 249 bytes man/jdata.Rd | 27 ++ man/noec_from_trend_test.Rd | 36 ++ man/print.noecTrendTest.Rd | 16 + man/report_dunnett_summary.Rd | 3 - man/stepDownTrendTest_NOEC.Rd | 37 ++ vignettes/articles/Trend-Testing.Rmd | 62 +++- vignettes/articles/Verification_which_JT.Rmd | 351 ++++++++++++------ .../articles/Verification_which_JT.Rmd.orig | 88 ----- vignettes/articles/Verification_which_JT.md | 165 ++++++++ 15 files changed, 800 insertions(+), 204 deletions(-) create mode 100644 R/stepDownTrendTest_wrapper.R create mode 100644 data/jdata.rda create mode 100644 man/jdata.Rd create mode 100644 man/noec_from_trend_test.Rd create mode 100644 man/print.noecTrendTest.Rd create mode 100644 man/stepDownTrendTest_NOEC.Rd delete mode 100644 vignettes/articles/Verification_which_JT.Rmd.orig create mode 100644 vignettes/articles/Verification_which_JT.md diff --git a/NAMESPACE b/NAMESPACE index 0c3c4ee..c5b6987 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -66,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) @@ -87,6 +88,7 @@ export(simplifyTreatment) export(simulate_dose_response) export(stepDownRSCABS) export(stepDownTrendTestBinom) +export(stepDownTrendTest_NOEC) export(stepKRSCABS) export(step_down_RSCABS) export(summaryZG) diff --git a/NEWS.md b/NEWS.md index 09b441f..f235026 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,8 +5,10 @@ ## 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 -- Added `compute_mdd_williams` and `compute_mdd_dunnett` for MDD calculations, `report_dunnett_summary` for the reporting of test result. +- 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 diff --git a/R/data_description.R b/R/data_description.R index 990a9ab..092c9d3 100644 --- a/R/data_description.R +++ b/R/data_description.R @@ -166,3 +166,26 @@ NULL #' @docType data #' @keywords datasets "metaldata" + + +#' Jonckheere's Synthetic Dataset +#' +#' This dataset is used by Jonckheere (1954) to illustrate a distribution-free k-sample test against ordered alternatives. +#' +#' @format A data frame with 16 rows and 2 columns: +#' \describe{ +#' \item{X}{Number of samples, indicating the group each observation belongs to.} +#' \item{Y}{Value of the observation.} +#' } +#' +#' @source Jonckheere, A. R. (1954). A Distribution-Free k-Sample Test Against Ordered Alternatives. \emph{Biometrika}, 41, 133-145. +#' +#' @usage data(jdata) +#' +#' @examples +#' data(jdata) +#' str(jdata) +#' summary(jdata) +#' +#' @name jdata +NULL diff --git a/R/stepDownTrendTest_wrapper.R b/R/stepDownTrendTest_wrapper.R new file mode 100644 index 0000000..9bb36be --- /dev/null +++ b/R/stepDownTrendTest_wrapper.R @@ -0,0 +1,187 @@ +#' Determine NOEC/LOEC using a Step-down Trend Test +#' +#' This function acts as a wrapper around PMCMRplus::stepDownTrendTest to +#' directly calculate the No-Observed-Effect Concentration (NOEC) and +#' Lowest-Observed-Effect Concentration (LOEC). +#' +#' It interprets the "step-up" ordered results from the underlying PMCMRplus +#' function by searching backwards from the full test to find the first +#' non-significant result (p >= alpha), as is standard for NOEC determination. +#' +#' @param formula A formula of the form `response ~ group`. +#' @param data A data frame containing the variables. The group variable should be a factor. +#' @param test The trend test to use (passed to PMCMRplus). Defaults to "jonckheereTest". +#' @param alpha The significance level to use for determining the NOEC. Defaults to 0.05. +#' @param ... Other arguments passed to `PMCMRplus::stepDownTrendTest`. +#' +#' @return A list of class `noecTrendTest` containing: +#' \item{NOEC}{The determined No-Observed-Effect Concentration.} +#' \item{LOEC}{The determined Lowest-Observed-Effect Concentration.} +#' \item{alpha}{The significance level used.} +#' \item{full_results}{The complete, original `trendPMCMR` object for inspection.} +#' @export +noec_from_trend_test <- function(formula, data, test = "jonckheereTest", alpha = 0.05, ...) { + + # --- Step 1: Input validation --- + group_var_name <- all.vars(formula)[2] + if (!is.factor(data[[group_var_name]])) { + warning(paste("Coercing group variable", group_var_name, "to a factor for correct ordering.")) + data[[group_var_name]] <- factor(data[[group_var_name]]) + } + + # --- Step 2: Run the underlying PMCMRplus test --- + # This gives us all the p-values, albeit in the "step-up" order. + results_obj <- PMCMRplus::stepDownTrendTest( + formula = formula, + data = data, + test = test, + ... + ) + + # --- Step 3: Extract p-values and dose levels --- + p_values <- results_obj$p.value + dose_levels <- rownames(p_values) + control_level <- colnames(p_values)[1] + all_levels <- c(control_level, dose_levels) + + # --- Step 4: Find the NOEC by searching backwards --- + noec <- NA + + # Iterate from the full test (last row) backwards to the first test (first row) + for (i in nrow(p_values):1) { + # Check if this test is NOT significant + if (p_values[i, 1] >= alpha) { + # If it's not significant, the highest dose in this test is the NOEC. + noec <- dose_levels[i] + break # We found the NOEC, so we can stop searching. + } + } + + # --- Step 5: Handle edge cases --- + loec <- NA + + # Case 1: The loop finished without finding a non-significant test. + # This means all tests were significant, so NOEC is the control. + if (is.na(noec)) { + noec <- control_level + loec <- dose_levels[1] # The lowest tested dose is the LOEC + } else { + # Case 2: A NOEC was found. The LOEC is the next dose level up. + noec_index <- which(all_levels == noec) + # Check if the NOEC is the highest tested dose + if (noec_index == length(all_levels)) { + loec <- "Not Detected" + } else { + loec <- all_levels[noec_index + 1] + } + } + + # --- Step 6: Create a clean output object --- + output <- list( + NOEC = noec, + LOEC = loec, + alpha = alpha, + full_results = results_obj + ) + + class(output) <- "noecTrendTest" + + return(output) +} + +#' Print method for noecTrendTest objects +#' @param x An object of class `noecTrendTest`. +#' @param ... Other arguments (not used). +#' @method print noecTrendTest +print.noecTrendTest <- function(x, ...) { + cat("--- Step-Down Trend Test for NOEC Determination ---\n") + cat("Test method:", x$full_results$method, "\n\n") + cat("NOEC:", x$NOEC, "\n") + cat("LOEC:", x$LOEC, "\n") + cat("Alpha level:", x$alpha, "\n\n") + cat("--- Full PMCMRplus Results ---\n") + full_results <- x$full_results + summary(full_results) + invisible(x) +} + + + + + + + + + + + +#' A Wrapper for PMCMRplus::stepDownTrendTest for NOEC Determination +#' +#' This function calls PMCMRplus::stepDownTrendTest but reorders the results +#' to follow a true "step-down" procedure, making it intuitive for finding a +#' NOEC (No-Observed-Effect Concentration). +#' +#' The standard ecotoxicology step-down procedure is: +#' 1. Test for a trend across all treatment groups vs. control. +#' 2. If significant, remove the highest dose and re-test. +#' 3. Repeat until the test is not significant. The highest dose in that +#' non-significant test is the NOEC. +#' +#' This wrapper presents the results in that exact sequence. +#' +#' @param formula A formula of the form `response ~ group`. +#' @param data A data frame containing the variables in the formula. +#' @param test The trend test to use (e.g., "jonckheereTest"). Passed to the original function. +#' @param ... Other arguments to be passed to PMCMRplus::stepDownTrendTest (e.g., alternative). +#' @return An object of class `trendPMCMR`, with p-values and statistics ordered +#' from the full dataset down to the smallest comparison (Control vs. T1). +#' @export +stepDownTrendTest_NOEC <- function(formula, data, test = "jonckheereTest", ...) { + + # --- Step 1: Call the original PMCMRplus function --- + # This gets us all the necessary calculations, even if they are in the "wrong" order. + original_result <- PMCMRplus::stepDownTrendTest( + formula = formula, + data = data, + test = test, + ... + ) + + # --- Step 2: Reverse the order of the results --- + # The original function outputs results in a "step-up" manner. We simply + # reverse the rows of the p-value and statistic matrices to get a "step-down" order. + + num_rows <- nrow(original_result$p.value) + reverse_order_indices <- num_rows:1 + + # Reverse the p-value matrix rows. `drop = FALSE` prevents R from simplifying + # the matrix to a vector if it only has one column. + p_reordered <- original_result$p.value[reverse_order_indices, , drop = FALSE] + + # Reverse the statistic matrix rows. + stat_reordered <- original_result$statistic[reverse_order_indices, , drop = FALSE] + + # --- Step 3: Create meaningful row names for the new order --- + # The original row names (e.g., "2", "3", "4") indicate the highest group included. + # We can create more descriptive names for our step-down procedure. + group_levels <- levels(data[[all.vars(formula)[2]]]) + new_rownames <- sapply(num_rows:1, function(i) { + # i=num_rows is the full test, i=1 is the control vs lowest dose test. + highest_group_name <- group_levels[i + 1] + paste0("Control through ", highest_group_name) + }) + + rownames(p_reordered) <- new_rownames + rownames(stat_reordered) <- new_rownames + + # --- Step 4: Build and return the new results object --- + # We copy the original object and just replace the parts we've changed. + new_result <- original_result + new_result$p.value <- p_reordered + new_result$statistic <- stat_reordered + + # Let's also update the method description to avoid confusion + new_result$method <- paste("Re-ordered Step-down", original_result$method) + + return(new_result) +} diff --git a/_pkgdown.yml b/_pkgdown.yml index d999a71..ebbc234 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -164,6 +164,7 @@ articles: - articles/Normality-Check - articles/Binomial_Extra_Variance - articles/Equivalence-Testing + - articles/A-Notes-on-Power - title: "Alternative Methods" desc: > diff --git a/data/jdata.rda b/data/jdata.rda new file mode 100644 index 0000000000000000000000000000000000000000..28128ec11d85166864b8e8cfc17e5ca84b80d4b8 GIT binary patch literal 249 zcmV#XNJ8l|HD^=5FS^#I)r140lHJ_|j?uKa05{oG3_Q-9T;t@4RTU literal 0 HcmV?d00001 diff --git a/man/jdata.Rd b/man/jdata.Rd new file mode 100644 index 0000000..655ae7e --- /dev/null +++ b/man/jdata.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data_description.R +\name{jdata} +\alias{jdata} +\title{Jonckheere's Synthetic Dataset} +\format{ +A data frame with 16 rows and 2 columns: +\describe{ +\item{X}{Number of samples, indicating the group each observation belongs to.} +\item{Y}{Value of the observation.} +} +} +\source{ +Jonckheere, A. R. (1954). A Distribution-Free k-Sample Test Against Ordered Alternatives. \emph{Biometrika}, 41, 133-145. +} +\usage{ +data(jdata) +} +\description{ +This dataset is used by Jonckheere (1954) to illustrate a distribution-free k-sample test against ordered alternatives. +} +\examples{ +data(jdata) +str(jdata) +summary(jdata) + +} diff --git a/man/noec_from_trend_test.Rd b/man/noec_from_trend_test.Rd new file mode 100644 index 0000000..e1475c1 --- /dev/null +++ b/man/noec_from_trend_test.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stepDownTrendTest_wrapper.R +\name{noec_from_trend_test} +\alias{noec_from_trend_test} +\title{Determine NOEC/LOEC using a Step-down Trend Test} +\usage{ +noec_from_trend_test(formula, data, test = "jonckheereTest", alpha = 0.05, ...) +} +\arguments{ +\item{formula}{A formula of the form \code{response ~ group}.} + +\item{data}{A data frame containing the variables. The group variable should be a factor.} + +\item{test}{The trend test to use (passed to PMCMRplus). Defaults to "jonckheereTest".} + +\item{alpha}{The significance level to use for determining the NOEC. Defaults to 0.05.} + +\item{...}{Other arguments passed to \code{PMCMRplus::stepDownTrendTest}.} +} +\value{ +A list of class \code{noecTrendTest} containing: +\item{NOEC}{The determined No-Observed-Effect Concentration.} +\item{LOEC}{The determined Lowest-Observed-Effect Concentration.} +\item{alpha}{The significance level used.} +\item{full_results}{The complete, original \code{trendPMCMR} object for inspection.} +} +\description{ +This function acts as a wrapper around PMCMRplus::stepDownTrendTest to +directly calculate the No-Observed-Effect Concentration (NOEC) and +Lowest-Observed-Effect Concentration (LOEC). +} +\details{ +It interprets the "step-up" ordered results from the underlying PMCMRplus +function by searching backwards from the full test to find the first +non-significant result (p >= alpha), as is standard for NOEC determination. +} diff --git a/man/print.noecTrendTest.Rd b/man/print.noecTrendTest.Rd new file mode 100644 index 0000000..d93db52 --- /dev/null +++ b/man/print.noecTrendTest.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stepDownTrendTest_wrapper.R +\name{print.noecTrendTest} +\alias{print.noecTrendTest} +\title{Print method for noecTrendTest objects} +\usage{ +\method{print}{noecTrendTest}(x, ...) +} +\arguments{ +\item{x}{An object of class \code{noecTrendTest}.} + +\item{...}{Other arguments (not used).} +} +\description{ +Print method for noecTrendTest objects +} diff --git a/man/report_dunnett_summary.Rd b/man/report_dunnett_summary.Rd index c5141ca..020b4ba 100644 --- a/man/report_dunnett_summary.Rd +++ b/man/report_dunnett_summary.Rd @@ -49,9 +49,6 @@ into a single, formatted output table. } \examples{ \dontrun{ -# Assuming `dunnett_test` and helper functions are defined in your environment -# and the 'multcomp' package is loaded. - # Generate sample data set.seed(42) my_data <- data.frame( diff --git a/man/stepDownTrendTest_NOEC.Rd b/man/stepDownTrendTest_NOEC.Rd new file mode 100644 index 0000000..fcc5761 --- /dev/null +++ b/man/stepDownTrendTest_NOEC.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stepDownTrendTest_wrapper.R +\name{stepDownTrendTest_NOEC} +\alias{stepDownTrendTest_NOEC} +\title{A Wrapper for PMCMRplus::stepDownTrendTest for NOEC Determination} +\usage{ +stepDownTrendTest_NOEC(formula, data, test = "jonckheereTest", ...) +} +\arguments{ +\item{formula}{A formula of the form \code{response ~ group}.} + +\item{data}{A data frame containing the variables in the formula.} + +\item{test}{The trend test to use (e.g., "jonckheereTest"). Passed to the original function.} + +\item{...}{Other arguments to be passed to PMCMRplus::stepDownTrendTest (e.g., alternative).} +} +\value{ +An object of class \code{trendPMCMR}, with p-values and statistics ordered +from the full dataset down to the smallest comparison (Control vs. T1). +} +\description{ +This function calls PMCMRplus::stepDownTrendTest but reorders the results +to follow a true "step-down" procedure, making it intuitive for finding a +NOEC (No-Observed-Effect Concentration). +} +\details{ +The standard ecotoxicology step-down procedure is: +\enumerate{ +\item Test for a trend across all treatment groups vs. control. +\item If significant, remove the highest dose and re-test. +\item Repeat until the test is not significant. The highest dose in that +non-significant test is the NOEC. +} + +This wrapper presents the results in that exact sequence. +} diff --git a/vignettes/articles/Trend-Testing.Rmd b/vignettes/articles/Trend-Testing.Rmd index b89cb34..4be2156 100644 --- a/vignettes/articles/Trend-Testing.Rmd +++ b/vignettes/articles/Trend-Testing.Rmd @@ -1,5 +1,7 @@ --- title: "Trend Testing" +editor_options: + chunk_output_type: console --- ```{r, include = FALSE} @@ -61,7 +63,65 @@ A step-down trend test is a sequential testing procedure used to identify the lo This approach helps identify the No Observed Effect Concentration (NOEC) or Lowest Observed Effect Concentration (LOEC) in dose-response studies. -The function `PMCMRplus::stepDownTrendTest` PMCMRplus implements this strategy for continuous data using various trend tests. For binomial survival data, we adapt this approach using the Cochran-Armitage test, which is specifically designed for binary outcomes across ordered groups, the function for performing step-down CA is called `drcHelper::stepDownTrendTestBinom`. There is also `drcHelper::step_down_RSCABS` and `drcHelper::stepDownRSCABS` for oridinal data. Please refere to the relevant data page for usage demo and explanations. +The function `PMCMRplus::stepDownTrendTest` PMCMRplus implements this strategy for continuous data using various trend tests. However, the stepDownTrendTest function does not follow the NOEC testing paradigm of starting with a full model and removing the highest dose. Instead, it builds the model upwards from the control. For binomial survival data, we adapt this approach using the Cochran-Armitage test, which is specifically designed for binary outcomes across ordered groups, the function for performing step-down CA is called `drcHelper::stepDownTrendTestBinom`. There is also `drcHelper::step_down_RSCABS` and `drcHelper::stepDownRSCABS` for ordinal data. Please refere to the relevant data page for usage demo and explanations. +```{r echo=FALSE} +data(jdata) +jdata$X <- factor(jdata$X) + +# --- 1. The Official PMCMRplus::stepDownTrendTest --- +# This is our reference. +official_result <- PMCMRplus::stepDownTrendTest(Y ~ X, data = jdata, test = "jonckheereTest") +cat("--- PMCMRplus::stepDownTrendTest Output ---\n") +print(official_result) +``` +```{r} +# --- 2. Manually Recreating the PMCMRplus "Step-Up" Logic --- +# This manually implements the logic inside the function to prove our understanding. +cat("\n--- Manual Recreation of PMCMRplus 'Step-Up' Logic ---\n") +levels_X <- levels(jdata$X) +# Loop from i=2 to 4 +for (i in 2:length(levels_X)) { + # Create a subset of data with groups from 1 to i + subset_data <- subset(jdata, X %in% levels_X[1:i]) + + # Run the test on this subset + res <- PMCMRplus::jonckheereTest(Y ~ X, data = subset_data) + + cat(sprintf("Test on Groups %s vs 1: p-value = %.6f\n", i, res$p.value)) +} + + +# --- 3. The Ecotoxicology "Step-Down" Logic (Your Manual Approach) --- +# This starts with all doses and removes the highest. +cat("\n--- Manual Recreation of Ecotoxicology 'Step-Down' Logic ---\n") +# Loop from i=4 down to 2 +for (i in length(levels_X):2) { + # Create a subset of data with groups from 1 to i + subset_data <- subset(jdata, X %in% levels_X[1:i]) + + # Run the test on this subset + res <- PMCMRplus::jonckheereTest(Y ~ X, data = subset_data) + + cat(sprintf("Test on Groups 1 to %s: p-value = %.6f\n", i, res$p.value)) +} +``` +Therefore, in the `drcHelper`, a wrapper function is implemented to give the results in correct meaningful order. + + +```{r} + +# --- Run our NEW WRAPPER function --- +cat("\n=========================================\n") +cat(" Output from our NEW WRAPPER function\n") +cat("=========================================\n") +wrapper_output <- stepDownTrendTest_NOEC(Y ~ X, data = jdata, test = "jonckheereTest") +print(wrapper_output) +``` + +```{r} +noec_output <- noec_from_trend_test(Response ~ Dose, data = dat_medium, test = "jonckheereTest") +print(noec_output) +``` diff --git a/vignettes/articles/Verification_which_JT.Rmd b/vignettes/articles/Verification_which_JT.Rmd index 0b95a62..c3f2de2 100644 --- a/vignettes/articles/Verification_which_JT.Rmd +++ b/vignettes/articles/Verification_which_JT.Rmd @@ -7,22 +7,24 @@ editor_options: ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, - comment = "#>" + comment = "#>", + eval = TRUE, + warning=FALSE, + message = FALSE ) -source("../knitr-setup.R") ``` -```r +```{r setup} library(drcHelper) ``` - -```r +```{r} library(tidyverse) ``` + ## Summary Jonckheere's test and Kendall's tau) are closely related. In some software implementations, the two test produce identical p-values. @@ -30,136 +32,265 @@ Jonckheere's test and Kendall's tau) are closely related. In some software imple `JtTest` from `npordtests` are doing the same as in `cor.test` with method being "kendall". Both function actually performs a kendall test and treat the dose column as numeric. On the other hand `PMCMRplus::jonckheereTest` and `DescTools::JonckheereTerpstraTest` do similar things. -## Data from Jonckheere (1954) +The differences in different implementations almost always come down to three things: +1. **How ties are handled** in the variance calculation. +2. **Whether a continuity correction** is applied to the z-score. +3. **Whether the p-value is derived from a normal approximation or an exact permutation test**. -```r -data(jdata) +Packages like `DescTools` and `PMCMRplus` are implementing the Jonckheere-Terpstra test directly, while `npordtests::JtTest` and `cor.test(method="kendall")` are leveraging the very close mathematical relationship between the JT statistic and Kendall's S statistic. + +### The Core Calculations + +Here is the standard procedure for the Jonckheere-Terpstra test using the normal approximation, which is what these packages are using for datasets of this size. + +Let's use your `jdata` example as the gold standard to walk through. +* k = 4 groups (doses) +* Group sizes: n1=4, n2=4, n3=4, n4=4 +* Total observations N = 16 + +#### Step 1: Calculate the Jonckheere-Terpstra (JT) Statistic + +The JT statistic is the sum of all Mann-Whitney U statistics for every pair of groups (i, j) where **group i comes before group j in the hypothesized order**. + +$$JT = \sum_{i=1}^{k-1} \sum_{j=i+1}^{k} U_{ij}$$ + +Where $U_{ij}$ is the number of pairs of observations (x, y) where x is from group i, y is from group j, and x < y. + +Using `jdata`, all packages agree on the statistic: **JT = 71**. + +#### Step 2: Calculate the Mean and Variance of JT under the Null Hypothesis + +This is where implementations can start to diverge. + +**Mean (E[JT]):** The formula is simple and universal. +$\mu_{JT} = \frac{N^2 - \sum_{i=1}^{k} n_i^2}{4}$ + +For `jdata`: +* N = 16 +* $$\sum n_i^2 = 4^2 + 4^2 + 4^2 + 4^2 = 64$$ +* $$\mu_{JT} = \frac{16^2 - 64}{4} = \frac{256 - 64}{4} = \frac{192}{4} = 48$$ + +This matches the mean `48` from `npordtests::JtTest`. Your implementation should also get this value. + +**Variance (Var(JT)):** The formula gets more complex, especially with ties. + +* **No Ties:** + $$\sigma^2_{JT} = \frac{N^2(2N+3) - \sum_{i=1}^{k} n_i^2(2n_i+3)}{72}$$ + +* **With Ties (The Correct Formula to Use):** + $$\sigma^2_{JT} = \frac{A}{72} - \frac{B}{12N(N-1)}$$ + Where: + $$A = N(N-1)(2N+5) - \sum n_i(n_i-1)(2n_i+5)$$ + $$B = [\sum n_i(n_i-1)(n_i-2)] \times [\sum t_j(t_j-1)(t_j-2)] + [\sum n_i(n_i-1)] \times [\sum t_j(t_j-1)]$$ + (t_j is the number of observations in a tied group of responses) + +Let's analyze the ties in `jdata$Y`: `{10, 20, 35, 38, 41, 43, 51, 55, 62, 70, 78, 85}`. There are no ties in the response variable `Y`. This simplifies the variance calculation greatly. The general formula for variance with ties in the *grouping* variable but not the *response* is: +$$\sigma^2_{JT} = \frac{N^2(2N+3) - \sum n_i^2(2n_i+3)}{72}$$ +(This is the same as the "no ties" formula, as "ties" usually refers to the response values). + +For `jdata`: +* N=16, n_i=4 +* $$\sum n_i^2(2n_i+3) = 4 \times [4^2(2 \cdot 4 + 3)] = 4 \times [16(11)] = 4 \times 176 = 704$$ +* $$N^2(2N+3) = 16^2(2 \cdot 16 + 3) = 256(35) = 8960$$ +* $$\sigma^2_{JT} = \frac{8960 - 704}{72} = \frac{8256}{72} = 114.6667$$ -prelimPlot3(jdata,dose_col = "X",response_col = "Y") +This **exactly matches** the variance from `npordtests::JtTest`. + +#### Step 3: Calculate the Z-statistic + +This is the second place where implementations differ. `npordtests` and `cor.test` use the Z-score calculated directly from the mean and variance above, without a continuity correction. `DescTools` and `PMCMRplus` are using a slightly different calculation invovling a continuity correction + +**Standard Z-statistic:** +$$Z = \frac{JT - \mu_{JT}}{\sigma_{JT}}$$ + +For `jdata`: +$$Z = \frac{71 - 48}{\sqrt{114.6667}} = \frac{23}{10.70825} = 2.147876$$ +This **exactly matches** `npordtests` and `cor.test`! + +```{r} +p_value <- 2 * pnorm(2.147876, lower.tail = FALSE) +p_value ``` -![plot of chunk unnamed-chunk-3](figure/unnamed-chunk-3-1.png) - -```r - -dunnett_test(jdata,"Y","X",include_random_effect = FALSE) -#> Dunnett Test Results -#> ------------------- -#> Model type: Fixed model with homoscedastic errors -#> Control level: 1 -#> Alpha level: 0.05 -#> -#> Results Table: -#> comparison estimate std.error statistic p.value conf.low conf.high significant -#> 2 - 1 15.50 34.05051 0.4552061 0.9401315 -75.78264 106.7826 FALSE -#> 3 - 1 39.75 34.05051 1.1673833 0.5307511 -51.53264 131.0326 FALSE -#> 4 - 1 60.25 34.05051 1.7694300 0.2322884 -31.03264 151.5326 FALSE -#> -#> NOEC Determination: -#> No significant effects detected at any dose. NOEC is at or above the highest tested dose. + +**Z-statistic with Continuity Correction:** +The correction is used for discrete distributions approximated by a continuous one. We subtract 0.5 from the numerator. +$$Z_{corr} = \frac{JT - \mu_{JT} - 0.5}{\sigma_{JT}}$$ + +For `jdata`: +$$Z_{corr} = \frac{71 - 48 - 0.5}{\sqrt{114.6667}} = \frac{22.5}{10.70825} = 2.10115$$ +```{r} +p_value <- 2 * pnorm(2.10115, lower.tail = FALSE) +p_value ``` +#### Step 4: Reconciling the p-values + +Now we can see why the p-values differ. + +* **`npordtests::JtTest` and `cor.test`:** They use the Z-score **without** continuity correction (Z = 2.147876). They report a **two-sided** p-value. + +* **`DescTools::JonckheereTerpstraTest`:** This package often uses continuity correction by default for its non-parametric tests. Let's check the p-value with our corrected Z-score (Z_corr = 2.10115). +* **`PMCMRplus::jonckheereTest`:** This package also reports a two-sided p-value (`p-value = 0.001726` for `lehmann` data). For the `jdata` set, the difference you noted (`-0.00196`) again points to a slightly different calculation, likely involving continuity correction. For the `lehmann` data, its Z-score of `3.1337` is slightly different from `npordtests`'s `3.1254`, indicating a different variance calculation (due to the many ties in the `lehmann` data) and/or a continuity correction. -```r + +### Compare the 4 Packages + +```{r eval=FALSE} +# Load necessary libraries library(npordtests) -res1 <- JtTest(Y~X,jdata) -#> --------------------------------------------------------- -#> Test : Jonckheere-Terpstra Test -#> data : Y and X -#> -#> Statistic = 71 -#> Mean = 48 -#> Variance = 114.6667 -#> Z = 2.147876 -#> Asymp. p-value = 0.0158618 -#> -#> Result : Null hypothesis is rejected. -#> --------------------------------------------------------- -res1$p.value *2 -#> [1] 0.03172359 -``` +library(DescTools) +library(PMCMRplus) +library(drcHelper) # For the jdata dataset +# --- Load Datasets --- -```r -res2 <- DescTools::JonckheereTerpstraTest(Y~X,jdata) -res2 -#> -#> Jonckheere-Terpstra test -#> -#> data: Y by X -#> JT = 71, p-value = 0.03368 -#> alternative hypothesis: two.sided -res2$statistic == res1$statistic -#> JT -#> TRUE -``` +# Data from Jonckheere (1954), provided in drcHelper +# The user confirmed these are the correct Y values +# Group (X): 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 +# Y : 19, 20, 60, 130, 21, 61, 80, 129, 40, 99, 100, 149, 49, 110, 151, 160 +data(jdata) +# Data from Lehmann (1975), provided in multiple packages +data(lehmann) -```r -res3 <- PMCMRplus::jonckheereTest(Y~X,jdata) -res3$statistic=res2$statistic -res3$p.value - res2$p.value -#> [1] -0.001960184 +# --- Reusable Summarization Function --- + +#' Summarize results from three different Jonckheere-Terpstra test implementations +#' +#' @param data The data frame to use. +#' @param formula A formula of the form response ~ group. +#' @return A data frame summarizing the results. +summarize_jt_tests <- function(data, formula) { + + results_list <- list() + + # 1. npordtests::JtTest + # Note: This returns a one-sided p-value by default. We double it for comparison. + # This implementation is mathematically equivalent to cor.test(..., method="kendall") + res1 <- npordtests::JtTest(formula, data = data) + results_list[[1]] <- data.frame( + Package = "npordtests", + Function = "JtTest", + JT_Statistic = as.numeric(res1$statistic), + Z_Statistic = NA, + P_Value = res1$p.value * 2, # Standardizing to two-sided + Alternative = "two.sided", + Notes = "P-value is doubled from one-sided output." + ) + + # 2. DescTools::JonckheereTerpstraTest + # Note: This implementation often includes a continuity correction. + res2 <- DescTools::JonckheereTerpstraTest(formula, data = data, alternative = "two.sided") + results_list[[2]] <- data.frame( + Package = "DescTools", + Function = "JonckheereTerpstraTest", + JT_Statistic = res2$statistic, + Z_Statistic = NA, # Not directly returned in the htest object + P_Value = res2$p.value, + Alternative = res2$alternative, + Notes = "Uses continuity correction. Z-stat not returned." + ) + + # 3. PMCMRplus::jonckheereTest + # Note: This also uses a continuity correction. + # The JT statistic is in 'estimate', the Z-statistic is in 'statistic'. + res3 <- PMCMRplus::jonckheereTest(formula, data = data, alternative = "two.sided") + results_list[[3]] <- data.frame( + Package = "PMCMRplus", + Function = "jonckheereTest", + JT_Statistic = res3$estimate, + Z_Statistic = res3$statistic, + P_Value = res3$p.value, + Alternative = res3$alternative, + Notes = "JT is in 'estimate'; Z is in 'statistic'." + ) + # For bonus validation: a quick check with cor.test + group_var <- all.vars(formula)[2] + response_var <- all.vars(formula)[1] + res_kendall <- cor.test(as.numeric(data[[group_var]]), as.numeric(data[[response_var]]), method = "kendall") + + results_list[[4]] <- data.frame( + Package = "stats", + Function = "cor.test", + JT_Statistic = NA, + Z_Statistic = res3$statistic, + P_Value = res3$p.value, + Alternative = res3$alternative, + Notes = "Z is in 'statistic'." + ) + # Combine all results into a single data frame + final_df <- do.call(rbind, results_list) + rownames(final_df) <- NULL # Clean up row names + + + return(final_df) +} + + +# --- Run and Print Summaries --- + +cat("==========================================\n") +cat(" Summary for 'jdata' Dataset\n") +cat("==========================================\n") +jdata_summary <- summarize_jt_tests(jdata, Y ~ X) +## print(jdata_summary, row.names = FALSE) +jdata_summary%>%knitr::kable(.,digits = 4) + +cat("\n\n==========================================\n") +cat(" Summary for 'lehmann' Dataset\n") +cat("==========================================\n") +lehmann_summary <- summarize_jt_tests(lehmann, Values ~ Group) +## print(lehmann_summary, row.names = FALSE) +lehmann_summary%>%knitr::kable(.,digits = 4) ``` +========================================== + Summary for 'jdata' Dataset +========================================== +|Package |Function | JT_Statistic| Z_Statistic| P_Value|Alternative |Notes | +|:----------|:----------------------|------------:|-----------:|-------:|:-----------|:------------------------------------------------| +|npordtests |JtTest | 71| | 0.0317|two.sided |P-value is doubled from one-sided output. | +|DescTools |JonckheereTerpstraTest | 71| | 0.0337|two.sided |Uses continuity correction. Z-stat not returned. | +|PMCMRplus |jonckheereTest | 71| 2.1479| 0.0317|two.sided |JT is in 'estimate'; Z is in 'statistic'. | +|stats |cor.test | | 2.1479| 0.0317|two.sided |Z is in 'statistic'. -```r -res4 <- cor.test(as.numeric(jdata$X),as.numeric(jdata$Y),method="kendall") -res4$statistic -#> z -#> 2.147876 -res4$p.value -#> [1] 0.03172359 -``` +========================================== + Summary for 'lehmann' Dataset +========================================== +|Package |Function | JT_Statistic| Z_Statistic| P_Value|Alternative |Notes | +|:----------|:----------------------|------------:|-----------:|-------:|:-----------|:------------------------------------------------| +|npordtests |JtTest | 1159| | 0.0018|two.sided |P-value is doubled from one-sided output. | +|DescTools |JonckheereTerpstraTest | 1159| | 0.0018|two.sided |Uses continuity correction. Z-stat not returned. | +|PMCMRplus |jonckheereTest | 1159| 3.1337| 0.0017|two.sided |JT is in 'estimate'; Z is in 'statistic'. | +|stats |cor.test | | 3.1337| 0.0017|two.sided |Z is in 'statistic'. +## Additional Difference -```r -data(lehmann) -res <- JtTest(Values~Group,lehmann) -#> --------------------------------------------------------- -#> Test : Jonckheere-Terpstra Test -#> data : Values and Group -#> -#> Statistic = 1159 -#> Mean = 857.5 -#> Variance = 9305.917 -#> Z = 3.125415 -#> Asymp. p-value = 0.0008877709 -#> -#> Result : Null hypothesis is rejected. -#> --------------------------------------------------------- -res4$p.value - 2*res1$p.value -#> [1] 1.387779e-17 +Another difference comes from the *dose group* column. When it is a factor or a numerical, the results are different using different testing function. + +```{r} +prelimPlot3(dat_medium) ``` -## Data from Lehmann (1975) - - -```r -data("lehmann") -DescTools::JonckheereTerpstraTest(Values~Group,lehmann) -#> -#> Jonckheere-Terpstra test -#> -#> data: Values by Group -#> JT = 1159, p-value = 0.001776 -#> alternative hypothesis: two.sided -PMCMRplus::jonckheereTest(Values~Group,lehmann) -#> -#> Jonckheere-Terpstra test -#> -#> data: Values by Group -#> z = 3.1337, p-value = 0.001726 -#> alternative hypothesis: two.sided -#> sample estimates: -#> JT -#> 1159 +```{r eval=FALSE} +dat1 <- dat_medium %>% mutate(Dose = as.numeric(as.character(Dose))) %>% filter(Dose <7) +dat_medium_summary <- summarize_jt_tests(dat1 %>% mutate(Dose = factor(Dose))%>% mutate(Response = -Response), Response ~ Dose) +dat_medium_summary %>%knitr::kable(.,digits = 4) +JtTest(Response ~ Dose,data = dat1 %>% mutate(Dose = factor(Dose)) %>% mutate(Response = -Response)) ``` +|Package |Function | JT_Statistic| Z_Statistic| P_Value|Alternative |Notes | +|:----------|:----------------------|------------:|-----------:|-------:|:-----------|:------------------------------------------------| +|npordtests |JtTest | 72| | 1e-03|two.sided |P-value is doubled from one-sided output. | +|DescTools |JonckheereTerpstraTest | 72| | 5e-04|two.sided |Uses continuity correction. Z-stat not returned. | +|PMCMRplus |jonckheereTest | 72| 3.2796| 1e-03|two.sided |JT is in 'estimate'; Z is in 'statistic'. | +|stats |cor.test | | 3.2796| 1e-03|two.sided |Z is in 'statistic'. + diff --git a/vignettes/articles/Verification_which_JT.Rmd.orig b/vignettes/articles/Verification_which_JT.Rmd.orig deleted file mode 100644 index ba663f5..0000000 --- a/vignettes/articles/Verification_which_JT.Rmd.orig +++ /dev/null @@ -1,88 +0,0 @@ ---- -title: "Verification: Which Jonckeere Terpstra Test to Use" -editor_options: - chunk_output_type: console ---- - -```{r, include = FALSE} -knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>", - eval = TRUE, - warning=FALSE, - message = FALSE -) -``` - -```{r setup} -library(drcHelper) -``` - - -```{r} - -library(tidyverse) -``` - -## Summary - -Jonckheere's test and Kendall's tau) are closely related. In some software implementations, the two test produce identical p-values. - -`JtTest` from `npordtests` are doing the same as in `cor.test` with method being "kendall". Both function actually performs a kendall test and treat the dose column as numeric. On the other hand `PMCMRplus::jonckheereTest` and `DescTools::JonckheereTerpstraTest` do similar things. - - -## Data from Jonckheere (1954) - -```{r} -data(jdata) - -prelimPlot3(jdata,dose_col = "X",response_col = "Y") - -dunnett_test(jdata,"Y","X",include_random_effect = FALSE) -``` - - -```{r} -library(npordtests) -res1 <- JtTest(Y~X,jdata) -res1$p.value *2 -``` - -```{r} -res2 <- DescTools::JonckheereTerpstraTest(Y~X,jdata) -res2 -res2$statistic == res1$statistic -``` - - -```{r} -res3 <- PMCMRplus::jonckheereTest(Y~X,jdata) -res3$statistic=res2$statistic -res3$p.value - res2$p.value - -``` - - -```{r} -res4 <- cor.test(as.numeric(jdata$X),as.numeric(jdata$Y),method="kendall") -res4$statistic -res4$p.value -``` - - - -```{r} -data(lehmann) -res <- JtTest(Values~Group,lehmann) -res4$p.value - 2*res1$p.value -``` - - -## Data from Lehmann (1975) - -```{r} -data("lehmann") -DescTools::JonckheereTerpstraTest(Values~Group,lehmann) -PMCMRplus::jonckheereTest(Values~Group,lehmann) -``` - diff --git a/vignettes/articles/Verification_which_JT.md b/vignettes/articles/Verification_which_JT.md new file mode 100644 index 0000000..0b95a62 --- /dev/null +++ b/vignettes/articles/Verification_which_JT.md @@ -0,0 +1,165 @@ +--- +title: "Verification: Which Jonckeere Terpstra Test to Use" +editor_options: + chunk_output_type: console +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +source("../knitr-setup.R") +``` + +```r +library(drcHelper) +``` + + + +```r + +library(tidyverse) +``` + +## Summary + +Jonckheere's test and Kendall's tau) are closely related. In some software implementations, the two test produce identical p-values. + +`JtTest` from `npordtests` are doing the same as in `cor.test` with method being "kendall". Both function actually performs a kendall test and treat the dose column as numeric. On the other hand `PMCMRplus::jonckheereTest` and `DescTools::JonckheereTerpstraTest` do similar things. + + +## Data from Jonckheere (1954) + + +```r +data(jdata) + +prelimPlot3(jdata,dose_col = "X",response_col = "Y") +``` + +![plot of chunk unnamed-chunk-3](figure/unnamed-chunk-3-1.png) + +```r + +dunnett_test(jdata,"Y","X",include_random_effect = FALSE) +#> Dunnett Test Results +#> ------------------- +#> Model type: Fixed model with homoscedastic errors +#> Control level: 1 +#> Alpha level: 0.05 +#> +#> Results Table: +#> comparison estimate std.error statistic p.value conf.low conf.high significant +#> 2 - 1 15.50 34.05051 0.4552061 0.9401315 -75.78264 106.7826 FALSE +#> 3 - 1 39.75 34.05051 1.1673833 0.5307511 -51.53264 131.0326 FALSE +#> 4 - 1 60.25 34.05051 1.7694300 0.2322884 -31.03264 151.5326 FALSE +#> +#> NOEC Determination: +#> No significant effects detected at any dose. NOEC is at or above the highest tested dose. +``` + + + +```r +library(npordtests) +res1 <- JtTest(Y~X,jdata) +#> --------------------------------------------------------- +#> Test : Jonckheere-Terpstra Test +#> data : Y and X +#> +#> Statistic = 71 +#> Mean = 48 +#> Variance = 114.6667 +#> Z = 2.147876 +#> Asymp. p-value = 0.0158618 +#> +#> Result : Null hypothesis is rejected. +#> --------------------------------------------------------- +res1$p.value *2 +#> [1] 0.03172359 +``` + + +```r +res2 <- DescTools::JonckheereTerpstraTest(Y~X,jdata) +res2 +#> +#> Jonckheere-Terpstra test +#> +#> data: Y by X +#> JT = 71, p-value = 0.03368 +#> alternative hypothesis: two.sided +res2$statistic == res1$statistic +#> JT +#> TRUE +``` + + + +```r +res3 <- PMCMRplus::jonckheereTest(Y~X,jdata) +res3$statistic=res2$statistic +res3$p.value - res2$p.value +#> [1] -0.001960184 +``` + + + +```r +res4 <- cor.test(as.numeric(jdata$X),as.numeric(jdata$Y),method="kendall") +res4$statistic +#> z +#> 2.147876 +res4$p.value +#> [1] 0.03172359 +``` + + + + +```r +data(lehmann) +res <- JtTest(Values~Group,lehmann) +#> --------------------------------------------------------- +#> Test : Jonckheere-Terpstra Test +#> data : Values and Group +#> +#> Statistic = 1159 +#> Mean = 857.5 +#> Variance = 9305.917 +#> Z = 3.125415 +#> Asymp. p-value = 0.0008877709 +#> +#> Result : Null hypothesis is rejected. +#> --------------------------------------------------------- +res4$p.value - 2*res1$p.value +#> [1] 1.387779e-17 +``` + + +## Data from Lehmann (1975) + + +```r +data("lehmann") +DescTools::JonckheereTerpstraTest(Values~Group,lehmann) +#> +#> Jonckheere-Terpstra test +#> +#> data: Values by Group +#> JT = 1159, p-value = 0.001776 +#> alternative hypothesis: two.sided +PMCMRplus::jonckheereTest(Values~Group,lehmann) +#> +#> Jonckheere-Terpstra test +#> +#> data: Values by Group +#> z = 3.1337, p-value = 0.001726 +#> alternative hypothesis: two.sided +#> sample estimates: +#> JT +#> 1159 +``` + From 8ba481ada1f40737b2949b25a8537ebabb1019a2 Mon Sep 17 00:00:00 2001 From: Zhenglei Gao Date: Fri, 3 Oct 2025 16:38:35 +0200 Subject: [PATCH 5/5] Increment version number to 0.0.5 --- DESCRIPTION | 2 +- NEWS.md | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0931312..49bb542 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: drcHelper Title: Collection and verification of helper functions for dose-response analysis -Version: 0.0.4.9001 +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")), diff --git a/NEWS.md b/NEWS.md index f235026..3b5d2ed 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,5 @@ +# drcHelper 0.0.5 + # drcHelper 0.0.4.9001() * Initial CRAN submission preparation.