|
1 | 1 | ## Function to calculate NOECs using multiple methods and broomed together |
| 2 | +## Unit Testing in # FILE: tests/testthat/test_quantal_categorical.R |
2 | 3 |
|
3 | 4 |
|
4 | 5 |
|
5 | | -## Unit Testing in # FILE: tests/testthat/test_quantal_categorical.R |
6 | | - |
7 | | -#' Tarone's Z Test |
| 6 | +#' Calculate NOEC Using Many-to-one Pairwise Tests |
8 | 7 | #' |
9 | | -#' Tests the goodness of fit of the binomial distribution. |
| 8 | +#' This function calculates the No Observed Effect Concentration (NOEC) from dose response data |
| 9 | +#' using pairwise comparison tests from the rstatix package. |
10 | 10 | #' |
11 | | -#' @param M Counts |
12 | | -#' @param N Trials |
| 11 | +#' @param data A data frame containing the dose response data |
| 12 | +#' @param response The name of the response variable (unquoted) |
| 13 | +#' @param dose The name of the dose variable (unquoted) |
| 14 | +#' @param control The level of the dose variable to be used as control |
| 15 | +#' @param test The statistical test to use: "t.test" or "wilcox.test" |
| 16 | +#' @param p_adjust_method Method for p-value adjustment for multiple comparisons |
| 17 | +#' @param alternative Direction of the alternative hypothesis: "two.sided", "greater", or "less" |
| 18 | +#' @param alpha Significance level (default: 0.05) |
13 | 19 | #' |
14 | | -#' @return a \code{htest} object |
15 | | -#' |
16 | | -#' @author \href{https://stats.stackexchange.com/users/173082/reinstate-monica}{Ben O'Neill} |
17 | | -#' @references \url{https://stats.stackexchange.com/a/410376/6378} and |
18 | | -#' 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} |
19 | | -#' @importFrom stats pnorm |
| 20 | +#' @return A list containing the NOEC value and the full test results |
20 | 21 | #' @export |
21 | | -#' @examples |
22 | | -#' #Generate example data |
23 | | -#' N <- c(30, 32, 40, 28, 29, 35, 30, 34, 31, 39) |
24 | | -#' M <- c( 9, 10, 22, 15, 8, 19, 16, 19, 15, 10) |
25 | | -#' Tarone.test(N, M) |
26 | | -Tarone.test <- function(N, M) { |
27 | | - |
28 | | - #Check validity of inputs |
29 | | - if(!(all(N == as.integer(N)))) { stop("Error: Number of trials should be integers"); } |
30 | | - if(min(N) < 1) { stop("Error: Number of trials should be positive"); } |
31 | | - if(!(all(M == as.integer(M)))) { stop("Error: Count values should be integers"); } |
32 | | - if(min(M) < 0) { stop("Error: Count values cannot be negative"); } |
33 | | - if(any(M > N)) { stop("Error: Observed count value exceeds number of trials"); } |
34 | | - |
35 | | - #Set description of test and data |
36 | | - method <- "Tarone's Z test"; |
37 | | - data.name <- paste0(deparse(substitute(M)), " successes from ", |
38 | | - deparse(substitute(N)), " trials"); |
39 | | - |
40 | | - #Set null and alternative hypotheses |
41 | | - null.value <- 0; |
42 | | - attr(null.value, "names") <- "dispersion parameter"; |
43 | | - alternative <- "greater"; |
44 | | - |
45 | | - #Calculate test statistics |
46 | | - estimate <- sum(M)/sum(N); |
47 | | - attr(estimate, "names") <- "proportion parameter"; |
48 | | - S <- ifelse(estimate == 1, sum(N), |
49 | | - sum((M - N*estimate)^2/(estimate*(1 - estimate)))); |
50 | | - statistic <- (S - sum(N))/sqrt(2*sum(N*(N-1))); |
51 | | - attr(statistic, "names") <- "z"; |
52 | | - |
53 | | - #Calculate p-value |
54 | | - p.value <- 2*pnorm(-abs(statistic), 0, 1); |
55 | | - attr(p.value, "names") <- NULL; |
56 | | - |
57 | | - #Create htest object |
58 | | - TEST <- list(method = method, data.name = data.name, |
59 | | - null.value = null.value, alternative = alternative, |
60 | | - estimate = estimate, statistic = statistic, p.value = p.value); |
61 | | - class(TEST) <- "htest"; |
62 | | - TEST; |
| 22 | +#' |
| 23 | +#' @importFrom dplyr mutate filter arrange %>% |
| 24 | +#' @importFrom rlang enquo quo_name .data |
| 25 | +#' @importFrom rstatix t_test wilcox_test |
| 26 | +#' @importFrom stats as.formula setNames |
| 27 | +calculate_noec_rstatix <- function(data, response, dose, control = "0", |
| 28 | + test = c("t.test", "wilcox.test"), |
| 29 | + p_adjust_method = "holm", |
| 30 | + alternative = "two.sided", |
| 31 | + alpha = 0.05) { |
| 32 | + |
| 33 | + # Match test argument |
| 34 | + test <- match.arg(test) |
| 35 | + |
| 36 | + # Capture variable names using rlang |
| 37 | + response_var <- rlang::enquo(response) |
| 38 | + dose_var <- rlang::enquo(dose) |
| 39 | + response_name <- rlang::quo_name(response_var) |
| 40 | + dose_name <- rlang::quo_name(dose_var) |
| 41 | + |
| 42 | + # Ensure dose is a factor |
| 43 | + if (!is.factor(data[[dose_name]])) { |
| 44 | + data[[dose_name]] <- as.factor(data[[dose_name]]) |
| 45 | + } |
| 46 | + |
| 47 | + # Create formula |
| 48 | + formula_obj <- stats::as.formula(paste(response_name, "~", dose_name)) |
| 49 | + |
| 50 | + # Run the appropriate test |
| 51 | + if (test == "t.test") { |
| 52 | + test_results <- rstatix::t_test( |
| 53 | + data = data, |
| 54 | + formula = formula_obj, |
| 55 | + ref.group = control, |
| 56 | + p.adjust.method = p_adjust_method, |
| 57 | + alternative = alternative |
| 58 | + ) |
| 59 | + } else if (test == "wilcox.test") { |
| 60 | + test_results <- rstatix::wilcox_test( |
| 61 | + data = data, |
| 62 | + formula = formula_obj, |
| 63 | + ref.group = control, |
| 64 | + p.adjust.method = p_adjust_method, |
| 65 | + alternative = alternative |
| 66 | + ) |
| 67 | + } |
| 68 | + |
| 69 | + # Extract dose levels and convert to numeric for sorting |
| 70 | + dose_levels <- levels(data[[dose_name]]) |
| 71 | + numeric_doses <- suppressWarnings(as.numeric(dose_levels)) |
| 72 | + |
| 73 | + # If conversion to numeric was successful, reorder the test results |
| 74 | + if (!any(is.na(numeric_doses))) { |
| 75 | + # Create a mapping from dose level to numeric value |
| 76 | + dose_mapping <- stats::setNames(numeric_doses, dose_levels) |
| 77 | + |
| 78 | + # Add numeric dose values to test results |
| 79 | + test_results <- test_results %>% |
| 80 | + dplyr::mutate( |
| 81 | + group1_numeric = dose_mapping[.data$group1], |
| 82 | + group2_numeric = dose_mapping[.data$group2] |
| 83 | + ) |
| 84 | + |
| 85 | + # Determine which column contains the non-control doses |
| 86 | + if (all(test_results$group1 == control)) { |
| 87 | + test_results <- test_results %>% |
| 88 | + dplyr::mutate(dose_numeric = .data$group2_numeric) |
| 89 | + } else if (all(test_results$group2 == control)) { |
| 90 | + test_results <- test_results %>% |
| 91 | + dplyr::mutate(dose_numeric = .data$group1_numeric) |
| 92 | + } else { |
| 93 | + # Mixed case - need to handle both possibilities |
| 94 | + test_results <- test_results %>% |
| 95 | + dplyr::mutate( |
| 96 | + dose_numeric = ifelse(.data$group1 == control, |
| 97 | + .data$group2_numeric, |
| 98 | + .data$group1_numeric) |
| 99 | + ) |
| 100 | + } |
| 101 | + |
| 102 | + # Sort by numeric dose |
| 103 | + test_results <- test_results %>% |
| 104 | + dplyr::arrange(.data$dose_numeric) |
| 105 | + } |
| 106 | + |
| 107 | + # Find the NOEC (highest dose with p > alpha) |
| 108 | + significant_results <- test_results %>% |
| 109 | + dplyr::filter(.data$p.adj <= alpha) |
| 110 | + |
| 111 | + if (nrow(significant_results) == 0) { |
| 112 | + # If no significant differences, NOEC is the highest tested dose |
| 113 | + all_doses <- as.numeric(as.character(unique(data[[dose_name]]))) |
| 114 | + all_doses <- all_doses[!is.na(all_doses) & all_doses != as.numeric(control)] |
| 115 | + |
| 116 | + if (length(all_doses) > 0) { |
| 117 | + noec_value <- max(all_doses, na.rm = TRUE) |
| 118 | + noec_message <- "No significant effects detected at any dose level" |
| 119 | + } else { |
| 120 | + noec_value <- NA |
| 121 | + noec_message <- "No valid dose levels found" |
| 122 | + } |
| 123 | + } else { |
| 124 | + # Get all doses with significant effects |
| 125 | + significant_doses <- c() |
| 126 | + |
| 127 | + for (i in 1:nrow(significant_results)) { |
| 128 | + row <- significant_results[i, ] |
| 129 | + if (row$group1 == control) { |
| 130 | + significant_doses <- c(significant_doses, as.character(row$group2)) |
| 131 | + } else { |
| 132 | + significant_doses <- c(significant_doses, as.character(row$group1)) |
| 133 | + } |
| 134 | + } |
| 135 | + |
| 136 | + # Convert to numeric for comparison |
| 137 | + significant_doses <- suppressWarnings(as.numeric(significant_doses)) |
| 138 | + |
| 139 | + # Get all tested doses |
| 140 | + all_doses <- suppressWarnings(as.numeric(as.character(unique(data[[dose_name]])))) |
| 141 | + all_doses <- all_doses[!is.na(all_doses) & all_doses != as.numeric(control)] |
| 142 | + all_doses <- sort(all_doses) |
| 143 | + |
| 144 | + # Find the highest non-significant dose |
| 145 | + non_significant_doses <- all_doses[!all_doses %in% significant_doses] |
| 146 | + |
| 147 | + if (length(non_significant_doses) == 0) { |
| 148 | + # If all doses show significant effects |
| 149 | + noec_value <- as.numeric(control) |
| 150 | + noec_message <- "All tested doses showed significant effects" |
| 151 | + } else { |
| 152 | + noec_value <- max(non_significant_doses) |
| 153 | + noec_message <- paste("NOEC determined as", noec_value) |
| 154 | + } |
| 155 | + } |
| 156 | + |
| 157 | + # Return both the NOEC value and the full test results |
| 158 | + return(list( |
| 159 | + noec = noec_value, |
| 160 | + noec_message = noec_message, |
| 161 | + test_results = test_results |
| 162 | + )) |
63 | 163 | } |
64 | 164 |
|
65 | 165 |
|
66 | 166 |
|
67 | 167 |
|
| 168 | + |
| 169 | + |
| 170 | + |
| 171 | + |
| 172 | + |
| 173 | + |
| 174 | + |
| 175 | + |
| 176 | + |
0 commit comments