diff --git a/_pkgdown.yml b/_pkgdown.yml index 9f8256b..7e41a44 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -55,6 +55,7 @@ navbar: - text: Example ECx Helper Functions Usage href: articles/drcHelper.html - text: Example NOEC Helper Functions Usage + href: articles/Dunnetts_Test_for_Data_with_Hierarchical_Structure.html - text: Using DRC and BMD href: articles/Examples_drc.html - text: Example DRC - OECD 201 @@ -107,6 +108,8 @@ navbar: menu: - text: System Testing href: articles/System_Testing.html + - text: Which JT test to use + href: articles/Verification_which_JT.html - text: Study Types and Templates href: articles/Test-Guidelines.html newissue: diff --git a/figure/unnamed-chunk-3-1.png b/figure/unnamed-chunk-3-1.png new file mode 100644 index 0000000..883a853 Binary files /dev/null and b/figure/unnamed-chunk-3-1.png differ diff --git a/figure/unnamed-chunk-8-1.png b/figure/unnamed-chunk-8-1.png new file mode 100644 index 0000000..883a853 Binary files /dev/null and b/figure/unnamed-chunk-8-1.png differ diff --git a/inst/DevNotes/github_actions.qmd b/inst/DevNotes/github_actions.qmd index c18819d..9f10cd4 100644 --- a/inst/DevNotes/github_actions.qmd +++ b/inst/DevNotes/github_actions.qmd @@ -12,17 +12,14 @@ r_build_args: --no-build-vignettes --no-manual --no-resave-data r_check_args: --no-build-vignettes --no-manual ``` -R buidling warnings: - -Found the following apparent S3 methods exported but not registered: - ED.plus ED.ZG plot.edList plot.modList - See section ‘Registering S3 methods’ in the ‘Writing R Extensions’ - manual. - - Undocumented code objects: - ‘dat_bcs1’ ‘dat_medium’ ‘dat_noED50’ ‘dat_shallow’ ‘dat_steep’ - ‘exampleHistData’ ‘getLineContrast’ ‘getQuadContrast’ ‘metaldata’ - ‘NTA_Ar_ext_mortality_expected’ ‘NTP_example_rate_expected’ - ‘rankTransform’ ‘test_cases_data’ ‘test_cases_res’ - ‘williamsTestLookUpTable’ -## Undocumented data sets: + +## external images + +The second way does not work! + +``` +resource_files: + - article_assets/quantal_NOEC.png + - ../assets/binomial_tank_effects_visualization.png +``` + diff --git a/inst/RStats/RStats_Sim_Power_Ana.qmd b/inst/RStats/RStats_Sim_Power_Ana.qmd index 2dd882f..cfaad53 100644 --- a/inst/RStats/RStats_Sim_Power_Ana.qmd +++ b/inst/RStats/RStats_Sim_Power_Ana.qmd @@ -27,6 +27,561 @@ library(tidyverse) library(drcHelper) ``` + +## Discussion Outcome + +Thanks for sharing your `simulate_dose_response` function. This gives me a good starting point for generating dose-response data for testing and power analysis. + +Let me now help you develop the `combine_multiple_comparisons` function and refine the existing functions to make them work together seamlessly. I'll also set up a framework for simulation-based power analysis. + +First, let's develop the `combine_multiple_comparisons` function: + +```{r} +#' Combine Multiple Comparison Test Results for Dose-Response Analysis +#' +#' This function applies various multiple comparison methods to dose-response data +#' and returns the results in a standardized format for easy comparison. +#' +#' @param data A data frame containing the dose-response data +#' @param response The name of the response variable +#' @param dose The name of the dose variable +#' @param group An optional grouping variable name +#' @param control The dose level to use as control (default: 0) +#' @param alternative Direction of the alternative hypothesis ("two.sided", "less", or "greater") +#' @param conf.level Confidence level for intervals +#' @param methods Vector of method names to apply. Available methods include: +#' "Dunnett_multcomp", "Dunnett_DescTools", "Dunnett_PMCMRplus", +#' "Williams_PMCMRplus", "Williams_JG", "stepDown_PMCMR" +#' +#' @return A tibble containing the combined results from all specified methods +#' @export +#' +#' @importFrom dplyr bind_rows +#' @importFrom stats aov as.formula +#' @importFrom tibble tibble +combine_multiple_comparisons <- function(data, response, dose, group = NULL, + control = 0, alternative = "two.sided", + conf.level = 0.95, + methods = c("Dunnett_multcomp", "Williams_PMCMRplus")) { + + # Check inputs + if (!is.data.frame(data)) { + stop("'data' must be a data frame") + } + + if (!response %in% names(data)) { + stop("'response' variable not found in data") + } + + if (!dose %in% names(data)) { + stop("'dose' variable not found in data") + } + + if (!is.null(group) && !group %in% names(data)) { + stop("'group' variable not found in data") + } + + # Ensure dose is treated as a factor + data[[dose]] <- as.factor(data[[dose]]) + + # Find the control level index + dose_levels <- levels(data[[dose]]) + control_idx <- which(dose_levels == as.character(control)) + + if (length(control_idx) == 0) { + stop("Control dose level not found in data") + } + + # Create formula based on whether a grouping variable is provided + if (is.null(group)) { + formula <- as.formula(paste(response, "~", dose)) + } else { + formula <- as.formula(paste(response, "~", dose, "+", group)) + } + + # Fit ANOVA model + aov_model <- aov(formula, data = data) + + # Initialize results list + results_list <- list() + + # Apply each requested method + for (method in methods) { + if (grepl("^Dunnett", method)) { + dunnett_method <- gsub("^Dunnett_", "", method) + method_name <- paste0("Dunnett_", dunnett_method) + + # Apply Dunnett test with the appropriate method + result <- broom_dunnett( + aov_model, + method = method_name, + control = control_idx, + alternative = alternative, + conf.level = conf.level + ) + + } else if (grepl("^Williams", method)) { + williams_method <- gsub("^Williams_", "", method) + method_name <- paste0("Williams_", williams_method) + + # Apply Williams test with the appropriate method + result <- broom_williams( + aov_model, + method = method_name, + control = control_idx, + alternative = alternative, + conf.level = conf.level + ) + + } else if (method == "stepDown_PMCMR") { + # Apply step-down test from PMCMR package + result <- try({ + sd_test <- PMCMR::stepDownTrendTest(data[[response]], data[[dose]]) + tibble::tibble( + comparison = paste(names(sd_test$p.value), "- control"), + estimate = NA_real_, # Step-down test doesn't provide estimates + p.value = sd_test$p.value, + conf.low = NA_real_, # Step-down test doesn't provide confidence intervals + conf.high = NA_real_, + method = "stepDown_PMCMR" + ) + }, silent = TRUE) + + if (inherits(result, "try-error")) { + warning("Step-down test failed. Error: ", attr(result, "condition")$message) + result <- tibble::tibble( + comparison = character(), + estimate = numeric(), + p.value = numeric(), + conf.low = numeric(), + conf.high = numeric(), + method = character() + ) + } + } else { + warning("Unknown method: ", method) + next + } + + # Add to results list + results_list[[method]] <- result + } + + # Combine all results + combined_results <- dplyr::bind_rows(results_list) + + # Add dose information to make results more interpretable + combined_results <- combined_results %>% + dplyr::mutate( + dose_level = sub(".*?([0-9.]+).*", "\\1", comparison), + control_level = as.character(control) + ) + + return(combined_results) +} +``` + +Now, let's add the `williamsTest_JG` function that was referenced in your `broom_williams` function: + +```{r} +#' Williams Test Implementation (Jones-Graul Method) +#' +#' This function implements the Williams test for monotonic dose-response relationships +#' using the Jones-Graul method. +#' +#' @param x An object of class formula, aov, lm, or a data frame or matrix +#' @param ... Additional arguments passed to methods +#' +#' @return A data frame with standardized output columns +#' @export +williamsTest_JG <- function(x, ...) { + UseMethod("williamsTest_JG") +} + +#' @export +williamsTest_JG.formula <- function(formula, data, ...) { + model <- aov(formula, data = data) + williamsTest_JG(model, ...) +} + +#' @export +williamsTest_JG.aov <- function(x, control = 1, alternative = c("greater", "less", "two.sided"), + conf.level = 0.95, ...) { + alternative <- match.arg(alternative) + + # Extract data from the model + model_data <- model.frame(x) + response_var <- model_data[, 1] + factor_var <- model_data[, 2] + + # Get group means and sample sizes + group_means <- tapply(response_var, factor_var, mean) + group_sizes <- tapply(response_var, factor_var, length) + + # Get residual variance from ANOVA + residual_var <- summary(x)[[1]]["Residuals", "Mean Sq"] + + # Number of groups + k <- length(group_means) + + if (control < 1 || control > k) { + stop("Control group index must be between 1 and k") + } + + # Calculate isotonic means (step-down approach) + isotonic_means <- numeric(k) + isotonic_means[control] <- group_means[control] + + # For doses higher than control + if (control < k) { + for (i in (control + 1):k) { + # Pool means if current mean is less than previous isotonic mean + if (i == control + 1) { + isotonic_means[i] <- max(group_means[i], group_means[control]) + } else { + isotonic_means[i] <- max(group_means[i], isotonic_means[i - 1]) + } + } + } + + # For doses lower than control (if any) + if (control > 1) { + for (i in (control - 1):1) { + # Pool means if current mean is greater than next isotonic mean + if (i == control - 1) { + isotonic_means[i] <- min(group_means[i], group_means[control]) + } else { + isotonic_means[i] <- min(group_means[i], isotonic_means[i + 1]) + } + } + } + + # Calculate test statistics + t_stats <- numeric(k - 1) + comparisons <- character(k - 1) + estimates <- numeric(k - 1) + p_values <- numeric(k - 1) + conf_low <- numeric(k - 1) + conf_high <- numeric(k - 1) + + idx <- 1 + for (i in setdiff(1:k, control)) { + # Calculate contrast weights + weights <- rep(0, k) + weights[control] <- -1 + weights[i] <- 1 + + # Calculate contrast + contrast <- sum(weights * isotonic_means) + + # Calculate standard error + se <- sqrt(residual_var * sum(weights^2 / group_sizes)) + + # Calculate t-statistic + t_stat <- contrast / se + t_stats[idx] <- t_stat + + # Calculate p-value based on alternative + if (alternative == "two.sided") { + p_value <- 2 * pt(-abs(t_stat), df = x$df.residual) + } else if (alternative == "greater") { + p_value <- pt(-t_stat, df = x$df.residual) + } else { + p_value <- pt(t_stat, df = x$df.residual) + } + + # Calculate confidence intervals + t_crit <- qt(1 - (1 - conf.level) / 2, df = x$df.residual) + margin <- t_crit * se + + # Store results + group_levels <- levels(factor_var) + comparisons[idx] <- paste(group_levels[i], "-", group_levels[control]) + estimates[idx] <- contrast + p_values[idx] <- p_value + conf_low[idx] <- contrast - margin + conf_high[idx] <- contrast + margin + + idx <- idx + 1 + } + + # Create result data frame + result <- tibble::tibble( + comparison = comparisons, + estimate = estimates, + p.value = p_values, + conf.low = conf_low, + conf.high = conf_high, + method = "Williams_JG" + ) + + return(result) +} + +#' @export +williamsTest_JG.default <- function(x, g, control = 1, ...) { + # Convert to formula interface + df <- data.frame(response = x, group = factor(g)) + williamsTest_JG(response ~ group, data = df, control = control, ...) +} +``` + +Now, let's create a function for simulation-based power analysis: + +```{r} +#' Perform Simulation-Based Power Analysis for Dose-Response Tests +#' +#' This function conducts a simulation-based power analysis for various multiple +#' comparison methods in dose-response analysis. +#' +#' @param n_sims Number of simulations to run +#' @param n_doses Number of dose levels +#' @param dose_range Vector of length 2 specifying the min and max dose values +#' @param m_tanks Number of tanks per dose +#' @param k_individuals Number of individuals per tank +#' @param var_tank Variance at the tank level +#' @param var_individual Variance at the individual level +#' @param response_function Function that calculates the response given a dose +#' @param alpha Significance level for hypothesis tests +#' @param methods Vector of method names to compare +#' @param include_individuals Whether to simulate individual-level data +#' @param ... Additional parameters passed to response_function +#' +#' @return A tibble containing power analysis results for each method and dose level +#' @export +#' +#' @importFrom dplyr bind_rows group_by summarize +#' @importFrom purrr map_dfr +#' @importFrom tibble tibble +power_analysis_dose_response <- function(n_sims = 100, + n_doses = 5, + dose_range = c(0, 20), + m_tanks = 3, + k_individuals = 10, + var_tank = 4, + var_individual = 2, + response_function = NULL, + alpha = 0.05, + methods = c("Dunnett_multcomp", "Williams_PMCMRplus"), + include_individuals = TRUE, + ...) { + + # Set default response function if not provided + if (is.null(response_function)) { + response_function <- function(dose, lower = 0, upper = 100, ED50 = 10, slope = 1) { + lower + (upper - lower) / (1 + exp(-slope * (dose - ED50))) + } + } + + # Function to run a single simulation + run_simulation <- function(sim_id) { + # Generate data + sim_data <- simulate_dose_response( + n_doses = n_doses, + dose_range = dose_range, + m_tanks = m_tanks, + k_individuals = k_individuals, + var_tank = var_tank, + var_individual = var_individual, + include_individuals = include_individuals, + response_function = response_function, + ... + ) + + # If we have individual-level data, aggregate to tank level + if (include_individuals) { + sim_data <- sim_data %>% + dplyr::group_by(Dose, Tank) %>% + dplyr::summarize(Response = mean(Response), .groups = "drop") + } + + # Run multiple comparison tests + test_results <- try({ + combine_multiple_comparisons( + data = sim_data, + response = "Response", + dose = "Dose", + control = min(sim_data$Dose), + alternative = "two.sided", + conf.level = 1 - alpha, + methods = methods + ) + }, silent = TRUE) + + if (inherits(test_results, "try-error")) { + warning("Error in simulation ", sim_id, ": ", attr(test_results, "condition")$message) + return(NULL) + } + + # Add simulation ID + test_results$simulation <- sim_id + + return(test_results) + } + + # Run all simulations + message("Running ", n_sims, " simulations...") + all_results <- purrr::map_dfr(1:n_sims, run_simulation) + + # Calculate power for each method and dose level + if (nrow(all_results) > 0) { + power_results <- all_results %>% + dplyr::group_by(method, comparison) %>% + dplyr::summarize( + n_significant = sum(p.value < alpha, na.rm = TRUE), + n_tests = sum(!is.na(p.value)), + power = n_significant / n_tests, + mean_estimate = mean(estimate, na.rm = TRUE), + sd_estimate = sd(estimate, na.rm = TRUE), + .groups = "drop" + ) + + return(power_results) + } else { + warning("No valid results from simulations") + return(tibble::tibble( + method = character(), + comparison = character(), + n_significant = integer(), + n_tests = integer(), + power = numeric(), + mean_estimate = numeric(), + sd_estimate = numeric() + )) + } +} +``` + +Let's also create a visualization function to help interpret the power analysis results: + +```{r} +#' Visualize Power Analysis Results for Dose-Response Tests +#' +#' This function creates visualizations of power analysis results from +#' multiple comparison methods in dose-response analysis. +#' +#' @param power_results A tibble containing power analysis results from power_analysis_dose_response +#' @param plot_type Type of plot to create: "power" (default), "estimates", or "both" +#' @param arrange_by How to arrange facets: "method" (default) or "dose" +#' +#' @return A ggplot2 object +#' @export +#' +#' @importFrom ggplot2 ggplot aes geom_point geom_line facet_wrap theme_minimal labs +#' @importFrom dplyr mutate +#' @importFrom tidyr extract +visualize_power_analysis <- function(power_results, + plot_type = c("power", "estimates", "both"), + arrange_by = c("method", "dose")) { + + plot_type <- match.arg(plot_type) + arrange_by <- match.arg(arrange_by) + + # Extract dose levels from comparison strings + power_results <- power_results %>% + tidyr::extract(comparison, + into = c("dose_level", "control_level"), + regex = "([0-9.]+) - ([0-9.]+)", + remove = FALSE, + convert = TRUE) + + if (plot_type == "power" || plot_type == "both") { + # Create power plot + power_plot <- ggplot2::ggplot(power_results, ggplot2::aes(x = dose_level, y = power, color = method)) + + ggplot2::geom_point() + + ggplot2::geom_line() + + ggplot2::geom_hline(yintercept = 0.8, linetype = "dashed", color = "gray50") + + ggplot2::theme_minimal() + + ggplot2::labs( + title = "Power Analysis for Dose-Response Tests", + x = "Dose Level", + y = "Statistical Power", + color = "Test Method" + ) + + if (arrange_by == "method") { + power_plot <- power_plot + ggplot2::facet_wrap(~ method) + } else { + power_plot <- power_plot + ggplot2::facet_wrap(~ dose_level) + } + } + + if (plot_type == "estimates" || plot_type == "both") { + # Create estimates plot + estimates_plot <- ggplot2::ggplot(power_results, + ggplot2::aes(x = dose_level, y = mean_estimate, color = method)) + + ggplot2::geom_point() + + ggplot2::geom_errorbar(ggplot2::aes(ymin = mean_estimate - sd_estimate, + ymax = mean_estimate + sd_estimate), width = 0.2) + + ggplot2::theme_minimal() + + ggplot2::labs( + title = "Effect Estimates for Dose-Response Tests", + x = "Dose Level", + y = "Mean Effect Estimate", + color = "Test Method" + ) + + if (arrange_by == "method") { + estimates_plot <- estimates_plot + ggplot2::facet_wrap(~ method) + } else { + estimates_plot <- estimates_plot + ggplot2::facet_wrap(~ dose_level) + } + } + + # Return appropriate plot(s) + if (plot_type == "power") { + return(power_plot) + } else if (plot_type == "estimates") { + return(estimates_plot) + } else { + return(list(power = power_plot, estimates = estimates_plot)) + } +} +``` + +Now let's create a simple example to demonstrate how these functions work together: + +```{r} +#' @examples +#' # Define a dose-response function with a threshold effect at dose 5 +#' threshold_response <- function(dose, threshold = 5, effect_size = 20) { +#' ifelse(dose > threshold, effect_size, 0) +#' } +#' +#' # Run a small power analysis +#' power_results <- power_analysis_dose_response( +#' n_sims = 50, +#' n_doses = 5, +#' dose_range = c(0, 10), +#' m_tanks = 3, +#' k_individuals = 8, +#' var_tank = 5, +#' var_individual = 2, +#' response_function = threshold_response, +#' threshold = 5, +#' effect_size = 15, +#' methods = c("Dunnett_multcomp", "Williams_PMCMRplus", "stepDown_PMCMR") +#' ) +#' +#' # Visualize the results +#' visualize_power_analysis(power_results, plot_type = "both") +``` + +These functions should provide a solid foundation for your dose-response analysis package. The `combine_multiple_comparisons` function integrates the results from various testing approaches, and the `power_analysis_dose_response` function allows you to compare the statistical power of these methods through simulation. + +Let me know if you'd like me to modify any aspects of these functions or if you need additional functionality! + + + + + + + + + + + + ## Simulate Data ```{r} diff --git a/vignettes/Dunnetts_Test_for_Data_with_Hierarchical_Structure_files/figure-html/unnamed-chunk-2-1.png b/vignettes/Dunnetts_Test_for_Data_with_Hierarchical_Structure_files/figure-html/unnamed-chunk-2-1.png deleted file mode 100644 index eb97913..0000000 Binary files a/vignettes/Dunnetts_Test_for_Data_with_Hierarchical_Structure_files/figure-html/unnamed-chunk-2-1.png and /dev/null differ diff --git a/vignettes/articles/NOEC_ECx_BMD.Rmd b/vignettes/articles/NOEC_ECx_BMD.Rmd index d812386..bf56db9 100644 --- a/vignettes/articles/NOEC_ECx_BMD.Rmd +++ b/vignettes/articles/NOEC_ECx_BMD.Rmd @@ -51,44 +51,44 @@ Current statistical practices typically involve the following steps: ### NOEC -– Definition: The highest tested concentration at which no statistically significant adverse effect is observed relative to the control. -– Regulatory Role: Historically used as the primary benchmark (e.g., for deriving acceptable daily intakes or environmental quality standards). +- Definition: The highest tested concentration at which no statistically significant adverse effect is observed relative to the control. +- Regulatory Role: Historically used as the primary benchmark (e.g., for deriving acceptable daily intakes or environmental quality standards). **Strengths:** -– Straightforward concept and quick to interpret -– Familiar to risk assessors and historically embedded in many regulatory frameworks +- Straightforward concept and quick to interpret +- Familiar to risk assessors and historically embedded in many regulatory frameworks **Weaknesses:** -– Strongly dependent on the dose spacing and study design -– Does not provide information on the magnitude of effects or their progression along a dose–response curve -– Binary interpretation ("effect" vs. "no effect") can oversimplify data -– Potentially influenced by low statistical power, especially in studies with small sample sizes +- Strongly dependent on the dose spacing and study design +- Does not provide information on the magnitude of effects or their progression along a dose–response curve +- Binary interpretation ("effect" vs. "no effect") can oversimplify data +- Potentially influenced by low statistical power, especially in studies with small sample sizes ### ECx and BMD -– Definition: The concentration at which a predefined percentage (x) of effect (e.g., 10\%, 20\% reduction in a biological endpoint) is observed relative to the control. +- Definition: The concentration at which a predefined percentage (x) of effect (e.g., 10\%, 20\% reduction in a biological endpoint) is observed relative to the control. -– Regulatory Role: Increasingly used to quantitatively describe the magnitude of adverse effects along the dose–response curve and considered more informative in deriving species sensitivity distributions. +- Regulatory Role: Increasingly used to quantitatively describe the magnitude of adverse effects along the dose–response curve and considered more informative in deriving species sensitivity distributions. **Strengths:** -– Utilizes the full data set by modeling the dose–response relationship -– Provides a more nuanced and continuous estimate along the effect continuum (e.g., EC10 is often seen as protective) -– Facilitates interpolation and extrapolation within the tested range -– Offers the possibility to quantify uncertainty via confidence intervals +- Utilizes the full data set by modeling the dose–response relationship +- Provides a more nuanced and continuous estimate along the effect continuum (e.g., EC10 is often seen as protective) +- Facilitates interpolation and extrapolation within the tested range +- Offers the possibility to quantify uncertainty via confidence intervals **Weaknesses:** -– Heavily dependent on the choice and fit of the dose–response model -– Requires more sophisticated statistical methods and expertise -– May be sensitive to violations of model assumptions if data quality or distribution is suboptimal -– In some cases, different models may yield different ECx estimates, raising issues of model selection and validation +- Heavily dependent on the choice and fit of the dose–response model +- Requires more sophisticated statistical methods and expertise +- May be sensitive to violations of model assumptions if data quality or distribution is suboptimal +- In some cases, different models may yield different ECx estimates, raising issues of model selection and validation ## Methodological Improvements -SETAC’s journals (Environmental Toxicology and Chemistry) regularly publish papers on methodological improvements in statistical evaluation of ecotoxicological studies. +Papers, scientific opinions, guidance are regularly published on methodological improvements in statistical evaluation of ecotoxicological studies. - We advocate for transparent reporting of uncertainties and assumptions in model selection or model averaging. - Integration of Bayesian methods and bootstrap techniques to better quantify uncertainty is emerging in research discussions. diff --git a/vignettes/articles/Quantal-Data.Rmd b/vignettes/articles/Quantal-Data.Rmd index 02ca03e..545543f 100644 --- a/vignettes/articles/Quantal-Data.Rmd +++ b/vignettes/articles/Quantal-Data.Rmd @@ -2,7 +2,7 @@ title: "Quantal Data" resource_files: - article_assets/quantal_NOEC.png - - ../assets/binomial_tank_effects_visualization.png + - article_assets/binomial_tank_effects_visualization.png editor_options: chunk_output_type: console --- @@ -23,7 +23,7 @@ library(drcHelper) Analysis of quantal data typically involves the following steps as shown in the flowchart below: -![Quantal NOEC Flowchart](https://bayer-group.github.io/drcHelper/articles/assets/quantal_NOEC.png) + ```{r echo=FALSE,out.width="90%"} #| fig.alt: > @@ -439,7 +439,7 @@ ggsave("binomial_tank_effects_visualization.png", combined_plot, width = 12, hei ``` -```{r echo=FALSE} +```{r echo=FALSE,eval=FALSE,include=FALSE} knitr::include_graphics("./../assets/binomial_tank_effects_visualization.png") ``` diff --git a/vignettes/articles/Test-Guidelines.Rmd b/vignettes/articles/Test-Guidelines.Rmd index 533ca8b..dffeef0 100644 --- a/vignettes/articles/Test-Guidelines.Rmd +++ b/vignettes/articles/Test-Guidelines.Rmd @@ -5,7 +5,9 @@ title: "Study Types and Templates" ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, - comment = "#>" + comment = "#>", + warning = FALSE, + message = FALSE ) library(tidyverse) ``` diff --git a/vignettes/articles/Verification_which_JT.Rmd b/vignettes/articles/Verification_which_JT.Rmd new file mode 100644 index 0000000..9327f7f --- /dev/null +++ b/vignettes/articles/Verification_which_JT.Rmd @@ -0,0 +1,158 @@ +--- +title: "Verification: Which Jonckeere Terpstra Test to Use" +editor_options: + chunk_output_type: console +--- + + + + +```r +library(drcHelper) +``` + + + +```r + +library(tidyverse) +``` + +## Summary + +`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 +``` + diff --git a/vignettes/articles/Verification_which_JT.Rmd.orig b/vignettes/articles/Verification_which_JT.Rmd.orig new file mode 100644 index 0000000..d0a2cfa --- /dev/null +++ b/vignettes/articles/Verification_which_JT.Rmd.orig @@ -0,0 +1,86 @@ +--- +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 + +`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) +``` +