diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
index d564d91..5ce0fa4 100644
--- a/.github/workflows/pkgdown.yaml
+++ b/.github/workflows/pkgdown.yaml
@@ -33,7 +33,7 @@ jobs:
 
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: any::pkgdown, any::tidyverse, any::drda, any::pander, any::metafor, any::kableExtra, any::plyr, any::ggprism, any::ggnewscale, any::GLMcat, any::quarto,any::Hmisc,any::outliers, any::multcomp, any::emmeans, local::.
+          extra-packages: any::pkgdown, any::tidyverse, any::drda, any::pander, any::metafor, any::kableExtra, any::plyr, any::ggprism, any::ggnewscale, any::GLMcat, any::quarto,any::Hmisc,any::outliers, any::multcomp, any::emmeans, any::patchwork, local::.
           needs: website
 
       - name: Build site
diff --git a/NAMESPACE b/NAMESPACE
index f64fca7..61d5fcc 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,8 +1,12 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(calpha.test,fisher)
+S3method(plot,StepDownRSCABS)
+S3method(print,RSCABS)
+S3method(print,StepDownRSCABS)
 S3method(print,drcComp)
 S3method(print,tskresult)
+S3method(summary,StepDownRSCABS)
 S3method(tsk,data.frame)
 S3method(tsk,numeric)
 export("%>%")
@@ -49,10 +53,13 @@ export(rankTransform)
 export(reshape_drcData)
 export(runRSCABS)
 export(run_RSCA)
+export(run_all_threshold_tests)
+export(run_threshold_RSCA)
 export(simDRdata)
 export(simplifyTreatment)
 export(stepDownRSCABS)
 export(stepKRSCABS)
+export(step_down_RSCABS)
 export(summaryZG)
 export(treatment2dose)
 export(tsk)
@@ -60,8 +67,31 @@ export(williamsTest_JG)
 import(dplyr)
 import(ggplot2)
 import(scales)
+importFrom(dplyr,across)
+importFrom(dplyr,arrange)
+importFrom(dplyr,bind_rows)
+importFrom(dplyr,filter)
+importFrom(dplyr,group_by)
+importFrom(dplyr,left_join)
+importFrom(dplyr,mutate)
+importFrom(dplyr,select)
+importFrom(dplyr,summarize)
+importFrom(dplyr,where)
+importFrom(ggplot2,aes)
+importFrom(ggplot2,element_blank)
+importFrom(ggplot2,element_text)
+importFrom(ggplot2,geom_text)
+importFrom(ggplot2,geom_tile)
+importFrom(ggplot2,ggplot)
+importFrom(ggplot2,labs)
+importFrom(ggplot2,scale_fill_gradient)
+importFrom(ggplot2,scale_y_reverse)
+importFrom(ggplot2,theme)
+importFrom(ggplot2,theme_minimal)
 importFrom(isotone,gpava)
 importFrom(magrittr,"%>%")
 importFrom(metafor,rma.mh)
 importFrom(purrr,map2)
+importFrom(stats,aggregate)
+importFrom(stats,fisher.test)
 importFrom(stats,pnorm)
diff --git a/R/RSCABS_AO.R b/R/RSCABS_AO.R
index 1675014..08df6fa 100644
--- a/R/RSCABS_AO.R
+++ b/R/RSCABS_AO.R
@@ -76,6 +76,8 @@ get_RS_adj_val <- function(group, replicate, affected, total) {
 #'
 #' @return The Z-statistic for the Cochran-Armitage trend test
 #'
+#' @author Originally by Allen Olmstead
+#'
 #' @details
 #' The Cochran-Armitage trend test examines whether there is a linear trend in
 #' proportions across ordered categories (treatment groups). This implementation
@@ -83,7 +85,7 @@ get_RS_adj_val <- function(group, replicate, affected, total) {
 #'
 #' The function assigns scores (1, 2, 3, ...) to the treatment groups and
 #' calculates the Z-statistic using the formula:
-#' Z = [sum(adj_x*d) - N*p_bar*d_bar] / sqrt[p_bar*(1-p_bar)*(sum(adj_n*d^2) - N*d_bar^2)]
+#' ``\deqn{Z = [sum(adj_x*d) - N*p_bar*d_bar] / sqrt[p_bar*(1-p_bar)*(sum(adj_n*d^2) - N*d_bar^2)]``
 #' where:
 #' - d are the scores (1, 2, 3, ...)
 #' - N is the total adjusted sample size
@@ -124,11 +126,12 @@ get_CA_Z <- function(adj_x, adj_n) {
 #' @param replicate Vector of replicate/tank identifiers within treatment groups
 #' @param affected Vector of counts of affected subjects (fish with injuries) in each replicate
 #' @param total Vector of total subjects (fish) in each replicate
+#' @param correction continuity correction when there is 1, default is 0, can be changed to 0.5.
 #'
 #' @return A list containing:
 #'   \item{interm_values}{A tibble with intermediate values from the Rao-Scott adjustment}
 #'   \item{Z}{The Z-statistic for the Cochran-Armitage trend test}
-#'
+#' @author Originally by Allen Olmstead
 #' @details
 #' This function combines the Rao-Scott adjustment and the Cochran-Armitage trend test
 #' to analyze dose-response relationships in clustered data. It first calculates
@@ -156,8 +159,770 @@ get_CA_Z <- function(adj_x, adj_n) {
 #' # Calculate p-value
 #' p_value <- 2 * (1 - pnorm(abs(result$Z)))
 #' print(paste("p-value:", p_value))
-run_RSCA <- function(group, replicate, affected, total) {
+run_RSCA <- function(group, replicate, affected, total,correction=0) {
+  affected <- affected + correction
+  total <- total + correction
   interm_values <- get_RS_adj_val(group, replicate, affected, total)
   Z <- get_CA_Z(interm_values$x_tilde, interm_values$n_tilde)
   list(interm_values = interm_values, Z = Z)
 }
+
+
+
+
+#' Run Rao-Scott Adjusted Cochran-Armitage Trend Test for a Specific Injury Threshold
+#'
+#' This function performs the Rao-Scott adjusted Cochran-Armitage trend test for
+#' a specific injury threshold, counting fish with injuries at or above that threshold
+#' as "affected".
+#'
+#' @param data A data frame containing fish injury data
+#' @param threshold Numeric value indicating the injury threshold (e.g., 1 for S1+)
+#' @param score_cols Character vector of column names containing injury scores (default: NULL, auto-detected)
+#' @param treatment_col Name of the column containing treatment groups (default: "tmt")
+#' @param replicate_col Name of the column containing tank/replicate IDs (default: "tank")
+#' @param total_col Name of the column containing total counts (default: "total")
+#' @param direction Character string indicating threshold direction: "greater" for ≥ threshold,
+#'        "less" for ≤ threshold (default: "greater")
+#' @param alternative Character string specifying the alternative hypothesis:
+#'        "two.sided" (default), "greater" (proportion increases with treatment level),
+#'        or "less" (proportion decreases with treatment level)
+#'
+#' @return A list containing:
+#'   \item{threshold}{The numeric threshold value}
+#'   \item{threshold_label}{Character label for the threshold (e.g., "S1+")}
+#'   \item{Z}{Z-statistic from the RSCA test}
+#'   \item{p_value}{P-value based on the specified alternative hypothesis}
+#'   \item{alternative}{The alternative hypothesis used}
+#'   \item{interm_values}{Intermediate values from the Rao-Scott adjustment}
+#'   \item{has_zero_counts}{Logical indicating if any treatment group had zero affected individuals}
+#'   \item{affected_counts}{Data frame showing affected counts by treatment group}
+#'
+#' @examples
+#' # Example data
+#' fish_data <- data.frame(
+#'   tmt = rep(c("Control", "Low", "Medium", "High"), each = 3),
+#'   tank = paste0(rep(c("Control", "Low", "Medium", "High"), each = 3),
+#'                 rep(1:3, times = 4)),
+#'   S0 = c(8,7,9, 6,5,7, 4,5,3, 2,3,1),
+#'   S1 = c(2,2,1, 3,4,2, 4,3,5, 4,3,5),
+#'   S2 = c(0,1,0, 1,1,1, 2,2,1, 3,3,3),
+#'   S3 = c(0,0,0, 0,0,0, 0,0,1, 1,1,1),
+#'   S4 = c(0,0,0, 0,0,0, 0,0,0, 0,0,0),
+#'   total = rep(10, 12)
+#' )
+#'
+#' # Two-sided test for trend in S2+ injuries
+#' result_two_sided <- run_threshold_RSCA(fish_data, threshold = 2)
+#'
+#' # One-sided test (proportion increases with treatment level)
+#' result_greater <- run_threshold_RSCA(fish_data, threshold = 2,
+#'                                     alternative = "greater")
+#'
+#' # One-sided test (proportion decreases with treatment level)
+#' result_less <- run_threshold_RSCA(fish_data, threshold = 2,
+#'                                  alternative = "less")
+#'
+#' @importFrom dplyr group_by summarize mutate select left_join
+#' @importFrom stats pnorm
+#' @export
+run_threshold_RSCA <- function(data, threshold,
+                               score_cols = NULL,
+                               treatment_col = "tmt",
+                               replicate_col = "tank",
+                               total_col = "total",
+                               direction = "greater",
+                               alternative = "two.sided") {
+
+  # Validate alternative parameter
+  if (!alternative %in% c("two.sided", "greater", "less")) {
+    stop("alternative must be one of: 'two.sided', 'greater', or 'less'")
+  }
+
+  # Auto-detect score columns if not provided
+  if (is.null(score_cols)) {
+    score_cols <- grep("^S[0-9]+$", names(data), value = TRUE)
+    if (length(score_cols) == 0) {
+      stop("No score columns found with pattern 'S[0-9]+'. Please specify score_cols manually.")
+    }
+  }
+
+  # Validate inputs
+  if (!treatment_col %in% names(data)) {
+    stop("Treatment column '", treatment_col, "' not found in data")
+  }
+  if (!replicate_col %in% names(data)) {
+    stop("Replicate column '", replicate_col, "' not found in data")
+  }
+  if (!total_col %in% names(data)) {
+    stop("Total column '", total_col, "' not found in data")
+  }
+  if (!all(score_cols %in% names(data))) {
+    missing <- setdiff(score_cols, names(data))
+    stop("Missing score columns: ", paste(missing, collapse = ", "))
+  }
+  if (!direction %in% c("greater", "less")) {
+    stop("Direction must be either 'greater' or 'less'")
+  }
+
+  # Extract numeric values from score column names
+  score_values <- as.numeric(gsub("\\D", "", score_cols))
+
+  # Determine which scores to include based on direction
+  if (direction == "greater") {
+    # Greater than or equal to threshold
+    selected_cols <- score_cols[score_values >= threshold]
+    threshold_label <- paste0("S", threshold, "+")
+  } else {
+    # Less than or equal to threshold
+    selected_cols <- score_cols[score_values <= threshold]
+    threshold_label <- paste0("S≤", threshold)
+  }
+
+  # Calculate affected counts
+  affected <- rowSums(data[, selected_cols, drop = FALSE])
+
+  # Check for zero counts in any treatment group
+  has_zero_counts <- any(tapply(affected, data[[treatment_col]], function(x) all(x == 0)))
+
+  # Calculate affected counts by treatment for reporting
+  affected_counts <- aggregate(
+    affected ~ data[[treatment_col]],
+    FUN = sum,
+    data = data.frame(affected = affected)
+  )
+  names(affected_counts) <- c("treatment", "affected")
+
+  # If any treatment group has all zeros, return NA for Z and p-value
+  if (has_zero_counts) {
+    warning("Some treatment groups have zero affected individuals at threshold ",
+            threshold_label, ". RSCA test may not be valid.")
+
+    # Try to run the test anyway, but catch errors
+    tryCatch({
+      # Run the RSCA test
+      result <- run_RSCA(
+        group = data[[treatment_col]],
+        replicate = data[[replicate_col]],
+        affected = affected,
+        total = data[[total_col]]
+      )
+
+      # Calculate p-value based on alternative hypothesis
+      Z <- result$Z
+      p_value <- switch(alternative,
+                        "two.sided" = 2 * pnorm(-abs(Z)),
+                        "greater" = pnorm(-Z),  # H1: proportion increases with treatment level
+                        "less" = pnorm(Z)       # H1: proportion decreases with treatment level
+      )
+
+      # Return results with warning flag
+      return(list(
+        threshold = threshold,
+        threshold_label = threshold_label,
+        Z = Z,
+        p_value = p_value,
+        alternative = alternative,
+        interm_values = result$interm_values,
+        has_zero_counts = TRUE,
+        affected_counts = affected_counts
+      ))
+    }, error = function(e) {
+      # Return NA results if test fails
+      return(list(
+        threshold = threshold,
+        threshold_label = threshold_label,
+        Z = NA,
+        p_value = NA,
+        alternative = alternative,
+        interm_values = NULL,
+        has_zero_counts = TRUE,
+        affected_counts = affected_counts,
+        error = e$message
+      ))
+    })
+  } else {
+    # Run the RSCA test
+    result <- run_RSCA(
+      group = data[[treatment_col]],
+      replicate = data[[replicate_col]],
+      affected = affected,
+      total = data[[total_col]]
+    )
+
+    # Calculate p-value based on alternative hypothesis
+    Z <- result$Z
+    p_value <- switch(alternative,
+                      "two.sided" = 2 * pnorm(-abs(Z)),
+                      "greater" = pnorm(-Z),  # H1: proportion increases with treatment level
+                      "less" = pnorm(Z)       # H1: proportion decreases with treatment level
+    )
+
+    # Return results
+    return(list(
+      threshold = threshold,
+      threshold_label = threshold_label,
+      Z = Z,
+      p_value = p_value,
+      alternative = alternative,
+      interm_values = result$interm_values,
+      has_zero_counts = FALSE,
+      affected_counts = affected_counts
+    ))
+  }
+}
+
+#' Run Rao-Scott Adjusted Cochran-Armitage Trend Tests for All Thresholds
+#'
+#' This function performs the Rao-Scott adjusted Cochran-Armitage trend test for
+#' multiple injury thresholds, providing a comprehensive analysis of dose-response
+#' relationships at different severity levels.
+#'
+#' @param data A data frame containing fish injury data
+#' @param max_score Maximum score value to consider (default: NULL, auto-detected)
+#' @param min_score Minimum score value to consider (default: 1)
+#' @param score_cols Character vector of column names containing injury scores (default: NULL, auto-detected)
+#' @param treatment_col Name of the column containing treatment groups (default: "tmt")
+#' @param replicate_col Name of the column containing tank/replicate IDs (default: "tank")
+#' @param total_col Name of the column containing total counts (default: "total")
+#' @param direction Character string indicating threshold direction: "greater" for ≥ threshold,
+#'        "less" for ≤ threshold (default: "greater")
+#' @param alternative Character string specifying the alternative hypothesis:
+#'        "two.sided" (default), "greater" (proportion increases with treatment level),
+#'        or "less" (proportion decreases with treatment level)
+#' @param include_fisher Logical indicating whether to use Fisher's exact test when RSCA fails (default: TRUE)
+#'
+#' @return A list containing:
+#'   \item{results}{Data frame with test results for each threshold}
+#'   \item{proportions}{Data frame with proportions of affected fish by treatment and threshold}
+#'   \item{detailed_results}{List of detailed results for each threshold}
+#'
+#' @examples
+#' # Example data
+#' fish_data <- data.frame(
+#'   tmt = rep(c("Control", "Low", "Medium", "High"), each = 3),
+#'   tank = paste0(rep(c("Control", "Low", "Medium", "High"), each = 3),
+#'                 rep(1:3, times = 4)),
+#'   S0 = c(8,7,9, 6,5,7, 4,5,3, 2,3,1),
+#'   S1 = c(2,2,1, 3,4,2, 4,3,5, 4,3,5),
+#'   S2 = c(0,1,0, 1,1,1, 2,2,1, 3,3,3),
+#'   S3 = c(0,0,0, 0,0,0, 0,0,1, 1,1,1),
+#'   S4 = c(0,0,0, 0,0,0, 0,0,0, 0,0,0),
+#'   total = rep(10, 12)
+#' )
+#'
+#' # Run two-sided tests for all thresholds
+#' all_results <- run_all_threshold_tests(fish_data)
+#'
+#' # Run one-sided tests (proportion increases with treatment level)
+#' all_results_greater <- run_all_threshold_tests(
+#'   fish_data,
+#'   alternative = "greater"
+#' )
+#'
+#' # View results table
+#' print(all_results$results)
+#' print(all_results_greater$results)
+#'
+#' @importFrom dplyr group_by summarize mutate select across where
+#' @importFrom stats fisher.test
+#' @export
+run_all_threshold_tests <- function(data,
+                                    max_score = NULL,
+                                    min_score = 1,
+                                    score_cols = NULL,
+                                    treatment_col = "tmt",
+                                    replicate_col = "tank",
+                                    total_col = "total",
+                                    direction = "greater",
+                                    alternative = "two.sided",
+                                    include_fisher = TRUE) {
+
+  # Validate alternative parameter
+  if (!alternative %in% c("two.sided", "greater", "less")) {
+    stop("alternative must be one of: 'two.sided', 'greater', or 'less'")
+  }
+
+  # Auto-detect score columns if not provided
+  if (is.null(score_cols)) {
+    score_cols <- grep("^S[0-9]+$", names(data), value = TRUE)
+    if (length(score_cols) == 0) {
+      stop("No score columns found with pattern 'S[0-9]+'. Please specify score_cols manually.")
+    }
+  }
+
+  # Auto-detect max_score if not provided
+  if (is.null(max_score)) {
+    score_values <- as.numeric(gsub("\\D", "", score_cols))
+    max_score <- max(score_values)
+  }
+
+  # Create a list to store results
+  results_list <- list()
+
+  # Run tests for each threshold
+  thresholds <- min_score:max_score
+  for (i in seq_along(thresholds)) {
+    threshold <- thresholds[i]
+    results_list[[i]] <- run_threshold_RSCA(
+      data = data,
+      threshold = threshold,
+      score_cols = score_cols,
+      treatment_col = treatment_col,
+      replicate_col = replicate_col,
+      total_col = total_col,
+      direction = direction,
+      alternative = alternative
+    )
+  }
+
+  # Compile results into a data frame
+  results_df <- data.frame(
+    Threshold = sapply(results_list, function(x) x$threshold_label),
+    Z_statistic = sapply(results_list, function(x) x$Z),
+    P_value = sapply(results_list, function(x) x$p_value),
+    Has_zero_counts = sapply(results_list, function(x) x$has_zero_counts),
+    Alternative = sapply(results_list, function(x) x$alternative),
+    stringsAsFactors = FALSE
+  )
+
+  # Add Fisher's exact test for thresholds with zero counts if requested
+  if (include_fisher) {
+    for (i in seq_along(thresholds)) {
+      if (results_list[[i]]$has_zero_counts && is.na(results_list[[i]]$p_value)) {
+        threshold <- thresholds[i]
+
+        # Determine which scores to include based on direction
+        score_values <- as.numeric(gsub("\\D", "", score_cols))
+        if (direction == "greater") {
+          selected_cols <- score_cols[score_values >= threshold]
+        } else {
+          selected_cols <- score_cols[score_values <= threshold]
+        }
+
+        # Calculate affected counts
+        affected <- rowSums(data[, selected_cols, drop = FALSE])
+        unaffected <- data[[total_col]] - affected
+
+        # Aggregate by treatment
+        agg_data <- aggregate(
+          cbind(affected, unaffected) ~ data[[treatment_col]],
+          FUN = sum
+        )
+        names(agg_data)[1] <- "treatment"
+
+        # Run Fisher's exact test
+        tryCatch({
+          cont_matrix <- as.matrix(agg_data[, c("affected", "unaffected")])
+          rownames(cont_matrix) <- agg_data$treatment
+
+          # Adjust Fisher's test based on alternative hypothesis
+          fisher_result <- switch(alternative,
+                                  "two.sided" = fisher.test(cont_matrix, alternative = "two.sided"),
+                                  "greater" = fisher.test(cont_matrix, alternative = "greater"),
+                                  "less" = fisher.test(cont_matrix, alternative = "less")
+          )
+
+          # Update results
+          results_df$P_value[i] <- fisher_result$p.value
+          results_df$Method[i] <- paste0("Fisher's Exact Test (", alternative, ")")
+        }, error = function(e) {
+          results_df$Method[i] <- "Not calculable"
+        })
+      } else {
+        results_df$Method[i] <- paste0("RSCA (", alternative, ")")
+      }
+    }
+  }
+
+  # Calculate proportions for each threshold by treatment
+  prop_by_treatment <- data.frame(Treatment = unique(data[[treatment_col]]))
+
+  for (threshold in thresholds) {
+    # Determine which scores to include based on direction
+    score_values <- as.numeric(gsub("\\D", "", score_cols))
+    if (direction == "greater") {
+      selected_cols <- score_cols[score_values >= threshold]
+      threshold_label <- paste0("S", threshold, "+")
+    } else {
+      selected_cols <- score_cols[score_values <= threshold]
+      threshold_label <- paste0("S≤", threshold)
+    }
+
+    # Calculate proportions by treatment
+    props <- aggregate(
+      cbind(affected = rowSums(data[, selected_cols, drop = FALSE]),
+            total = data[[total_col]]) ~ data[[treatment_col]],
+      FUN = sum
+    )
+    names(props)[1] <- "Treatment"
+    props$proportion <- props$affected / props$total
+
+    # Add to the data frame
+    prop_by_treatment[[threshold_label]] <-
+      props$proportion[match(prop_by_treatment$Treatment, props$Treatment)]
+  }
+
+  # Return both results and proportions
+  return(structure(list(
+    results = results_df,
+    proportions = prop_by_treatment,
+    detailed_results = results_list,
+    alternative = alternative
+  ),class="RSCABS"))
+}
+
+
+#' Printing method for run_all_threshold_tests results
+#'
+#' @param x
+#'
+#' @return printed results from run_all_threshold_tests
+#' @export
+#'
+print.RSCABS <- function(x){
+  cat("\nresult table: \n\n")
+  print(x$results)
+  cat("\n Summarized Proportions: \n")
+  print(x$proportions)
+  cat("\n Alternative Hypothesis: \n")
+  print(x$alternative)
+  invisible(x$detailed_results)
+}
+
+
+
+
+
+#' Perform Step-Down Rao-Scott Adjusted Cochran-Armitage Trend Test Procedure
+#'
+#' This function performs a step-down procedure using Rao-Scott adjusted Cochran-Armitage
+#' trend tests (RSCABS). The procedure systematically excludes the highest treatment groups
+#' to identify at which treatment level the dose-response relationship becomes significant.
+#'
+#' @param data A data frame containing fish injury data
+#' @param treatment_col Name of the column containing treatment groups (default: "tmt")
+#' @param treatment_order Optional vector specifying the order of treatment groups from
+#'        lowest to highest dose. If NULL, alphabetical order is used (default: NULL)
+#' @param max_score Maximum score value to consider (default: NULL, auto-detected)
+#' @param min_score Minimum score value to consider (default: 1)
+#' @param score_cols Character vector of column names containing injury scores (default: NULL, auto-detected)
+#' @param replicate_col Name of the column containing tank/replicate IDs (default: "tank")
+#' @param total_col Name of the column containing total counts (default: "total")
+#' @param direction Character string indicating threshold direction: "greater" for ≥ threshold,
+#'        "less" for ≤ threshold (default: "greater")
+#' @param alternative Character string specifying the alternative hypothesis:
+#'        "two.sided", "greater" (proportion increases with treatment level),
+#'        or "less" (proportion decreases with treatment level) (default: "greater")
+#' @param include_fisher Logical indicating whether to use Fisher's exact test when RSCA fails (default: TRUE)
+#'
+#' @return A list of class "StepDownRSCABS" containing:
+#'   \item{combined_results}{Data frame with test results for all steps and thresholds}
+#'   \item{step_results}{List of RSCABS objects for each step}
+#'   \item{summary}{Summary of significant findings by step}
+#'   \item{lowest_significant}{Information about the lowest significant treatment level}
+#'   \item{parameters}{List of parameters used for the analysis}
+#'
+#' @examples
+#' # Example data
+#' fish_data <- data.frame(
+#'   tmt = c(rep("Control", 8), rep("Low", 4), rep("Medium", 4), rep("High", 4)),
+#'   tank = c(paste0("C", 1:8), paste0("L", 1:4), paste0("M", 1:4), paste0("H", 1:4)),
+#'   S0 = c(3, 2, 3, 2, 2, 1, 2, 3, 3, 3, 4, 2, 2, 3, 2, 2, 2, 3, 1, 2),
+#'   S1 = c(1, 2, 0, 1, 2, 2, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 2, 0),
+#'   S2 = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1),
+#'   S3 = c(0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1),
+#'   total = c(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
+#' )
+#'
+#' # Run step-down procedure with default parameters
+#' result <- step_down_RSCABS(fish_data, treatment_col = "tmt",
+#'                           treatment_order = c("Control", "Low", "Medium", "High"))
+#'
+#' # Print results
+#' print(result)
+#'
+#' # Plot results
+#' plot(result)
+#'
+#' @importFrom dplyr filter mutate select bind_rows arrange
+#' @importFrom stats aggregate
+#' @export
+step_down_RSCABS <- function(data,
+                             treatment_col = "tmt",
+                             treatment_order = NULL,
+                             max_score = NULL,
+                             min_score = 1,
+                             score_cols = NULL,
+                             replicate_col = "tank",
+                             total_col = "total",
+                             direction = "greater",
+                             alternative = "greater",
+                             include_fisher = TRUE) {
+
+  # Input validation
+  if (!treatment_col %in% names(data)) {
+    stop("Treatment column '", treatment_col, "' not found in data")
+  }
+
+  # Get treatment levels
+  if (is.null(treatment_order)) {
+    treatments <- sort(unique(data[[treatment_col]]))
+  } else {
+    # Validate treatment_order
+    if (!all(unique(data[[treatment_col]]) %in% treatment_order)) {
+      missing <- setdiff(unique(data[[treatment_col]]), treatment_order)
+      stop("Some treatments in data are not in treatment_order: ",
+           paste(missing, collapse = ", "))
+    }
+    if (!all(treatment_order %in% unique(data[[treatment_col]]))) {
+      extra <- setdiff(treatment_order, unique(data[[treatment_col]]))
+      stop("Some treatments in treatment_order are not in data: ",
+           paste(extra, collapse = ", "))
+    }
+    treatments <- treatment_order
+  }
+
+  # Ensure treatment is a factor with the correct order
+  data[[treatment_col]] <- factor(data[[treatment_col]], levels = treatments)
+
+  # Number of treatment levels
+  n_treatments <- length(treatments)
+
+  # Check if there are at least 2 treatment levels
+  if (n_treatments < 2) {
+    stop("At least 2 treatment levels are required for step-down procedure")
+  }
+
+  # Initialize list to store results
+  step_results <- list()
+  combined_results <- data.frame()
+
+  # Perform tests for each step
+  for (i in 1:(n_treatments-1)) {
+    # Select treatments for this step
+    included_treatments <- treatments[1:(n_treatments-i+1)]
+    ## need to droplevels so that the step step does not calculate for non-existing treatment levels.
+    step_data <- droplevels(data[data[[treatment_col]] %in% included_treatments, ])
+
+    # Run tests
+    results <- run_all_threshold_tests(
+      data = step_data,
+      max_score = max_score,
+      min_score = min_score,
+      score_cols = score_cols,
+      treatment_col = treatment_col,
+      replicate_col = replicate_col,
+      total_col = total_col,
+      direction = direction,
+      alternative = alternative,
+      include_fisher = include_fisher
+    )
+
+    # Add step information to results dataframe
+    results$results$Step <- i
+    results$results$Included_Treatments <- paste(included_treatments, collapse = ", ")
+    results$results$Highest_Treatment <- included_treatments[length(included_treatments)]
+
+    # Store results
+    step_results[[i]] <- results
+    combined_results <- bind_rows(combined_results, results$results)
+  }
+
+  # Create summary of significant findings by step
+  summary <- list()
+  lowest_significant <- list(
+    step = 0,
+    treatment = NA,
+    threshold = NA,
+    p_value = NA
+  )
+  ##browser()
+  for (i in 1:(n_treatments-1)) {
+    step_data <- combined_results[combined_results$Step == i, ]
+    sig_thresholds <- step_data[step_data$P_value < 0.05, "Threshold"]
+
+    summary[[i]] <- list(
+      step = i,
+      included_treatments = unique(step_data$Included_Treatments),
+      highest_treatment = unique(step_data$Highest_Treatment),
+      significant_thresholds = sig_thresholds,
+      has_significant_findings = length(sig_thresholds) > 0
+    )
+
+    # Update lowest significant treatment if this step has significant findings
+    if (length(sig_thresholds) > 0 && i > lowest_significant$step) {
+      lowest_sig_row <- step_data[step_data$P_value < 0.05, ][
+        which.min(step_data[step_data$P_value < 0.05, "P_value"]), ]
+
+      lowest_significant$step <- i
+      lowest_significant$treatment <- unique(step_data$Highest_Treatment)
+      lowest_significant$threshold <- lowest_sig_row$Threshold
+      lowest_significant$p_value <- lowest_sig_row$P_value
+    }
+  }
+
+  # If no significant findings were found, update lowest_significant
+  if (lowest_significant$step==0) {
+    lowest_significant$step <- NA
+  }
+
+  # Return results
+  result <- list(
+    combined_results = combined_results,
+    step_results = step_results,
+    summary = summary,
+    lowest_significant = lowest_significant,
+    parameters = list(
+      treatment_col = treatment_col,
+      treatment_order = treatments,
+      max_score = max_score,
+      min_score = min_score,
+      score_cols = score_cols,
+      replicate_col = replicate_col,
+      total_col = total_col,
+      direction = direction,
+      alternative = alternative,
+      include_fisher = include_fisher
+    )
+  )
+
+  class(result) <- "StepDownRSCABS"
+  return(result)
+}
+
+#' Print method for StepDownRSCABS objects
+#'
+#' @param x A StepDownRSCABS object
+#' @param ... Additional arguments (not used)
+#'
+#' @return Invisibly returns the object
+#' @export
+print.StepDownRSCABS <- function(x, printLowest=FALSE,...) {
+  cat("Step-Down RSCABS Analysis\n")
+  cat("========================\n\n")
+
+  cat("Parameters:\n")
+  cat("  Direction:", x$parameters$direction, "\n")
+  cat("  Alternative hypothesis:", x$parameters$alternative, "\n")
+  cat("  Treatment levels:", paste(x$parameters$treatment_order, collapse = ", "), "\n\n")
+
+  cat("Summary of findings:\n")
+  for (i in seq_along(x$summary)) {
+    cat("Step", i, ": ")
+    cat("Included treatments:", x$summary[[i]]$included_treatments, "\n")
+
+    if (x$summary[[i]]$has_significant_findings) {
+      cat("  Significant thresholds:",
+          paste(x$summary[[i]]$significant_thresholds, collapse = ", "), "\n")
+    } else {
+      cat("  No significant findings\n")
+    }
+  }
+  if(printLowest){
+    cat("\nLowest treatment level with significant findings:\n")
+    if (!is.na(x$lowest_significant$step)) {
+      cat("  Treatment:", x$lowest_significant$treatment, "\n")
+      cat("  Threshold:", x$lowest_significant$threshold, "\n")
+      cat("  P-value:", round(x$lowest_significant$p_value, 4), "\n")
+    } else {
+      cat("  No significant findings at any treatment level\n")
+    }
+
+  }
+
+  invisible(x)
+}
+
+#' Plot method for StepDownRSCABS objects
+#'
+#' @param x A StepDownRSCABS object
+#' @param ... Additional arguments (not used)
+#'
+#' @return A ggplot object
+#' @importFrom ggplot2 ggplot aes geom_tile geom_text scale_fill_gradient scale_y_reverse labs theme_minimal theme element_text element_blank
+#' @export
+plot.StepDownRSCABS <- function(x, ...) {
+  # Prepare data for plotting
+  plot_data <- x$combined_results
+
+  # Create step labels
+  step_labels <- sapply(seq_along(x$summary), function(i) {
+    paste("Step", i, ":", x$summary[[i]]$highest_treatment)
+  })
+
+  # Create the plot
+  p <- ggplot(plot_data, aes(x = Threshold, y = Step, fill = -log10(P_value))) +
+    geom_tile(color = "white") +
+    geom_text(aes(label = sprintf("%.4f", round(P_value, 4))), color = "black", size = 3) +
+    scale_fill_gradient(low = "white", high = "steelblue",
+                        name = "-log10(p-value)") +
+    scale_y_reverse(breaks = seq_along(step_labels),
+                    labels = step_labels) +
+    labs(
+      title = "Step-Down RSCABS Results",
+      subtitle = "P-values for each threshold at each step",
+      x = "Injury Threshold",
+      y = "Step"
+    ) +
+    theme_minimal() +
+    theme(
+      plot.title = element_text(size = 14, face = "bold"),
+      axis.text = element_text(size = 10),
+      panel.grid.major = element_blank(),
+      panel.grid.minor = element_blank(),
+      legend.position = "right"
+    )
+
+  return(p)
+}
+
+#' Summary method for StepDownRSCABS objects
+#'
+#' @param object A StepDownRSCABS object
+#' @param ... Additional arguments (not used)
+#'
+#' @return A list with summary information
+#' @export
+summary.StepDownRSCABS <- function(object, ...) {
+  # Create a summary data frame of results
+  summary_df <- data.frame(
+    Step = integer(),
+    Included_Treatments = character(),
+    Highest_Treatment = character(),
+    Significant_Thresholds = character(),
+    Min_P_Value = numeric(),
+    stringsAsFactors = FALSE
+  )
+
+  for (i in seq_along(object$summary)) {
+    step_summary <- object$summary[[i]]
+
+    if (step_summary$has_significant_findings) {
+      sig_thresholds <- paste(step_summary$significant_thresholds, collapse = ", ")
+      min_p <- min(object$combined_results$P_value[
+        object$combined_results$Step == i &
+          object$combined_results$Threshold %in% step_summary$significant_thresholds
+      ])
+    } else {
+      sig_thresholds <- "None"
+      min_p <- NA
+    }
+
+    summary_df <- rbind(summary_df, data.frame(
+      Step = i,
+      Included_Treatments = step_summary$included_treatments,
+      Highest_Treatment = step_summary$highest_treatment,
+      Significant_Thresholds = sig_thresholds,
+      Min_P_Value = min_p,
+      stringsAsFactors = FALSE
+    ))
+  }
+
+  # Return summary information
+  list(
+    summary_table = summary_df,
+    lowest_significant = object$lowest_significant,
+    parameters = object$parameters
+  )
+}
diff --git a/R/drc_Helper.R b/R/drc_Helper.R
index 9b00e00..66fd647 100644
--- a/R/drc_Helper.R
+++ b/R/drc_Helper.R
@@ -20,6 +20,14 @@
 #'
 #' @examples
 #' \dontrun{
+#' data("dat_medium")
+#' dat_medium <- dat_medium %>% mutate(Treatment=factor(Dose,levels=unique(Dose)))
+#' dat_medium$Response[dat_medium$Response < 0] <- 0
+#' mod <- drm(Response~Dose,data=dat_medium,fct=LN.4())
+#' p1 <- plot.modList(list(mod))
+#' addECxCI(p1,object=mod,EDres=NULL,trend="Decrease",endpoint="EC", respLev=c(10,20,50),
+#' textAjust.x=0.01,textAjust.y=0.3,useObsCtr=FALSE,d0=NULL,textsize = 4,lineheight = 0.5,xmin=0.012)+
+#' ylab("Response Variable [unit]") + xlab("Concentration [µg a.s./L]")
 #' }
 addECxCI <- function(p = NULL, object, EDres = NULL, trend = "Decrease", endpoint = "ErC", respLev = c(10, 20, 50),
                      textAjust.x = 0.1, textAjust.y = 0.05, useObsCtr = FALSE, d0 = NULL, textsize = 2, lineheight = 1, xmin = 0.05, ...) {
@@ -364,7 +372,6 @@ mselect.ZG <- mselect.plus
 #'
 #' @return The function prints the `Comparison` and `EFSA` components
 #'  of the object by default.
-#' @S3method
 #' @export
 print.drcComp <- function(x, ...) {
   x <- x[c("Comparison", "EFSA")]
@@ -388,7 +395,7 @@ print.drcComp <- function(x, ...) {
 #' user input when the use does not want to use the
 #' data provided in the modList
 #'
-#' @return
+#' @return a ggplot object
 #' @export plot.modList
 #'
 #' @examples
@@ -582,15 +589,14 @@ calcSteepnessOverlap <- function(mod = NULL, obj = NULL, trend = "Decrease", CI
 #' ECx calculation together with normalized width proposed by EFSA SO.
 #'
 #' @param modList list of models
-#' @param respLev
-#' @param trend
-#' @param ...
+#' @param respLev response levels or ECx for example, 10, 20, 50
+#' @param trend "Decrease" or "Incrrease".
+#' @param ... optional to pass to ED function.
 #'
 #' @return ED result table
 #' @export mselect.ED
 #' @importFrom purrr map2
 #'
-#' @examples
 mselect.ED <- function(modList, respLev = c(10, 20, 50),
                        trend = "Decrease", ...) {
   edList <- lapply(modList, function(obj) {
@@ -618,16 +624,49 @@ mselect.ED <- function(modList, respLev = c(10, 20, 50),
   return(edResTab)
 }
 
-## Normalised width
-#' Title
+
+#' Calculate Normalized Width for Effect Concentration Estimates
 #'
-#' @param x
-#' @param ED
+#' @description
+#' Calculates the normalized width of confidence intervals for effect concentration
+#' (EC) estimates by dividing the interval width by the estimate value. Also
+#' provides a rating for the calculated normalized width.
+#'
+#' @param x A data frame containing EC estimates, typically output from ED
+#'   function. Must contain columns 'Upper', 'Lower', and 'Estimate'.
+#' @param ED Character string indicating the source of EC estimates. Must be either
+#'   "ZG" (default) or "drc".
+#'
+#' @return A data frame with two columns:
+#'   \itemize{
+#'     \item NW: Normalized width calculated as (Upper - Lower) / Estimate
+#'     \item Rating: Categorical rating of the normalized width
+#'   }
+#'   Row names are preserved from input for "ZG" or modified to "EC_x" format
+#'   for "drc".
 #'
-#' @return
 #' @export
 #'
 #' @examples
+#' # Example with ZG format
+#' ec_data <- data.frame(
+#'   Lower = c(1.2, 2.3),
+#'   Upper = c(1.8, 3.1),
+#'   Estimate = c(1.5, 2.7),
+#'   row.names = c("EC10", "EC50")
+#' )
+#' calcNW(ec_data, ED = "ZG")
+#'
+#' # Example with drc format
+#' drc_data <- data.frame(
+#'   Lower = c(1.2, 2.3),
+#'   Upper = c(1.8, 3.1),
+#'   Estimate = c(1.5, 2.7),
+#'   row.names = c("Dose:1:10", "Dose:1:50")
+#' )
+#' calcNW(drc_data, ED = "drc")
+#'
+#' @seealso \code{\link{ECx_rating}} for the rating system used
 calcNW <- function(x, ED = "ZG") {
   x <- as.data.frame(x) ## x is the ED object from ED function.
   out <- as.data.frame((x$Upper - x$Lower) / x$Estimate)
@@ -640,14 +679,33 @@ calcNW <- function(x, ED = "ZG") {
   return(out)
 }
 
-#' Title
+#' Rate Effect Concentration Estimates Based on Normalized Width
+#'
+#' @description
+#' Assigns qualitative ratings to normalized width values of confidence intervals
+#' based on predefined thresholds.
 #'
-#' @param x
+#' @param x Numeric vector of normalized width values to be rated
+#'
+#' @return Character vector of ratings with the following categories:
+#'   \itemize{
+#'     \item "Excellent" for NW < 0.2
+#'     \item "Good" for 0.2 <= NW < 0.5
+#'     \item "Fair" for 0.5 <= NW < 1
+#'     \item "Poor" for 1 <= NW < 2
+#'     \item "Bad" for NW >= 2
+#'     \item "Not defined" for NA or NaN values
+#'   }
 #'
-#' @return
 #' @export
 #'
 #' @examples
+#' # Basic examples
+#' ECx_rating(c(0.1, 0.3, 0.7, 1.5, 2.5, NA))
+#'
+#' # Example with typical normalized widths
+#' nw_values <- c(0.15, 0.45, 0.95, 1.8, 2.5)
+#' ECx_rating(nw_values)
 ECx_rating <- function(x) {
   ## normalised width (NW) of CI
   ## NW   |  Rating # nolint: commented_code_linter.
@@ -680,17 +738,45 @@ ECx_rating <- function(x) {
   return(outRating)
 }
 
-#' Compare models with multiple criteria
-#'
-#' @param modRes
-#' @param modList
-#' @param trend
-#' @param ...
+#' Compare Dose-Response Models Using Multiple Criteria
+#'
+#' @description
+#' Performs comprehensive comparison of dose-response models using multiple
+#' criteria including model fit, certainty of protection, and steepness.
+#'
+#' @param modRes Optional. Result object containing model comparison information.
+#'   If provided, modList parameter is ignored.
+#' @param modList Optional. List of fitted dose-response models to compare.
+#'   Required if modRes is NULL.
+#' @param trend Character string specifying the expected trend direction.
+#'   Must be "Decrease" or "Increase".
+#' @param CI Character string specifying the confidence interval method.
+#'   One of "delta", "inv", or "bmd-inv". Defaults to "delta".
+#' @param ... Additional arguments passed to internal functions
+#'
+#' @return A data frame containing model comparison results with columns:
+#'   \itemize{
+#'     \item Original comparison metrics
+#'     \item Certainty_Protection: Measure of protection certainty
+#'     \item Steepness: Model steepness evaluation
+#'     \item No Effect p-val: P-value for test against no-effect model
+#'   }
 #'
-#' @return
 #' @export
 #'
 #' @examples
+#' \dontrun{
+#' data("dat_medium")
+#' dat_medium <- dat_medium %>% mutate(Treatment=factor(Dose,levels=unique(Dose)))
+#' dat_medium$Response[dat_medium$Response < 0] <- 0
+#' mod <- drm(Response~Dose,data=dat_medium,fct=LL.3())
+#' fctList <- list(LN.4(),LL.4(),W1.3(),LL2.2())
+#' res <- mselect.plus(mod,fctList = fctList )
+#' modList <- res$modList
+#'
+#' # Compare models using existing comparison results
+#' drcCompare(modRes = res, trend = "Decrease")
+#' }
 drcCompare <- function(modRes = NULL, modList = NULL, trend = "Decrease",
                        CI = c("delta", "inv", "bmd-inv"), ...) { # nolint
   CI <- match.arg(CI)
@@ -699,7 +785,10 @@ drcCompare <- function(modRes = NULL, modList = NULL, trend = "Decrease",
     modList <- modRes$modList
     comparison <- modRes$Comparison
   } else {
-    if (is.null(modList)) stop("Need the model output list from previous step!")
+    if (is.null(modList)) stop("Need the model output list from previous step!") else{
+      ## use modList directly
+      stop("A direct model comparison has not yet implemented, this is a placeholder.")
+    }
   }
   # A significance test is provided for the comparison of the dose-response model considered and the simple linear regression model with slope 0 (a horizontal regression line corresponding to no dose effect) # nolint: line_length_linter.
   # The likelihood ratio test statistic and the corresponding degrees of freedom and p-value are reported. # nolint: line_length_linter.
diff --git a/README.Rmd b/README.Rmd
index dba6ca8..b8fdcc4 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -1,5 +1,7 @@
 ---
 output: github_document
+editor_options: 
+  chunk_output_type: console
 ---
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
@@ -112,10 +114,33 @@ p
 ## Adding ECx and ECx CI's to the plots
 
 ```{r}
+p1 <- plot.modList(modList[1])
+addECxCI(p1,object=modList[[1]],EDres=NULL,trend="Decrease",endpoint="EC", respLev=c(10,20,50),
+                     textAjust.x=0.01,textAjust.y=0.3,useObsCtr=FALSE,d0=NULL,textsize = 4,lineheight = 0.5,xmin=0.012)+ ylab("Response Variable [unit]") + xlab("Concentration [µg a.s./L]")
 ## addECxCI(p)
 ```
 
+## Report ECx
 
+```{r}
+resED <- t(edResTab[1:3, c(2,4,5,6)])
+colnames(resED) <- paste("EC", c(10,20,50))
+knitr::kable(resED,caption = "Response Variable at day N",digits = 3)
+```
+**Calculate specific ECx:  **
+```{r}
+mod <-modList[[1]]
+edres <- ED.plus(mod,c(5,10,20,50),trend="Decrease")
+edres%>%knitr::kable(.,digits = 3)
+```
+
+
+## Model Output
+
+```{r}
+modsum <- summary(mod)
+knitr::kable(coef(modsum),digits = 3)
+```
 
 
 ## ToDo
diff --git a/README.md b/README.md
index cd3e15e..2dc77c5 100644
--- a/README.md
+++ b/README.md
@@ -158,10 +158,64 @@ p
 
 ## Adding ECx and ECx CI’s to the plots
 
+``` r
+p1 <- plot.modList(modList[1])
+addECxCI(p1,object=modList[[1]],EDres=NULL,trend="Decrease",endpoint="EC", respLev=c(10,20,50),
+                     textAjust.x=0.01,textAjust.y=0.3,useObsCtr=FALSE,d0=NULL,textsize = 4,lineheight = 0.5,xmin=0.012)+ ylab("Response Variable [unit]") + xlab("Concentration [µg a.s./L]")
+```
+
+<img src="man/figures/README-unnamed-chunk-7-1.png" width="100%" />
+
 ``` r
 ## addECxCI(p)
 ```
 
+## Report ECx
+
+``` r
+resED <- t(edResTab[1:3, c(2,4,5,6)])
+colnames(resED) <- paste("EC", c(10,20,50))
+knitr::kable(resED,caption = "Response Variable at day N",digits = 3)
+```
+
+|          | EC 10 | EC 20 | EC 50 |
+|:---------|------:|------:|------:|
+| Estimate | 1.699 | 2.067 | 3.034 |
+| Lower    | 1.465 | 1.817 | 2.786 |
+| Upper    | 1.990 | 2.321 | 3.284 |
+| NW       | 0.309 | 0.244 | 0.164 |
+
+Response Variable at day N
+
+**Calculate specific ECx: **
+
+``` r
+mod <-modList[[1]]
+edres <- ED.plus(mod,c(5,10,20,50),trend="Decrease")
+edres%>%knitr::kable(.,digits = 3)
+```
+
+|       | Estimate | Std. Error | Lower | Upper |
+|:------|---------:|-----------:|------:|------:|
+| EC 5  |    1.447 |      0.163 | 1.107 | 1.787 |
+| EC 10 |    1.699 |      0.159 | 1.367 | 2.032 |
+| EC 20 |    2.067 |      0.151 | 1.753 | 2.382 |
+| EC 50 |    3.034 |      0.152 | 2.716 | 3.352 |
+
+## Model Output
+
+``` r
+modsum <- summary(mod)
+knitr::kable(coef(modsum),digits = 3)
+```
+
+|               | Estimate | Std. Error | t-value | p-value |
+|:--------------|---------:|-----------:|--------:|--------:|
+| b:(Intercept) |   -2.300 |      0.309 |  -7.441 |   0.000 |
+| c:(Intercept) |    0.532 |      0.177 |   3.005 |   0.007 |
+| d:(Intercept) |    7.719 |      0.174 |  44.474 |   0.000 |
+| e:(Intercept) |    2.914 |      0.148 |  19.750 |   0.000 |
+
 ## ToDo
 
 - [ ] Develop all test cases for NOEC functions
diff --git a/inst/DevNotes/function_dev/plot_oridinal.R b/inst/DevNotes/function_dev/plot_oridinal.R
new file mode 100644
index 0000000..443250e
--- /dev/null
+++ b/inst/DevNotes/function_dev/plot_oridinal.R
@@ -0,0 +1,693 @@
+# Load necessary libraries
+library(ggplot2)
+library(dplyr)
+library(tidyr)
+library(patchwork)  # For combining plots
+
+
+
+# Load necessary libraries
+library(ggplot2)
+library(dplyr)
+library(tidyr)
+library(patchwork)  # For combining plots
+
+# Create the simulated fish data
+simulated_fish_data <- tibble::tibble(
+  treatment = c(rep("Control", 8), rep("Low", 4), rep("Medium", 4), rep("High", 4)),
+  tank = c(paste0("C", 1:8), paste0("L", 1:4), paste0("M", 1:4), paste0("H", 1:4)),
+  S0 = c(3, 2, 4, 3, 2, 1, 2, 3, 3, 3, 4, 2, 2, 3, 2, 2, 2, 3, 1, 2),
+  S1 = c(1, 2, 0, 1, 2, 2, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 2, 0),
+  S2 = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1),
+  S3 = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1),
+  total = c(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
+)
+
+# Ensure treatment is an ordered factor for proper plotting
+simulated_fish_data$treatment <- factor(simulated_fish_data$treatment,
+                                        levels = c("Control", "Low", "Medium", "High"))
+
+# Convert to long format for easier plotting
+fish_data_long <- simulated_fish_data %>%
+  pivot_longer(
+    cols = starts_with("S"),
+    names_to = "score",
+    values_to = "count"
+  ) %>%
+  # Ensure score is an ordered factor
+  mutate(score = factor(score, levels = c("S0", "S1", "S2", "S3")))
+
+# 1. Stacked Bar Plot by Treatment Group
+p1 <- fish_data_long %>%
+  group_by(treatment, score) %>%
+  summarize(total_count = sum(count), .groups = "drop") %>%
+  ggplot(aes(x = treatment, y = total_count, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  labs(
+    title = "Distribution of Injury Scores by Treatment Group",
+    x = "Treatment Group",
+    y = "Number of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "right"
+  )
+
+# 2. Proportion Plot by Treatment Group
+p2 <- fish_data_long %>%
+  group_by(treatment, score) %>%
+  summarize(total_count = sum(count), .groups = "drop") %>%
+  group_by(treatment) %>%
+  mutate(proportion = total_count / sum(total_count)) %>%
+  ggplot(aes(x = treatment, y = proportion, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  labs(
+    title = "Proportion of Injury Scores by Treatment Group",
+    x = "Treatment Group",
+    y = "Proportion of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "right"
+  )
+
+# 3. Boxplot of Injury Proportions by Treatment
+p3 <- simulated_fish_data %>%
+  mutate(
+    injury_prop = (S1 + S2 + S3) / total,
+    severe_injury_prop = (S2 + S3) / total
+  ) %>%
+  pivot_longer(
+    cols = c(injury_prop, severe_injury_prop),
+    names_to = "injury_type",
+    values_to = "proportion"
+  ) %>%
+  mutate(injury_type = factor(injury_type,
+                              levels = c("injury_prop", "severe_injury_prop"),
+                              labels = c("Any Injury (S1-S3)", "Severe Injury (S2-S3)"))) %>%
+  ggplot(aes(x = treatment, y = proportion, fill = treatment)) +
+  geom_boxplot(alpha = 0.7) +
+  geom_jitter(width = 0.2, height = 0, alpha = 0.7) +
+  facet_wrap(~ injury_type) +
+  labs(
+    title = "Distribution of Injury Proportions by Treatment Group",
+    x = "Treatment Group",
+    y = "Proportion of Fish",
+    fill = "Treatment"
+  ) +
+  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
+  scale_fill_brewer(palette = "Set2") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "none"
+  )
+
+# 4. Heat Map of Injury Counts by Tank
+p4 <- fish_data_long %>%
+  ggplot(aes(x = tank, y = score, fill = count)) +
+  geom_tile(color = "white") +
+  geom_text(aes(label = count), color = "black") +
+  facet_grid(. ~ treatment, scales = "free_x", space = "free_x") +
+  labs(
+    title = "Heat Map of Injury Counts by Tank",
+    x = "Tank",
+    y = "Injury Score",
+    fill = "Count"
+  ) +
+  scale_fill_gradient(low = "white", high = "#fc8d62") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    axis.text.x = element_text(angle = 45, hjust = 1),
+    panel.grid = element_blank(),
+    legend.position = "right"
+  )
+
+# 5. Dot plot showing proportion of each injury type by treatment
+p5 <- fish_data_long %>%
+  group_by(treatment, score) %>%
+  summarize(
+    total_count = sum(count),
+    .groups = "drop"
+  ) %>%
+  group_by(treatment) %>%
+  mutate(proportion = total_count / sum(total_count)) %>%
+  ggplot(aes(x = treatment, y = proportion, color = score)) +
+  geom_point(size = 4) +
+  geom_line(aes(group = score), linewidth = 1) +
+  facet_wrap(~ score, nrow = 1) +
+  labs(
+    title = "Trend in Proportion of Each Injury Score Across Treatments",
+    x = "Treatment Group",
+    y = "Proportion of Fish",
+    color = "Injury Score"
+  ) +
+  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
+  scale_color_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "none",
+    panel.spacing = unit(1, "lines")
+  )
+
+# 6. Stacked bar plot for each tank
+p6 <- fish_data_long %>%
+  ggplot(aes(x = tank, y = count, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  facet_grid(. ~ treatment, scales = "free_x", space = "free_x") +
+  labs(
+    title = "Distribution of Injury Scores by Individual Tank",
+    x = "Tank",
+    y = "Number of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    axis.text.x = element_text(angle = 45, hjust = 1),
+    legend.position = "right"
+  )
+
+# Combine plots using patchwork
+top_row <- p1 + p2
+middle_row <- p3
+bottom_row <- p4 / p5
+last_row <- p6
+
+combined_plot <- (top_row) / (middle_row) / (bottom_row) / (last_row) +
+  plot_layout(heights = c(1, 1, 2, 1)) +
+  plot_annotation(
+    title = "Visualization of Simulated Fish Injury Data",
+    subtitle = "Multiple visualizations showing the distribution of injury scores across treatment groups and tanks",
+    theme = theme(
+      plot.title = element_text(size = 16, face = "bold"),
+      plot.subtitle = element_text(size = 12)
+    )
+  )
+
+# Print the combined plot
+print(combined_plot)
+
+# Statistical summary table
+summary_table <- simulated_fish_data %>%
+  group_by(treatment) %>%
+  summarize(
+    n_tanks = n(),
+    total_fish = sum(total),
+    s0_count = sum(S0),
+    s0_percent = 100 * s0_count / total_fish,
+    s1_count = sum(S1),
+    s1_percent = 100 * s1_count / total_fish,
+    s2_count = sum(S2),
+    s2_percent = 100 * s2_count / total_fish,
+    s3_count = sum(S3),
+    s3_percent = 100 * s3_count / total_fish,
+    any_injury = sum(S1 + S2 + S3),
+    any_injury_percent = 100 * any_injury / total_fish,
+    .groups = "drop"
+  ) %>%
+  mutate(across(ends_with("percent"), ~ round(., 1)))
+
+# Print the summary table
+print(summary_table)
+
+# Calculate injury severity index (weighted average of injury scores)
+severity_by_treatment <- simulated_fish_data %>%
+  mutate(
+    severity_index = (0*S0 + 1*S1 + 2*S2 + 3*S3) / total
+  ) %>%
+  group_by(treatment) %>%
+  summarize(
+    mean_severity = mean(severity_index),
+    sd_severity = sd(severity_index),
+    .groups = "drop"
+  )
+
+# Print severity index
+print(severity_by_treatment)
+
+
+
+# Example data with increasing trend in injury rates
+fish_data <- data.frame(
+  tmt = rep(c("Control", "Low", "Medium", "High"), each = 3),
+  tank = paste0(rep(c("Control", "Low", "Medium", "High"), each = 3),
+                rep(1:3, times = 4)),
+  S0 = c(9,8,9, 7,6,7, 5,4,5, 3,2,3),
+  S1 = c(1,2,1, 3,4,3, 5,6,5, 7,8,7),
+  total = rep(10, 12)
+)
+
+# Ensure treatment is an ordered factor for proper plotting
+fish_data$tmt <- factor(fish_data$tmt, levels = c("Control", "Low", "Medium", "High"))
+
+# Convert to long format for easier plotting
+fish_data_long <- fish_data %>%
+  pivot_longer(
+    cols = starts_with("S"),
+    names_to = "score",
+    values_to = "count"
+  ) %>%
+  # Ensure score is an ordered factor
+  mutate(score = factor(score, levels = paste0("S", 0:1)))
+
+# 1. Stacked Bar Plot by Treatment Group
+p1 <- ggplot(fish_data_long, aes(x = tmt, y = count, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  labs(
+    title = "Distribution of Injury Scores by Treatment Group",
+    x = "Treatment Group",
+    y = "Number of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_manual(values = c("S0" = "#66c2a5", "S1" = "#fc8d62")) +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "bottom"
+  )
+
+# 2. Proportion Plot by Treatment Group
+p2 <- fish_data_long %>%
+  group_by(tmt, score) %>%
+  summarize(total_count = sum(count), .groups = "drop") %>%
+  group_by(tmt) %>%
+  mutate(proportion = total_count / sum(total_count)) %>%
+  ggplot(aes(x = tmt, y = proportion, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  labs(
+    title = "Proportion of Injury Scores by Treatment Group",
+    x = "Treatment Group",
+    y = "Proportion of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_manual(values = c("S0" = "#66c2a5", "S1" = "#fc8d62")) +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "bottom"
+  )
+
+# 3. Individual Tank Data with Treatment Averages
+p3 <- fish_data %>%
+  mutate(S1_prop = S1 / total) %>%
+  ggplot(aes(x = tmt, y = S1_prop)) +
+  # Add individual tank points
+  geom_point(aes(color = tank), position = position_jitter(width = 0.2), size = 3, alpha = 0.7) +
+  # Add treatment group means
+  stat_summary(fun = mean, geom = "point", shape = 18, size = 5, color = "black") +
+  stat_summary(fun = mean, geom = "line", group = 1, linewidth = 1, color = "black") +
+  labs(
+    title = "Proportion of S1 Injuries by Treatment Group",
+    subtitle = "Individual tanks shown as points, group means as diamonds",
+    x = "Treatment Group",
+    y = "Proportion of S1 Injuries",
+    color = "Tank"
+  ) +
+  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "bottom"
+  )
+
+# 4. Boxplot of S1 Proportions
+p4 <- fish_data %>%
+  mutate(S1_prop = S1 / total) %>%
+  ggplot(aes(x = tmt, y = S1_prop)) +
+  geom_boxplot(aes(fill = tmt), alpha = 0.7) +
+  geom_point(position = position_jitter(width = 0.2), size = 2) +
+  labs(
+    title = "Boxplot of S1 Injury Proportions by Treatment",
+    x = "Treatment Group",
+    y = "Proportion of S1 Injuries"
+  ) +
+  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
+  scale_fill_brewer(palette = "Set2") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "none"
+  )
+
+# 5. Heat Map of Injury Counts by Tank
+p5 <- fish_data_long %>%
+  ggplot(aes(x = tank, y = score, fill = count)) +
+  geom_tile(color = "white") +
+  geom_text(aes(label = count), color = "black") +
+  facet_grid(. ~ tmt, scales = "free_x", space = "free_x") +
+  labs(
+    title = "Heat Map of Injury Counts by Tank",
+    x = "Tank",
+    y = "Injury Score",
+    fill = "Count"
+  ) +
+  scale_fill_gradient(low = "white", high = "#fc8d62") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    axis.text.x = element_text(angle = 45, hjust = 1),
+    panel.grid = element_blank(),
+    legend.position = "bottom"
+  )
+
+# 6. Line Plot of S1+ Proportion Trend
+p6 <- fish_data %>%
+  group_by(tmt) %>%
+  summarize(
+    total_fish = sum(total),
+    total_s1 = sum(S1),
+    prop_s1 = total_s1 / total_fish,
+    .groups = "drop"
+  ) %>%
+  ggplot(aes(x = tmt, y = prop_s1, group = 1)) +
+  geom_line(linewidth = 1.5, color = "#fc8d62") +
+  geom_point(size = 4, color = "#fc8d62") +
+  labs(
+    title = "Trend in S1+ Injury Proportion",
+    x = "Treatment Group",
+    y = "Proportion of Fish with S1+ Injuries"
+  ) +
+  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12)
+  )
+
+# Combine plots using patchwork
+top_row <- p1 + p2
+middle_row <- p3 + p4
+bottom_row <- p5 + p6
+
+combined_plot <- top_row / middle_row / bottom_row +
+  plot_layout(heights = c(1, 1, 1.2)) +
+  plot_annotation(
+    title = "Visualization of Fish Injury Data Across Treatment Groups",
+    subtitle = "Multiple visualization approaches showing the increasing trend in injury rates",
+    theme = theme(
+      plot.title = element_text(size = 16, face = "bold"),
+      plot.subtitle = element_text(size = 12)
+    )
+  )
+
+# Print the combined plot
+combined_plot
+
+# Save the plot if needed
+# ggsave("fish_injury_visualizations.png", combined_plot, width = 12, height = 14, dpi = 300)
+
+# Statistical summary table
+summary_table <- fish_data %>%
+  group_by(tmt) %>%
+  summarize(
+    n_tanks = n(),
+    total_fish = sum(total),
+    s0_count = sum(S0),
+    s0_percent = 100 * s0_count / total_fish,
+    s1_count = sum(S1),
+    s1_percent = 100 * s1_count / total_fish,
+    .groups = "drop"
+  ) %>%
+  mutate(across(ends_with("percent"), ~ round(., 1)))
+
+# Print the summary table
+print(summary_table)
+
+# Run RSCA tests with different alternative hypotheses
+run_threshold_RSCA <- function(data, threshold = 1, alternative = "two.sided") {
+  # Simple implementation for demonstration
+  affected <- data$S1  # For S1+ threshold
+
+  # Calculate Z-statistic (simplified for demonstration)
+  # This is a simplified calculation and not the full RSCA implementation
+  tmt_numeric <- as.numeric(factor(data$tmt))
+  prop_affected <- affected / data$total
+
+  # Fit linear model for demonstration
+  model <- lm(prop_affected ~ tmt_numeric)
+  z_stat <- summary(model)$coefficients[2, 3]
+
+  # Calculate p-value based on alternative
+  p_value <- switch(alternative,
+                    "two.sided" = 2 * pnorm(-abs(z_stat)),
+                    "greater" = pnorm(-z_stat),
+                    "less" = pnorm(z_stat))
+
+  return(list(Z = z_stat, p_value = p_value, alternative = alternative))
+}
+
+# Run tests with different alternatives
+two_sided_result <- run_threshold_RSCA(fish_data, alternative = "two.sided")
+greater_result <- run_threshold_RSCA(fish_data, alternative = "greater")
+less_result <- run_threshold_RSCA(fish_data, alternative = "less")
+
+# Print results
+cat("\nRSCA Test Results for S1+ Threshold:\n")
+cat("Two-sided test: Z =", round(two_sided_result$Z, 3),
+    ", p-value =", round(two_sided_result$p_value, 4), "\n")
+cat("One-sided test (greater): Z =", round(greater_result$Z, 3),
+    ", p-value =", round(greater_result$p_value, 4), "\n")
+cat("One-sided test (less): Z =", round(less_result$Z, 3),
+    ", p-value =", round(less_result$p_value, 4), "\n")
+
+
+# Load necessary libraries
+library(dplyr)
+library(knitr)
+
+# Create the simulated fish data
+simulated_fish_data <- tibble::tibble(
+  treatment = c(rep("Control", 8), rep("Low", 4), rep("Medium", 4), rep("High", 4)),
+  tank = c(paste0("C", 1:8), paste0("L", 1:4), paste0("M", 1:4), paste0("H", 1:4)),
+  S0 = c(3, 2, 3, 2, 2, 1, 2, 3, 3, 3, 4, 2, 2, 3, 2, 2, 2, 3, 1, 2),
+  S1 = c(1, 2, 0, 1, 2, 2, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 2, 0),
+  S2 = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1),
+  S3 = c(0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1),
+  total = c(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
+)
+
+# Ensure treatment is an ordered factor for proper filtering
+simulated_fish_data$treatment <- factor(simulated_fish_data$treatment,
+                                        levels = c("Control", "Low", "Medium", "High"))
+
+# Define a function to implement the run_all_threshold_tests function
+# This is a simplified implementation for demonstration
+run_all_threshold_tests <- function(data, min_score = 1, max_score = 3,
+                                    alternative = "greater", direction = "greater") {
+  # Find score columns
+  score_cols <- grep("^S[0-9]$", names(data), value = TRUE)
+
+  # Initialize results dataframe
+  results <- data.frame(
+    Threshold = character(),
+    Z_statistic = numeric(),
+    P_value = numeric(),
+    stringsAsFactors = FALSE
+  )
+
+  # For each threshold
+  for (threshold in min_score:max_score) {
+    # Determine which scores to include based on direction
+    if (direction == "greater") {
+      selected_cols <- score_cols[as.numeric(gsub("S", "", score_cols)) >= threshold]
+      threshold_label <- paste0("S", threshold, "+")
+    } else {
+      selected_cols <- score_cols[as.numeric(gsub("S", "", score_cols)) <= threshold]
+      threshold_label <- paste0("S≤", threshold)
+    }
+
+    # Calculate affected counts for each tank
+    affected <- rowSums(data[, selected_cols, drop = FALSE])
+    total <- data$total
+
+    # Calculate proportions by treatment
+    props <- aggregate(
+      cbind(affected = affected, total = total) ~ treatment,
+      data = data,
+      FUN = sum
+    )
+    props$proportion <- props$affected / props$total
+
+    # Simple trend test (linear model on proportions)
+    # This is a simplified calculation, not the full RSCA implementation
+    trt_numeric <- as.numeric(props$treatment)
+    model <- lm(proportion ~ trt_numeric, data = props)
+    z_stat <- summary(model)$coefficients[2, 3]  # t-statistic as proxy for Z
+
+    # Calculate p-value based on alternative
+    p_value <- switch(alternative,
+                      "two.sided" = 2 * pt(-abs(z_stat), df = nrow(props) - 2),
+                      "greater" = pt(-z_stat, df = nrow(props) - 2),
+                      "less" = pt(z_stat, df = nrow(props) - 2))
+
+    # Add to results
+    results <- rbind(results, data.frame(
+      Threshold = threshold_label,
+      Z_statistic = z_stat,
+      P_value = p_value,
+      stringsAsFactors = FALSE
+    ))
+  }
+
+  return(results)
+}
+
+# Perform step-down procedure
+perform_step_down <- function(data, min_score = 1, max_score = 3,
+                              alternative = "greater", direction = "greater") {
+  # Get ordered treatment levels
+  treatments <- levels(data$treatment)
+  n_treatments <- length(treatments)
+
+  # Initialize list to store results
+  step_results <- list()
+
+  # Perform tests for each step
+  for (i in 1:(n_treatments-1)) {
+    # Select treatments for this step
+    included_treatments <- treatments[1:(n_treatments-i+1)]
+    step_data <- data %>% filter(treatment %in% included_treatments)
+
+    # Run tests
+    results <- run_all_threshold_tests(
+      step_data,
+      min_score = min_score,
+      max_score = max_score,
+      alternative = alternative,
+      direction = direction
+    )
+
+    # Add step information
+    results$Step <- i
+    results$Included_Treatments <- paste(included_treatments, collapse = ", ")
+    results$Highest_Treatment <- included_treatments[length(included_treatments)]
+
+    # Store results
+    step_results[[i]] <- results
+  }
+
+  # Combine all results
+  combined_results <- do.call(rbind, step_results)
+
+  return(combined_results)
+}
+
+# Run the step-down procedure
+step_down_results <- perform_step_down(
+  simulated_fish_data,
+  min_score = 1,
+  max_score = 3,
+  alternative = "greater",
+  direction = "greater"
+)
+
+# Format the results for better display
+formatted_results <- step_down_results %>%
+  mutate(
+    Z_statistic = round(Z_statistic, 3),
+    P_value = round(P_value, 4),
+    Significance = case_when(
+      P_value < 0.001 ~ "***",
+      P_value < 0.01 ~ "**",
+      P_value < 0.05 ~ "*",
+      P_value < 0.1 ~ ".",
+      TRUE ~ ""
+    )
+  ) %>%
+  select(Step, Highest_Treatment, Included_Treatments, Threshold, Z_statistic, P_value, Significance)
+
+# Print formatted results
+cat("## Step-Down Procedure Results\n\n")
+print(kable(formatted_results, caption = "RSCA Test Results by Step and Threshold"))
+
+# Create a visualization of the step-down results
+library(ggplot2)
+
+# Plot p-values by step and threshold
+ggplot(step_down_results, aes(x = Threshold, y = Step, fill = -log10(P_value))) +
+  geom_tile(color = "white") +
+  geom_text(aes(label = sprintf("%.4f", round(P_value, 4))), color = "black", size = 3) +
+  scale_fill_gradient(low = "white", high = "steelblue",
+                      name = "-log10(p-value)") +
+  scale_y_reverse(breaks = 1:3,
+                  labels = c("All treatments",
+                             "Excluding High",
+                             "Control vs Low only")) +
+  labs(
+    title = "Step-Down Procedure Results",
+    subtitle = "P-values for each threshold at each step",
+    x = "Injury Threshold",
+    y = "Step"
+  ) +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 14, face = "bold"),
+    axis.text = element_text(size = 10),
+    panel.grid.major = element_blank(),
+    panel.grid.minor = element_blank(),
+    legend.position = "right"
+  )
+
+# Create a summary of the step-down procedure
+cat("\n\n## Step-Down Procedure Summary\n\n")
+cat("The step-down procedure systematically excludes the highest treatment groups to identify at which level the dose-response relationship becomes significant.\n\n")
+
+cat("### Step 1: All Treatment Groups\n")
+step1 <- step_down_results %>% filter(Step == 1)
+cat("Included treatments: ", unique(step1$Included_Treatments), "\n")
+sig_thresholds1 <- step1 %>% filter(P_value < 0.05) %>% pull(Threshold)
+if(length(sig_thresholds1) > 0) {
+  cat("Significant thresholds: ", paste(sig_thresholds1, collapse = ", "), "\n\n")
+} else {
+  cat("No significant thresholds found.\n\n")
+}
+
+cat("### Step 2: Excluding Highest Treatment Group\n")
+step2 <- step_down_results %>% filter(Step == 2)
+cat("Included treatments: ", unique(step2$Included_Treatments), "\n")
+sig_thresholds2 <- step2 %>% filter(P_value < 0.05) %>% pull(Threshold)
+if(length(sig_thresholds2) > 0) {
+  cat("Significant thresholds: ", paste(sig_thresholds2, collapse = ", "), "\n\n")
+} else {
+  cat("No significant thresholds found.\n\n")
+}
+
+cat("### Step 3: Control vs. Lowest Treatment Group Only\n")
+step3 <- step_down_results %>% filter(Step == 3)
+cat("Included treatments: ", unique(step3$Included_Treatments), "\n")
+sig_thresholds3 <- step3 %>% filter(P_value < 0.05) %>% pull(Threshold)
+if(length(sig_thresholds3) > 0) {
+  cat("Significant thresholds: ", paste(sig_thresholds3, collapse = ", "), "\n\n")
+} else {
+  cat("No significant thresholds found.\n\n")
+}
+
+cat("### Overall Conclusion\n")
+lowest_sig_step <- min(step_down_results %>% filter(P_value < 0.05) %>% pull(Step), Inf)
+if(is.infinite(lowest_sig_step)) {
+  cat("No significant dose-response relationship was detected at any step.\n")
+} else {
+  lowest_sig_treatment <- step_down_results %>%
+    filter(Step == lowest_sig_step) %>%
+    pull(Highest_Treatment) %>%
+    unique()
+  lowest_sig_threshold <- step_down_results %>%
+    filter(Step == lowest_sig_step, P_value < 0.05) %>%
+    arrange(P_value) %>%
+    pull(Threshold) %>%
+    first()
+
+  cat("The lowest treatment level at which a significant dose-response relationship was detected is: ",
+      lowest_sig_treatment, " (threshold: ", lowest_sig_threshold, ").\n", sep = "")
+}
+
+
diff --git a/inst/DevNotes/function_dev/test_RSCABS.R b/inst/DevNotes/function_dev/test_RSCABS.R
new file mode 100644
index 0000000..187ae18
--- /dev/null
+++ b/inst/DevNotes/function_dev/test_RSCABS.R
@@ -0,0 +1,19 @@
+dat_bcs1 <- dat_bcs1%>% dplyr::select(-c(S0_1,S0_2)) %>% mutate(tmt = ifelse(tmt == "SC", "C", tmt))
+
+res1 <- step_down_RSCABS(
+  dat_bcs1,
+  treatment_col = "tmt",replicate_col = "tank",
+  treatment_order = c("C","T1","T2","T3"),
+  alternative = "greater"
+)
+res1$combined_results
+res1
+plot(res1)
+print(res1,printLowest = T)
+
+dat1 <- convert_fish_data(dat_bcs1,direction="to_individual",treatment_col = "tmt",replicate_col = "tank" )
+dat1 <- (dat1) %>% mutate(score1=as.numeric(factor(score))-1,score2=score1+5) %>% dplyr::select(-c(score)) %>% as.data.frame
+
+## Note runRSCABS expects more than one endpoints.
+res <- runRSCABS(dat1,'tmt','tank',test.type='RS')
+res[1:3,] %>% mutate(Effect = gsub("score1","S",Effect) )
diff --git a/man/ECx_rating.Rd b/man/ECx_rating.Rd
index ec7fe4c..9072d00 100644
--- a/man/ECx_rating.Rd
+++ b/man/ECx_rating.Rd
@@ -2,13 +2,33 @@
 % Please edit documentation in R/drc_Helper.R
 \name{ECx_rating}
 \alias{ECx_rating}
-\title{Title}
+\title{Rate Effect Concentration Estimates Based on Normalized Width}
 \usage{
 ECx_rating(x)
 }
 \arguments{
-\item{x}{}
+\item{x}{Numeric vector of normalized width values to be rated}
+}
+\value{
+Character vector of ratings with the following categories:
+\itemize{
+\item "Excellent" for NW < 0.2
+\item "Good" for 0.2 <= NW < 0.5
+\item "Fair" for 0.5 <= NW < 1
+\item "Poor" for 1 <= NW < 2
+\item "Bad" for NW >= 2
+\item "Not defined" for NA or NaN values
+}
 }
 \description{
-Title
+Assigns qualitative ratings to normalized width values of confidence intervals
+based on predefined thresholds.
+}
+\examples{
+# Basic examples
+ECx_rating(c(0.1, 0.3, 0.7, 1.5, 2.5, NA))
+
+# Example with typical normalized widths
+nw_values <- c(0.15, 0.45, 0.95, 1.8, 2.5)
+ECx_rating(nw_values)
 }
diff --git a/man/addECxCI.Rd b/man/addECxCI.Rd
index 4ef67c8..0fc8cb4 100644
--- a/man/addECxCI.Rd
+++ b/man/addECxCI.Rd
@@ -58,5 +58,13 @@ adding ECx estimation and interval to the model output plot
 }
 \examples{
 \dontrun{
+data("dat_medium")
+dat_medium <- dat_medium \%>\% mutate(Treatment=factor(Dose,levels=unique(Dose)))
+dat_medium$Response[dat_medium$Response < 0] <- 0
+mod <- drm(Response~Dose,data=dat_medium,fct=LN.4())
+p1 <- plot.modList(list(mod))
+addECxCI(p1,object=mod,EDres=NULL,trend="Decrease",endpoint="EC", respLev=c(10,20,50),
+textAjust.x=0.01,textAjust.y=0.3,useObsCtr=FALSE,d0=NULL,textsize = 4,lineheight = 0.5,xmin=0.012)+
+ylab("Response Variable [unit]") + xlab("Concentration [µg a.s./L]")
 }
 }
diff --git a/man/calcNW.Rd b/man/calcNW.Rd
index d3af25a..c1fadac 100644
--- a/man/calcNW.Rd
+++ b/man/calcNW.Rd
@@ -2,13 +2,51 @@
 % Please edit documentation in R/drc_Helper.R
 \name{calcNW}
 \alias{calcNW}
-\title{Title}
+\title{Calculate Normalized Width for Effect Concentration Estimates}
 \usage{
 calcNW(x, ED = "ZG")
 }
 \arguments{
-\item{ED}{}
+\item{x}{A data frame containing EC estimates, typically output from ED
+function. Must contain columns 'Upper', 'Lower', and 'Estimate'.}
+
+\item{ED}{Character string indicating the source of EC estimates. Must be either
+"ZG" (default) or "drc".}
+}
+\value{
+A data frame with two columns:
+\itemize{
+\item NW: Normalized width calculated as (Upper - Lower) / Estimate
+\item Rating: Categorical rating of the normalized width
+}
+Row names are preserved from input for "ZG" or modified to "EC_x" format
+for "drc".
 }
 \description{
-Title
+Calculates the normalized width of confidence intervals for effect concentration
+(EC) estimates by dividing the interval width by the estimate value. Also
+provides a rating for the calculated normalized width.
+}
+\examples{
+# Example with ZG format
+ec_data <- data.frame(
+  Lower = c(1.2, 2.3),
+  Upper = c(1.8, 3.1),
+  Estimate = c(1.5, 2.7),
+  row.names = c("EC10", "EC50")
+)
+calcNW(ec_data, ED = "ZG")
+
+# Example with drc format
+drc_data <- data.frame(
+  Lower = c(1.2, 2.3),
+  Upper = c(1.8, 3.1),
+  Estimate = c(1.5, 2.7),
+  row.names = c("Dose:1:10", "Dose:1:50")
+)
+calcNW(drc_data, ED = "drc")
+
+}
+\seealso{
+\code{\link{ECx_rating}} for the rating system used
 }
diff --git a/man/drcCompare.Rd b/man/drcCompare.Rd
index 0320757..2c9327b 100644
--- a/man/drcCompare.Rd
+++ b/man/drcCompare.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/drc_Helper.R
 \name{drcCompare}
 \alias{drcCompare}
-\title{Compare models with multiple criteria}
+\title{Compare Dose-Response Models Using Multiple Criteria}
 \usage{
 drcCompare(
   modRes = NULL,
@@ -13,8 +13,40 @@ drcCompare(
 )
 }
 \arguments{
-\item{...}{}
+\item{modRes}{Optional. Result object containing model comparison information.
+If provided, modList parameter is ignored.}
+
+\item{modList}{Optional. List of fitted dose-response models to compare.
+Required if modRes is NULL.}
+
+\item{trend}{Character string specifying the expected trend direction.
+Must be "Decrease" or "Increase".}
+
+\item{CI}{Character string specifying the confidence interval method.
+One of "delta", "inv", or "bmd-inv". Defaults to "delta".}
+
+\item{...}{Additional arguments passed to internal functions}
+}
+\value{
+A data frame containing model comparison results with columns:
+\itemize{
+\item Original comparison metrics
+\item Certainty_Protection: Measure of protection certainty
+\item Steepness: Model steepness evaluation
+\item No Effect p-val: P-value for test against no-effect model
+}
 }
 \description{
-Compare models with multiple criteria
+Performs comprehensive comparison of dose-response models using multiple
+criteria including model fit, certainty of protection, and steepness.
+}
+\examples{
+\dontrun{
+# Compare models using modList
+models <- list(mod1 = fit1, mod2 = fit2)
+drcCompare(modList = models, trend = "Decrease")
+
+# Compare models using existing comparison results
+drcCompare(modRes = previous_results, trend = "Decrease")
+}
 }
diff --git a/man/figures/README-unnamed-chunk-7-1.png b/man/figures/README-unnamed-chunk-7-1.png
new file mode 100644
index 0000000..24c2184
Binary files /dev/null and b/man/figures/README-unnamed-chunk-7-1.png differ
diff --git a/man/get_CA_Z.Rd b/man/get_CA_Z.Rd
index e197199..baedad7 100644
--- a/man/get_CA_Z.Rd
+++ b/man/get_CA_Z.Rd
@@ -25,7 +25,7 @@ uses adjusted values to account for clustering in the data.
 
 The function assigns scores (1, 2, 3, ...) to the treatment groups and
 calculates the Z-statistic using the formula:
-Z =  / sqrt
+\verb{\deqn\{Z = [sum(adj_x*d) - N*p_bar*d_bar] / sqrt[p_bar*(1-p_bar)*(sum(adj_n*d^2) - N*d_bar^2)]}
 where:
 \itemize{
 \item d are the scores (1, 2, 3, ...)
@@ -45,3 +45,6 @@ adj_vals <- get_RS_adj_val(
 # Calculate Z-statistic
 Z <- get_CA_Z(adj_vals$x_tilde, adj_vals$n_tilde)
 }
+\author{
+Originally by Allen Olmstead
+}
diff --git a/man/mselect.ED.Rd b/man/mselect.ED.Rd
index b3f572a..991a9c7 100644
--- a/man/mselect.ED.Rd
+++ b/man/mselect.ED.Rd
@@ -9,7 +9,11 @@ mselect.ED(modList, respLev = c(10, 20, 50), trend = "Decrease", ...)
 \arguments{
 \item{modList}{list of models}
 
-\item{...}{}
+\item{respLev}{response levels or ECx for example, 10, 20, 50}
+
+\item{trend}{"Decrease" or "Incrrease".}
+
+\item{...}{optional to pass to ED function.}
 }
 \value{
 ED result table
diff --git a/man/plot.StepDownRSCABS.Rd b/man/plot.StepDownRSCABS.Rd
new file mode 100644
index 0000000..e1cfe7e
--- /dev/null
+++ b/man/plot.StepDownRSCABS.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RSCABS_AO.R
+\name{plot.StepDownRSCABS}
+\alias{plot.StepDownRSCABS}
+\title{Plot method for StepDownRSCABS objects}
+\usage{
+\method{plot}{StepDownRSCABS}(x, ...)
+}
+\arguments{
+\item{x}{A StepDownRSCABS object}
+
+\item{...}{Additional arguments (not used)}
+}
+\value{
+A ggplot object
+}
+\description{
+Plot method for StepDownRSCABS objects
+}
diff --git a/man/plot.modList.Rd b/man/plot.modList.Rd
index 8314c0b..d077150 100644
--- a/man/plot.modList.Rd
+++ b/man/plot.modList.Rd
@@ -44,6 +44,9 @@ confidence bands}
 user input when the use does not want to use the
 data provided in the modList}
 }
+\value{
+a ggplot object
+}
 \description{
 Plot a list of models together.
 }
diff --git a/man/print.RSCABS.Rd b/man/print.RSCABS.Rd
new file mode 100644
index 0000000..ad00f10
--- /dev/null
+++ b/man/print.RSCABS.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RSCABS_AO.R
+\name{print.RSCABS}
+\alias{print.RSCABS}
+\title{Printing method for run_all_threshold_tests results}
+\usage{
+\method{print}{RSCABS}(x)
+}
+\arguments{
+\item{x}{}
+}
+\value{
+printed results from run_all_threshold_tests
+}
+\description{
+Printing method for run_all_threshold_tests results
+}
diff --git a/man/print.StepDownRSCABS.Rd b/man/print.StepDownRSCABS.Rd
new file mode 100644
index 0000000..2d800ed
--- /dev/null
+++ b/man/print.StepDownRSCABS.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RSCABS_AO.R
+\name{print.StepDownRSCABS}
+\alias{print.StepDownRSCABS}
+\title{Print method for StepDownRSCABS objects}
+\usage{
+\method{print}{StepDownRSCABS}(x, printLowest = FALSE, ...)
+}
+\arguments{
+\item{x}{A StepDownRSCABS object}
+
+\item{...}{Additional arguments (not used)}
+}
+\value{
+Invisibly returns the object
+}
+\description{
+Print method for StepDownRSCABS objects
+}
diff --git a/man/run_RSCA.Rd b/man/run_RSCA.Rd
index 88d1f0e..a0ff1f8 100644
--- a/man/run_RSCA.Rd
+++ b/man/run_RSCA.Rd
@@ -4,7 +4,7 @@
 \alias{run_RSCA}
 \title{Run Rao-Scott Adjusted Cochran-Armitage Trend Test}
 \usage{
-run_RSCA(group, replicate, affected, total)
+run_RSCA(group, replicate, affected, total, correction = 0)
 }
 \arguments{
 \item{group}{Vector of treatment group identifiers}
@@ -14,6 +14,8 @@ run_RSCA(group, replicate, affected, total)
 \item{affected}{Vector of counts of affected subjects (fish with injuries) in each replicate}
 
 \item{total}{Vector of total subjects (fish) in each replicate}
+
+\item{correction}{continuity correction when there is 1, default is 0, can be changed to 0.5.}
 }
 \value{
 A list containing:
@@ -52,3 +54,6 @@ print(result$Z)
 p_value <- 2 * (1 - pnorm(abs(result$Z)))
 print(paste("p-value:", p_value))
 }
+\author{
+Originally by Allen Olmstead
+}
diff --git a/man/run_all_threshold_tests.Rd b/man/run_all_threshold_tests.Rd
new file mode 100644
index 0000000..8da1e4b
--- /dev/null
+++ b/man/run_all_threshold_tests.Rd
@@ -0,0 +1,82 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RSCABS_AO.R
+\name{run_all_threshold_tests}
+\alias{run_all_threshold_tests}
+\title{Run Rao-Scott Adjusted Cochran-Armitage Trend Tests for All Thresholds}
+\usage{
+run_all_threshold_tests(
+  data,
+  max_score = NULL,
+  min_score = 1,
+  score_cols = NULL,
+  treatment_col = "tmt",
+  replicate_col = "tank",
+  total_col = "total",
+  direction = "greater",
+  alternative = "two.sided",
+  include_fisher = TRUE
+)
+}
+\arguments{
+\item{data}{A data frame containing fish injury data}
+
+\item{max_score}{Maximum score value to consider (default: NULL, auto-detected)}
+
+\item{min_score}{Minimum score value to consider (default: 1)}
+
+\item{score_cols}{Character vector of column names containing injury scores (default: NULL, auto-detected)}
+
+\item{treatment_col}{Name of the column containing treatment groups (default: "tmt")}
+
+\item{replicate_col}{Name of the column containing tank/replicate IDs (default: "tank")}
+
+\item{total_col}{Name of the column containing total counts (default: "total")}
+
+\item{direction}{Character string indicating threshold direction: "greater" for ≥ threshold,
+"less" for ≤ threshold (default: "greater")}
+
+\item{alternative}{Character string specifying the alternative hypothesis:
+"two.sided" (default), "greater" (proportion increases with treatment level),
+or "less" (proportion decreases with treatment level)}
+
+\item{include_fisher}{Logical indicating whether to use Fisher's exact test when RSCA fails (default: TRUE)}
+}
+\value{
+A list containing:
+\item{results}{Data frame with test results for each threshold}
+\item{proportions}{Data frame with proportions of affected fish by treatment and threshold}
+\item{detailed_results}{List of detailed results for each threshold}
+}
+\description{
+This function performs the Rao-Scott adjusted Cochran-Armitage trend test for
+multiple injury thresholds, providing a comprehensive analysis of dose-response
+relationships at different severity levels.
+}
+\examples{
+# Example data
+fish_data <- data.frame(
+  tmt = rep(c("Control", "Low", "Medium", "High"), each = 3),
+  tank = paste0(rep(c("Control", "Low", "Medium", "High"), each = 3),
+                rep(1:3, times = 4)),
+  S0 = c(8,7,9, 6,5,7, 4,5,3, 2,3,1),
+  S1 = c(2,2,1, 3,4,2, 4,3,5, 4,3,5),
+  S2 = c(0,1,0, 1,1,1, 2,2,1, 3,3,3),
+  S3 = c(0,0,0, 0,0,0, 0,0,1, 1,1,1),
+  S4 = c(0,0,0, 0,0,0, 0,0,0, 0,0,0),
+  total = rep(10, 12)
+)
+
+# Run two-sided tests for all thresholds
+all_results <- run_all_threshold_tests(fish_data)
+
+# Run one-sided tests (proportion increases with treatment level)
+all_results_greater <- run_all_threshold_tests(
+  fish_data,
+  alternative = "greater"
+)
+
+# View results table
+print(all_results$results)
+print(all_results_greater$results)
+
+}
diff --git a/man/run_threshold_RSCA.Rd b/man/run_threshold_RSCA.Rd
new file mode 100644
index 0000000..a1413b2
--- /dev/null
+++ b/man/run_threshold_RSCA.Rd
@@ -0,0 +1,79 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RSCABS_AO.R
+\name{run_threshold_RSCA}
+\alias{run_threshold_RSCA}
+\title{Run Rao-Scott Adjusted Cochran-Armitage Trend Test for a Specific Injury Threshold}
+\usage{
+run_threshold_RSCA(
+  data,
+  threshold,
+  score_cols = NULL,
+  treatment_col = "tmt",
+  replicate_col = "tank",
+  total_col = "total",
+  direction = "greater",
+  alternative = "two.sided"
+)
+}
+\arguments{
+\item{data}{A data frame containing fish injury data}
+
+\item{threshold}{Numeric value indicating the injury threshold (e.g., 1 for S1+)}
+
+\item{score_cols}{Character vector of column names containing injury scores (default: NULL, auto-detected)}
+
+\item{treatment_col}{Name of the column containing treatment groups (default: "tmt")}
+
+\item{replicate_col}{Name of the column containing tank/replicate IDs (default: "tank")}
+
+\item{total_col}{Name of the column containing total counts (default: "total")}
+
+\item{direction}{Character string indicating threshold direction: "greater" for ≥ threshold,
+"less" for ≤ threshold (default: "greater")}
+
+\item{alternative}{Character string specifying the alternative hypothesis:
+"two.sided" (default), "greater" (proportion increases with treatment level),
+or "less" (proportion decreases with treatment level)}
+}
+\value{
+A list containing:
+\item{threshold}{The numeric threshold value}
+\item{threshold_label}{Character label for the threshold (e.g., "S1+")}
+\item{Z}{Z-statistic from the RSCA test}
+\item{p_value}{P-value based on the specified alternative hypothesis}
+\item{alternative}{The alternative hypothesis used}
+\item{interm_values}{Intermediate values from the Rao-Scott adjustment}
+\item{has_zero_counts}{Logical indicating if any treatment group had zero affected individuals}
+\item{affected_counts}{Data frame showing affected counts by treatment group}
+}
+\description{
+This function performs the Rao-Scott adjusted Cochran-Armitage trend test for
+a specific injury threshold, counting fish with injuries at or above that threshold
+as "affected".
+}
+\examples{
+# Example data
+fish_data <- data.frame(
+  tmt = rep(c("Control", "Low", "Medium", "High"), each = 3),
+  tank = paste0(rep(c("Control", "Low", "Medium", "High"), each = 3),
+                rep(1:3, times = 4)),
+  S0 = c(8,7,9, 6,5,7, 4,5,3, 2,3,1),
+  S1 = c(2,2,1, 3,4,2, 4,3,5, 4,3,5),
+  S2 = c(0,1,0, 1,1,1, 2,2,1, 3,3,3),
+  S3 = c(0,0,0, 0,0,0, 0,0,1, 1,1,1),
+  S4 = c(0,0,0, 0,0,0, 0,0,0, 0,0,0),
+  total = rep(10, 12)
+)
+
+# Two-sided test for trend in S2+ injuries
+result_two_sided <- run_threshold_RSCA(fish_data, threshold = 2)
+
+# One-sided test (proportion increases with treatment level)
+result_greater <- run_threshold_RSCA(fish_data, threshold = 2,
+                                    alternative = "greater")
+
+# One-sided test (proportion decreases with treatment level)
+result_less <- run_threshold_RSCA(fish_data, threshold = 2,
+                                 alternative = "less")
+
+}
diff --git a/man/step_down_RSCABS.Rd b/man/step_down_RSCABS.Rd
new file mode 100644
index 0000000..257b919
--- /dev/null
+++ b/man/step_down_RSCABS.Rd
@@ -0,0 +1,83 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RSCABS_AO.R
+\name{step_down_RSCABS}
+\alias{step_down_RSCABS}
+\title{Perform Step-Down Rao-Scott Adjusted Cochran-Armitage Trend Test Procedure}
+\usage{
+step_down_RSCABS(
+  data,
+  treatment_col = "tmt",
+  treatment_order = NULL,
+  max_score = NULL,
+  min_score = 1,
+  score_cols = NULL,
+  replicate_col = "tank",
+  total_col = "total",
+  direction = "greater",
+  alternative = "greater",
+  include_fisher = TRUE
+)
+}
+\arguments{
+\item{data}{A data frame containing fish injury data}
+
+\item{treatment_col}{Name of the column containing treatment groups (default: "tmt")}
+
+\item{treatment_order}{Optional vector specifying the order of treatment groups from
+lowest to highest dose. If NULL, alphabetical order is used (default: NULL)}
+
+\item{max_score}{Maximum score value to consider (default: NULL, auto-detected)}
+
+\item{min_score}{Minimum score value to consider (default: 1)}
+
+\item{score_cols}{Character vector of column names containing injury scores (default: NULL, auto-detected)}
+
+\item{replicate_col}{Name of the column containing tank/replicate IDs (default: "tank")}
+
+\item{total_col}{Name of the column containing total counts (default: "total")}
+
+\item{direction}{Character string indicating threshold direction: "greater" for ≥ threshold,
+"less" for ≤ threshold (default: "greater")}
+
+\item{alternative}{Character string specifying the alternative hypothesis:
+"two.sided", "greater" (proportion increases with treatment level),
+or "less" (proportion decreases with treatment level) (default: "greater")}
+
+\item{include_fisher}{Logical indicating whether to use Fisher's exact test when RSCA fails (default: TRUE)}
+}
+\value{
+A list of class "StepDownRSCABS" containing:
+\item{combined_results}{Data frame with test results for all steps and thresholds}
+\item{step_results}{List of RSCABS objects for each step}
+\item{summary}{Summary of significant findings by step}
+\item{lowest_significant}{Information about the lowest significant treatment level}
+\item{parameters}{List of parameters used for the analysis}
+}
+\description{
+This function performs a step-down procedure using Rao-Scott adjusted Cochran-Armitage
+trend tests (RSCABS). The procedure systematically excludes the highest treatment groups
+to identify at which treatment level the dose-response relationship becomes significant.
+}
+\examples{
+# Example data
+fish_data <- data.frame(
+  tmt = c(rep("Control", 8), rep("Low", 4), rep("Medium", 4), rep("High", 4)),
+  tank = c(paste0("C", 1:8), paste0("L", 1:4), paste0("M", 1:4), paste0("H", 1:4)),
+  S0 = c(3, 2, 3, 2, 2, 1, 2, 3, 3, 3, 4, 2, 2, 3, 2, 2, 2, 3, 1, 2),
+  S1 = c(1, 2, 0, 1, 2, 2, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 2, 0),
+  S2 = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1),
+  S3 = c(0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1),
+  total = c(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
+)
+
+# Run step-down procedure with default parameters
+result <- step_down_RSCABS(fish_data, treatment_col = "tmt",
+                          treatment_order = c("Control", "Low", "Medium", "High"))
+
+# Print results
+print(result)
+
+# Plot results
+plot(result)
+
+}
diff --git a/man/summary.StepDownRSCABS.Rd b/man/summary.StepDownRSCABS.Rd
new file mode 100644
index 0000000..7619d37
--- /dev/null
+++ b/man/summary.StepDownRSCABS.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RSCABS_AO.R
+\name{summary.StepDownRSCABS}
+\alias{summary.StepDownRSCABS}
+\title{Summary method for StepDownRSCABS objects}
+\usage{
+\method{summary}{StepDownRSCABS}(object, ...)
+}
+\arguments{
+\item{object}{A StepDownRSCABS object}
+
+\item{...}{Additional arguments (not used)}
+}
+\value{
+A list with summary information
+}
+\description{
+Summary method for StepDownRSCABS objects
+}
diff --git a/tests/testthat/test_RSCABS.R b/tests/testthat/test_RSCABS.R
index d72e005..761b908 100644
--- a/tests/testthat/test_RSCABS.R
+++ b/tests/testthat/test_RSCABS.R
@@ -1,4 +1,485 @@
 library(testthat)
+## testthat::test_file("tests/testthat/test_RSCABS.R")
+
+
+# Unit tests for step_down_RSCABS function
+
+describe("step_down_RSCABS", {
+  # Create a sample dataset
+  test_data <- data.frame(
+    tmt = c(rep("Control", 8), rep("Low", 4), rep("Medium", 4), rep("High", 4)),
+    tank = c(paste0("C", 1:8), paste0("L", 1:4), paste0("M", 1:4), paste0("H", 1:4)),
+    S0 = c(3, 2, 3, 2, 2, 1, 2, 3, 3, 3, 4, 2, 2, 3, 2, 2, 2, 3, 1, 2),
+    S1 = c(1, 2, 0, 1, 2, 2, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 2, 0),
+    S2 = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1),
+    S3 = c(0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1),
+    total = c(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
+  )
+  it("returns an object of class StepDownRSCABS", {
+    result <- suppressWarnings({step_down_RSCABS(
+      test_data,
+      treatment_order = c("Control", "Low", "Medium", "High")
+    )
+    })
+
+    expect_true(inherits(result, "StepDownRSCABS"))
+    expect_equal(names(result), c("combined_results", "step_results", "summary",
+                                  "lowest_significant", "parameters"))
+  })
+
+  it("performs the correct number of steps", {
+    result <- suppressWarnings({step_down_RSCABS(
+      test_data,
+      treatment_order = c("Control", "Low", "Medium", "High")
+    )})
+
+    # Should have 3 steps: All, No High, Control+Low only
+    expect_equal(length(result$step_results), 3)
+    expect_equal(length(result$summary), 3)
+    expect_equal(unique(result$combined_results$Step), 1:3)
+  })
+
+  it("correctly identifies the lowest significant treatment level", {
+    result <- suppressWarnings({step_down_RSCABS(
+      test_data,
+      treatment_order = c("Control", "Low", "Medium", "High")
+    )})
+
+
+    expect_equal(result$lowest_significant$step, NA)
+    expect_equal(result$lowest_significant$treatment, NA)
+  })
+
+    it("handles custom treatment order", {
+      # Test with reversed treatment order
+      result <- suppressWarnings({step_down_RSCABS(
+        test_data,
+        treatment_order = c("High", "Medium", "Low", "Control")
+      )})
+
+      # Check that the steps follow the specified order
+      expect_equal(result$summary[[1]]$highest_treatment, "Control")
+      expect_equal(result$summary[[2]]$highest_treatment, "Low")
+      expect_equal(result$summary[[3]]$highest_treatment, "Medium")
+    })
+
+
+  it("validates input parameters", {
+    # Missing treatment column
+    expect_error(
+      step_down_RSCABS(test_data, treatment_col = "treatment"),
+      "Treatment column 'treatment' not found in data"
+    )
+
+    # Invalid treatment order
+    expect_error(
+      step_down_RSCABS(test_data, treatment_order = c("Control", "Low", "High")),
+      "Some treatments in data are not in treatment_order"
+    )
+
+    expect_error(
+      step_down_RSCABS(test_data, treatment_order = c("Control", "Low", "Medium", "High", "Extra")),
+      "Some treatments in treatment_order are not in data"
+    )
+
+    # Only one treatment level
+    single_treatment_data <- test_data[test_data$tmt == "Control", ]
+    expect_error(
+      step_down_RSCABS(single_treatment_data),
+      "At least 2 treatment levels are required"
+    )
+  })
+})
+
+describe("print.StepDownRSCABS", {
+  # Create a simple StepDownRSCABS object for testing
+  test_object <- list(
+    combined_results = data.frame(
+      Threshold = rep(c("S1+", "S2+", "S3+"), 3),
+      Z_statistic = rep(2:4, each = 3),
+      P_value = c(0.03, 0.04, 0.06, 0.05, 0.07, 0.08, 0.1, 0.2, 0.3),
+      Step = rep(1:3, each = 3),
+      Included_Treatments = rep(c("Control, Low, Medium, High",
+                                  "Control, Low, Medium",
+                                  "Control, Low"), each = 3),
+      Highest_Treatment = rep(c("High", "Medium", "Low"), each = 3)
+    ),
+    step_results = list(),
+    summary = list(
+      list(
+        step = 1,
+        included_treatments = "Control, Low, Medium, High",
+        highest_treatment = "High",
+        significant_thresholds = c("S1+", "S2+"),
+        has_significant_findings = TRUE
+      ),
+      list(
+        step = 2,
+        included_treatments = "Control, Low, Medium",
+        highest_treatment = "Medium",
+        significant_thresholds = "S1+",
+        has_significant_findings = TRUE
+      ),
+      list(
+        step = 3,
+        included_treatments = "Control, Low",
+        highest_treatment = "Low",
+        significant_thresholds = character(0),
+        has_significant_findings = FALSE
+      )
+    ),
+    lowest_significant = list(
+      step = 2,
+      treatment = "Medium",
+      threshold = "S1+",
+      p_value = 0.05
+    ),
+    parameters = list(
+      treatment_col = "tmt",
+      treatment_order = c("Control", "Low", "Medium", "High"),
+      direction = "greater",
+      alternative = "greater"
+    )
+  )
+  class(test_object) <- "StepDownRSCABS"
+
+  it("prints without error", {
+    expect_output(print(test_object), "Step-Down RSCABS Analysis")
+    expect_output(print(test_object,printLowest = T), "Lowest treatment level with significant findings")
+    expect_output(print(test_object,printLowest = T), "Treatment: Medium")
+  })
+})
+
+describe("summary.StepDownRSCABS", {
+  # Create a simple StepDownRSCABS object for testing
+  test_object <- list(
+    combined_results = data.frame(
+      Threshold = rep(c("S1+", "S2+", "S3+"), 3),
+      Z_statistic = rep(2:4, each = 3),
+      P_value = c(0.03, 0.04, 0.06, 0.05, 0.07, 0.08, 0.1, 0.2, 0.3),
+      Step = rep(1:3, each = 3),
+      Included_Treatments = rep(c("Control, Low, Medium, High",
+                                  "Control, Low, Medium",
+                                  "Control, Low"), each = 3),
+      Highest_Treatment = rep(c("High", "Medium", "Low"), each = 3)
+    ),
+    step_results = list(),
+    summary = list(
+      list(
+        step = 1,
+        included_treatments = "Control, Low, Medium, High",
+        highest_treatment = "High",
+        significant_thresholds = c("S1+", "S2+"),
+        has_significant_findings = TRUE
+      ),
+      list(
+        step = 2,
+        included_treatments = "Control, Low, Medium",
+        highest_treatment = "Medium",
+        significant_thresholds = "S1+",
+        has_significant_findings = TRUE
+      ),
+      list(
+        step = 3,
+        included_treatments = "Control, Low",
+        highest_treatment = "Low",
+        significant_thresholds = character(0),
+        has_significant_findings = FALSE
+      )
+    ),
+    lowest_significant = list(
+      step = 2,
+      treatment = "Medium",
+      threshold = "S1+",
+      p_value = 0.05
+    ),
+    parameters = list(
+      treatment_col = "tmt",
+      treatment_order = c("Control", "Low", "Medium", "High"),
+      direction = "greater",
+      alternative = "greater"
+    )
+  )
+  class(test_object) <- "StepDownRSCABS"
+
+  it("returns a summary with the correct structure", {
+    result <- summary(test_object)
+
+    expect_true(is.list(result))
+    expect_equal(names(result), c("summary_table", "lowest_significant", "parameters"))
+    expect_equal(nrow(result$summary_table), 3)
+    expect_equal(result$lowest_significant$treatment, "Medium")
+  })
+})
+
+describe("plot.StepDownRSCABS", {
+  # Skip if ggplot2 is not available
+  skip_if_not_installed("ggplot2")
+
+  # Create a simple StepDownRSCABS object for testing
+  test_object <- list(
+    combined_results = data.frame(
+      Threshold = rep(c("S1+", "S2+", "S3+"), 3),
+      Z_statistic = rep(2:4, each = 3),
+      P_value = c(0.03, 0.04, 0.06, 0.05, 0.07, 0.08, 0.1, 0.2, 0.3),
+      Step = rep(1:3, each = 3),
+      Included_Treatments = rep(c("Control, Low, Medium, High",
+                                  "Control, Low, Medium",
+                                  "Control, Low"), each = 3),
+      Highest_Treatment = rep(c("High", "Medium", "Low"), each = 3)
+    ),
+    step_results = list(),
+    summary = list(
+      list(
+        step = 1,
+        included_treatments = "Control, Low, Medium, High",
+        highest_treatment = "High",
+        significant_thresholds = c("S1+", "S2+"),
+        has_significant_findings = TRUE
+      ),
+      list(
+        step = 2,
+        included_treatments = "Control, Low, Medium",
+        highest_treatment = "Medium",
+        significant_thresholds = "S1+",
+        has_significant_findings = TRUE
+      ),
+      list(
+        step = 3,
+        included_treatments = "Control, Low",
+        highest_treatment = "Low",
+        significant_thresholds = character(0),
+        has_significant_findings = FALSE
+      )
+    ),
+    lowest_significant = list(
+      step = 2,
+      treatment = "Medium",
+      threshold = "S1+",
+      p_value = 0.05
+    ),
+    parameters = list(
+      treatment_col = "tmt",
+      treatment_order = c("Control", "Low", "Medium", "High"),
+      direction = "greater",
+      alternative = "greater"
+    )
+  )
+  class(test_object) <- "StepDownRSCABS"
+
+  it("returns a ggplot object", {
+    result <- plot(test_object)
+
+    expect_true(inherits(result, "gg"))
+    expect_true(inherits(result, "ggplot"))
+  })
+})
+
+
+
+# Unit tests for RSCA threshold functions
+
+# Test for run_threshold_RSCA
+describe("run_threshold_RSCA", {
+  # Create a sample dataset
+  test_data <- data.frame(
+    tmt = rep(c("Control", "Low", "Medium", "High"), each = 3),
+    tank = paste0(rep(c("Control", "Low", "Medium", "High"), each = 3),
+                  rep(1:3, times = 4)),
+    S0 = c(8,7,9, 6,5,7, 4,5,3, 2,3,1),
+    S1 = c(2,2,1, 3,4,2, 4,3,5, 4,3,5),
+    S2 = c(0,1,0, 1,1,1, 2,2,1, 3,3,3),
+    S3 = c(0,0,0, 0,0,0, 0,0,1, 1,1,1),
+    S4 = c(0,0,0, 0,0,0, 0,0,0, 0,0,0),
+    total = rep(10, 12)
+  )
+
+  # Dataset with zero counts
+  zero_data <- data.frame(
+    tmt = rep(c("Control", "Low", "High"), each = 3),
+    tank = paste0(rep(c("Control", "Low", "High"), each = 3),
+                  rep(1:3, times = 3)),
+    S0 = c(10,10,10, 8,9,8, 5,6,5),
+    S1 = c(0,0,0, 2,1,2, 3,2,3),
+    S2 = c(0,0,0, 0,0,0, 2,2,2),
+    total = rep(10, 9)
+  )
+
+  it("correctly calculates RSCA test for a valid threshold", {
+    # Test S1+ threshold (all groups have non-zero counts)
+    result <- run_threshold_RSCA(test_data, threshold = 1)
+
+    # Check structure
+    expect_equal(result$threshold, 1)
+    expect_equal(result$threshold_label, "S1+")
+    expect_false(is.na(result$Z))
+    expect_false(is.na(result$p_value))
+
+    expect_false(result$has_zero_counts) ## the run_threshold_RSCA function does not return has_zero_counts directly.
+
+    # Check affected counts
+    expect_equal(nrow(result$interm_values), 4)  # 4 treatment groups
+    expect_true(all(result$interm_values$x > 0))  # All groups have affected fish
+  })
+
+  it("handles thresholds with zero counts", {
+    # Test S3+ threshold (some groups have zero counts)
+    expect_warning(
+      result <- run_threshold_RSCA(test_data, threshold = 3),
+      "Some treatment groups have zero affected individuals"
+    )
+
+    # Check that it flags zero counts
+    expect_true(result$has_zero_counts)
+
+    # Check affected counts
+    expect_equal(sum(result$affected_counts$affected == 0), 2)  # Control and Low have zero S3+ fish
+  })
+
+  it("works with 'less than or equal to' direction", {
+    # Test S≤1 threshold
+    result <- run_threshold_RSCA(test_data, threshold = 1, direction = "less")
+
+    # Check structure
+    expect_equal(result$threshold_label, "S≤1")
+    expect_false(is.na(result$Z))
+    expect_false(is.na(result$p_value))
+
+    # Check that counts make sense
+    total_s0_s1 <- sum(test_data$S0) + sum(test_data$S1)
+    expect_equal(sum(result$affected_counts$affected), total_s0_s1)
+  })
+
+  it("auto-detects score columns", {
+    # Create data with non-standard column names but keep S columns
+    modified_data <- test_data
+    modified_data$treatment <- modified_data$tmt
+    modified_data$replicate <- modified_data$tank
+    modified_data$count <- modified_data$total
+    modified_data$tmt <- NULL
+    modified_data$tank <- NULL
+    modified_data$total <- NULL
+
+    # Should still find S0-S4 columns
+    result <- run_threshold_RSCA(
+      modified_data,
+      threshold = 1,
+      treatment_col = "treatment",
+      replicate_col = "replicate",
+      total_col = "count"
+    )
+
+    expect_false(is.na(result$Z))
+  })
+
+  it("raises error with invalid inputs", {
+    # Missing treatment column
+    expect_error(
+      run_threshold_RSCA(test_data, threshold = 1, treatment_col = "treatment"),
+      "Treatment column 'treatment' not found in data"
+    )
+
+    # Invalid direction
+    expect_error(
+      run_threshold_RSCA(test_data, threshold = 1, direction = "invalid"),
+      "Direction must be either 'greater' or 'less'"
+    )
+
+    # Invalid score columns
+    expect_error(
+      run_threshold_RSCA(test_data, threshold = 1, score_cols = c("S0", "S1", "S99")),
+      "Missing score columns"
+    )
+  })
+})
+
+# Test for run_all_threshold_tests
+describe("run_all_threshold_tests", {
+  # Create a sample dataset
+  test_data <- data.frame(
+    tmt = rep(c("Control", "Low", "Medium", "High"), each = 3),
+    tank = paste0(rep(c("Control", "Low", "Medium", "High"), each = 3),
+                  rep(1:3, times = 4)),
+    S0 = c(8,7,9, 6,5,7, 4,5,3, 2,3,1),
+    S1 = c(2,2,1, 3,4,2, 4,3,5, 4,3,5),
+    S2 = c(0,1,0, 1,1,1, 2,2,1, 3,3,3),
+    S3 = c(0,0,0, 0,0,0, 0,0,1, 1,1,1),
+    S4 = c(0,0,0, 0,0,0, 0,0,0, 0,0,0),
+    total = rep(10, 12)
+  )
+
+  it("runs tests for all thresholds", {
+    # Run tests for S1+ through S4+
+    result <- suppressWarnings(run_all_threshold_tests(test_data)) ## warning is expected
+
+    # Check structure
+    expect_equal(nrow(result$results), 4)  # 4 thresholds (S1+ through S4+)
+    expect_equal(result$results$Threshold, c("S1+", "S2+", "S3+", "S4+"))
+
+    # Check that proportions table has correct structure
+    expect_equal(nrow(result$proportions), 4)  # 4 treatment groups
+    expect_equal(ncol(result$proportions), 5)  # Treatment + 4 thresholds
+
+    # Check that detailed results are included
+    expect_equal(length(result$detailed_results), 4)
+  })
+
+  it("handles custom min and max scores", {
+    # Run tests only for S2+ and S3+
+    result <- suppressWarnings(run_all_threshold_tests(test_data, min_score = 2, max_score = 3))
+
+    # Check structure
+    expect_equal(nrow(result$results), 2)  # 2 thresholds
+    expect_equal(result$results$Threshold, c("S2+", "S3+"))
+  })
+
+  it("works with 'less than or equal to' direction", {
+    # Run tests for S≤0 through S≤4
+    result <- run_all_threshold_tests(
+      test_data,
+      min_score = 0,
+      max_score = 4,
+      direction = "less"
+    )
+
+    # Check structure
+    expect_equal(nrow(result$results), 5)  # 5 thresholds
+    expect_equal(result$results$Threshold, c("S≤0", "S≤1", "S≤2", "S≤3", "S≤4"))
+
+    # Check proportions make sense
+    # S≤4 should be 1.0 for all treatments (all fish have S≤4)
+    expect_equal(unique(result$proportions$`S≤4`), 1)
+  })
+
+  it("includes Fisher's exact test for thresholds with zero counts", {
+    # Run with Fisher's test for thresholds with zero counts
+    result <- suppressWarnings(run_all_threshold_tests(test_data, include_fisher = TRUE))
+
+    # Check that methods are assigned correctly
+    expect_true("Method" %in% names(result$results))
+    expect_true(any( grepl("RSCA",result$results$Method)))
+
+    # For thresholds with zero counts, should use Fisher's test or mark as not calculable
+    zero_count_rows <- which(result$results$Has_zero_counts)
+    if (length(zero_count_rows) > 0) {
+      expect_true(all(grepl("Fisher's Exact Test",result$results$Method[zero_count_rows]) |
+                        grepl("Not calculable",result$results$Method[zero_count_rows])) )
+    }
+  })
+
+  it("auto-detects max score", {
+    # Create data with only S0-S3 columns
+    modified_data <- test_data[, c("tmt", "tank", "S0", "S1", "S2", "S3", "total")]
+
+    # Should detect max score of 3
+    result <- suppressWarnings(run_all_threshold_tests(modified_data))
+
+    # Check structure
+    expect_equal(nrow(result$results), 3)  # 3 thresholds (S1+ through S3+)
+    expect_equal(result$results$Threshold, c("S1+", "S2+", "S3+"))
+  })
+})
+
 describe("RSCABS",{
   # Mock data for testing
   set.seed(123)
@@ -50,3 +531,8 @@ describe("RSCABS",{
   })
 
 })
+
+
+
+
+
diff --git a/vignettes/articles/Count_Data.Rmd b/vignettes/articles/Count_Data.Rmd
index fdb8f8a..9ede6f3 100644
--- a/vignettes/articles/Count_Data.Rmd
+++ b/vignettes/articles/Count_Data.Rmd
@@ -6,6 +6,24 @@ date: December 4, 2024
 author: "Zhenglei Gao"
 ---
 
+*To be written. Contribution is welcome.*
+
+## Testing Approaches
+
+###  Fisher's Exact Test for Count Data
+
+```{r}
+
+```
+
+### Data Transformation
+
+
+## Modelling Approaches
+
+
+
+
 
 ## Poisson, Quasi-Poisson, Negative-Binomial
 
diff --git a/vignettes/articles/Example_RSCABS.Rmd b/vignettes/articles/Example_RSCABS.Rmd
index 13c4d2b..d42244d 100644
--- a/vignettes/articles/Example_RSCABS.Rmd
+++ b/vignettes/articles/Example_RSCABS.Rmd
@@ -8,12 +8,358 @@ editor_options:
 ---
 
 
+## Example Usage of the RSCABS functions in drcHelper
 
 
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  message = FALSE,
+  warning = FALSE
+)
+```
+
+```{r setup}
+library(drcHelper)
+library(tidyverse)
+```
+
+
+The function in *drcHelper* is `run_threshold_RSCA` and `run_all_threshold_tests`. The latter performs the Rao-Scott adjusted Cochran-Armitage trend test by slices (for multiple injury thresholds) and 
+providing a comprehensive analysis of dose-response relationships at different severity levels. It also outputs an invisible components including all the detailed interim results for each tested threshold level. A step-down procedure can be performed using `step_down_RSCABS`. 
+
+
+
+```{r}
+# Example data with increasing trend in injury rates
+# Create the simulated fish data
+sim_data_1<- tibble::tibble(
+  treatment = c(rep("Control", 8), rep("Low", 4), rep("Medium", 4), rep("High", 4)),
+  tank = c(paste0("C", 1:8), paste0("L", 1:4), paste0("M", 1:4), paste0("H", 1:4)),
+  S0 = c(3, 2, 4, 3, 2, 1, 2, 3, 3, 3, 4, 2, 2, 3, 2, 2, 2, 3, 1, 2),
+  S1 = c(1, 2, 0, 1, 2, 2, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 2, 0),
+  S2 = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1),
+  S3 = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1),
+  total = c(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
+)
+
+# Ensure treatment is an ordered factor for proper plotting
+sim_data_1$treatment <- factor(sim_data_1$treatment,
+                                        levels = c("Control", "Low", "Medium", "High"))
+
+
+# Run two-sided test
+two_sided_result <- run_threshold_RSCA(sim_data_1, threshold = 2,treatment_col = "treatment",replicate_col = "tank")
+print(paste("Two-sided p-value:", round(two_sided_result$p_value, 4)))
+
+# Run one-sided test (greater)
+greater_result <- run_threshold_RSCA(sim_data_1, threshold = 2, alternative = "greater",treatment_col = "treatment",replicate_col = "tank")
+print(paste("One-sided (greater) p-value:", round(greater_result$p_value, 4)))
+
+# Run one-sided test (less)
+less_result <- run_threshold_RSCA(sim_data_1, threshold = 1, alternative = "less",treatment_col = "treatment",replicate_col = "tank")
+print(paste("One-sided (less) p-value:", round(less_result$p_value, 4))) 
+```
+
+```{r}
+run_all_threshold_tests(sim_data_1,min_score = 0,max_score = 2,direction="less",alternative ="less",treatment_col = "treatment",replicate_col = "tank")
+run_all_threshold_tests(sim_data_1,min_score = 1,max_score = 3,direction="greater",alternative ="greater",treatment_col = "treatment",replicate_col = "tank")
+run_all_threshold_tests(sim_data_1%>%filter(treatment!="High")%>%droplevels(.),min_score = 1,max_score = 3,direction="greater",alternative ="greater",treatment_col = "treatment",replicate_col = "tank")
+```
 
+```{r}
+# Run step-down procedure
+result <- step_down_RSCABS(
+  sim_data_1,
+  treatment_col = "treatment",replicate_col = "tank",
+  treatment_order = c("Control", "Low", "Medium", "High"),
+  alternative = "greater"
+)
+result
+summary_result <- summary(result)
+print(summary_result$summary_table)
+plot(result)
 
+print(result,printLowest = T)
+```
 
-## RSCABS
+### Comparison with Joe's Implementation
+
+Note the p-values are the same for the critical treatment and score combinations. I do think Joe's implementation has some fault in terms of alternative hypothesis and p-value calculation (It can return 1-p sometimes, not sure why). But I haven't tried to figured that out.  
+
+```{r}
+dat1 <- convert_fish_data(sim_data_1,direction="to_individual",treatment_col = "treatment",replicate_col = "tank" )
+dat1 <- (dat1) %>% mutate(score1=as.numeric(factor(score))-1,score2=score1+5) %>% dplyr::select(-c(score)) %>% as.data.frame
+
+## Note runRSCABS expects more than one endpoints. 
+res <- runRSCABS(dat1,'treatment','tank',test.type='RS')
+res[1:5,] %>% mutate(Effect = gsub("score1","S",Effect) )
+```
+
+
+
+
+Note that you need to specify both direction and alternative according to your test hypotheses. The drcHelper function provides alternative specification by 
+
+
+- Two-sided test (alternative = "two.sided"):
+      -Tests whether there is any trend (increasing or decreasing) in the proportion of affected fish across treatment groups
+      - Null hypothesis (H0): No trend in proportions
+      - Alternative hypothesis (H1): There is a trend (either increasing or decreasing)
+      - P-value calculation: 2 * pnorm(-abs(Z))
+      - Use when you have no prior expectation about the direction of the effect
+
+- One-sided test (greater) (alternative = "greater"):
+      - Tests whether the proportion of affected fish increases with treatment level
+      - Null hypothesis (H0): No increasing trend
+      - Alternative hypothesis (H1): Proportion increases with treatment level
+      - P-value calculation: pnorm(-Z)
+      - Use when you expect higher treatment levels to have more affected fish
+
+- One-sided test (less) (alternative = "less"):
+     - Tests whether the proportion of affected fish decreases with treatment level
+     - Null hypothesis (H0): No decreasing trend
+     - Alternative hypothesis (H1): Proportion decreases with treatment level
+     - P-value calculation: pnorm(Z)
+     - Use when you expect higher treatment levels to have fewer affected fish
+
+In toxicology studies, alternative = "greater" in combination with direction = "greater" is often most appropriate since higher treatment levels are typically expected to cause more adverse effects.
+
+In comparison with the conventional RSCABS test by `runRSCABS`, the `run_all_threshold_tests` provides:
+
+- Flexible threshold direction: Supports both "greater than or equal to" and "less than or equal to" thresholds
+- Automatic score column detection: Identifies score columns based on naming pattern
+- Zero count handling: Provides warnings and alternative approaches for thresholds with zero counts
+- Comprehensive results: Returns detailed information including proportions, test statistics, and intermediate values
+- Fisher's exact test fallback: Optionally uses Fisher's exact test when RSCA is not valid due to zero counts
+
+
+## Visulaize the data
+
+
+```{r echo=FALSE,fig.width=15,fig.height=10}
+library(ggplot2)
+library(dplyr)
+library(tidyr)
+library(patchwork)  # For combining plots
+simulated_fish_data <- sim_data_1
+
+# Convert to long format for easier plotting
+fish_data_long <- simulated_fish_data %>%
+  pivot_longer(
+    cols = starts_with("S"),
+    names_to = "score",
+    values_to = "count"
+  ) %>%
+  # Ensure score is an ordered factor
+  mutate(score = factor(score, levels = c("S0", "S1", "S2", "S3")))
+
+# 1. Stacked Bar Plot by Treatment Group
+p1 <- fish_data_long %>%
+  group_by(treatment, score) %>%
+  summarize(total_count = sum(count), .groups = "drop") %>%
+  ggplot(aes(x = treatment, y = total_count, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  labs(
+    title = "Distribution of Injury Scores by Treatment Group",
+    x = "Treatment Group",
+    y = "Number of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "right"
+  )
+
+# 2. Proportion Plot by Treatment Group
+p2 <- fish_data_long %>%
+  group_by(treatment, score) %>%
+  summarize(total_count = sum(count), .groups = "drop") %>%
+  group_by(treatment) %>%
+  mutate(proportion = total_count / sum(total_count)) %>%
+  ggplot(aes(x = treatment, y = proportion, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  labs(
+    title = "Proportion of Injury Scores by Treatment Group",
+    x = "Treatment Group",
+    y = "Proportion of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "right"
+  )
+
+# 3. Boxplot of Injury Proportions by Treatment
+p3 <- simulated_fish_data %>%
+  mutate(
+    injury_prop = (S1 + S2 + S3) / total,
+    severe_injury_prop = (S2 + S3) / total
+  ) %>%
+  pivot_longer(
+    cols = c(injury_prop, severe_injury_prop),
+    names_to = "injury_type",
+    values_to = "proportion"
+  ) %>%
+  mutate(injury_type = factor(injury_type,
+                              levels = c("injury_prop", "severe_injury_prop"),
+                              labels = c("Any Injury (S1-S3)", "Severe Injury (S2-S3)"))) %>%
+  ggplot(aes(x = treatment, y = proportion, fill = treatment)) +
+  geom_boxplot(alpha = 0.7) +
+  geom_jitter(width = 0.2, height = 0, alpha = 0.7) +
+  facet_wrap(~ injury_type) +
+  labs(
+    title = "Distribution of Injury Proportions by Treatment Group",
+    x = "Treatment Group",
+    y = "Proportion of Fish",
+    fill = "Treatment"
+  ) +
+  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
+  scale_fill_brewer(palette = "Set2") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "none"
+  )
+
+# 4. Heat Map of Injury Counts by Tank
+p4 <- fish_data_long %>%
+  ggplot(aes(x = tank, y = score, fill = count)) +
+  geom_tile(color = "white") +
+  geom_text(aes(label = count), color = "black") +
+  facet_grid(. ~ treatment, scales = "free_x", space = "free_x") +
+  labs(
+    title = "Heat Map of Injury Counts by Tank",
+    x = "Tank",
+    y = "Injury Score",
+    fill = "Count"
+  ) +
+  scale_fill_gradient(low = "white", high = "#fc8d62") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    axis.text.x = element_text(angle = 45, hjust = 1),
+    panel.grid = element_blank(),
+    legend.position = "right"
+  )
+
+# 5. Dot plot showing proportion of each injury type by treatment
+p5 <- fish_data_long %>%
+  group_by(treatment, score) %>%
+  summarize(
+    total_count = sum(count),
+    .groups = "drop"
+  ) %>%
+  group_by(treatment) %>%
+  mutate(proportion = total_count / sum(total_count)) %>%
+  ggplot(aes(x = treatment, y = proportion, color = score)) +
+  geom_point(size = 4) +
+  geom_line(aes(group = score), linewidth = 1) +
+  facet_wrap(~ score, nrow = 1) +
+  labs(
+    title = "Trend in Proportion of Each Injury Score Across Treatments",
+    x = "Treatment Group",
+    y = "Proportion of Fish",
+    color = "Injury Score"
+  ) +
+  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
+  scale_color_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    legend.position = "none",
+    panel.spacing = unit(1, "lines")
+  )
+
+# 6. Stacked bar plot for each tank
+p6 <- fish_data_long %>%
+  ggplot(aes(x = tank, y = count, fill = score)) +
+  geom_bar(stat = "identity", position = "stack") +
+  facet_grid(. ~ treatment, scales = "free_x", space = "free_x") +
+  labs(
+    title = "Distribution of Injury Scores by Individual Tank",
+    x = "Tank",
+    y = "Number of Fish",
+    fill = "Injury Score"
+  ) +
+  scale_fill_brewer(palette = "YlOrRd") +
+  theme_minimal() +
+  theme(
+    plot.title = element_text(size = 12),
+    axis.text.x = element_text(angle = 45, hjust = 1),
+    legend.position = "right"
+  )
+
+# Combine plots using patchwork
+top_row <- p1 + p2
+middle_row <- p3
+bottom_row <- p4 / p5
+last_row <- p6
+
+combined_plot <- (top_row) / (middle_row) / (bottom_row) / (last_row) +
+  plot_layout(heights = c(1, 1, 2, 1)) +
+  plot_annotation(
+    title = "Visualization of Simulated Fish Injury Data",
+    subtitle = "Multiple visualizations showing the distribution of injury scores across treatment groups and tanks",
+    theme = theme(
+      plot.title = element_text(size = 16, face = "bold"),
+      plot.subtitle = element_text(size = 12)
+    )
+  )
+
+# Print the combined plot
+print(combined_plot)
+```
+
+### Smmary table of the data
+
+```{r echo=FALSE}
+# Statistical summary table
+summary_table <- simulated_fish_data %>%
+  group_by(treatment) %>%
+  summarize(
+    n_tanks = n(),
+    total_fish = sum(total),
+    s0_count = sum(S0),
+    s0_percent = 100 * s0_count / total_fish,
+    s1_count = sum(S1),
+    s1_percent = 100 * s1_count / total_fish,
+    s2_count = sum(S2),
+    s2_percent = 100 * s2_count / total_fish,
+    s3_count = sum(S3),
+    s3_percent = 100 * s3_count / total_fish,
+    any_injury = sum(S1 + S2 + S3),
+    any_injury_percent = 100 * any_injury / total_fish,
+    .groups = "drop"
+  ) %>%
+  mutate(across(ends_with("percent"), ~ round(., 1)))
+
+# Print the summary table
+print(summary_table)
+
+# Calculate injury severity index (weighted average of injury scores)
+severity_by_treatment <- simulated_fish_data %>%
+  mutate(
+    severity_index = (0*S0 + 1*S1 + 2*S2 + 3*S3) / total
+  ) %>%
+  group_by(treatment) %>%
+  summarize(
+    mean_severity = mean(severity_index),
+    sd_severity = sd(severity_index),
+    .groups = "drop"
+  )
+
+# Print severity index
+print(severity_by_treatment)
+```
+
+
+## Understanding RSCABS
 
 RSCABS (Rao-Scott Cochran-Armitage by slice) is designed to analyze histopathological results from standard toxicology experiments, for example the MEOGRT.
 
@@ -34,19 +380,7 @@ The Z-statistic calculation breaks down. Possible solutions could be using alter
 
 
 
-```{r, include = FALSE}
-knitr::opts_chunk$set(
-  collapse = TRUE,
-  comment = "#>",
-  message = FALSE,
-  warning = FALSE
-)
-```
 
-```{r setup}
-library(drcHelper)
-library(tidyverse)
-```
 
 
 ### Some Backgrounds
@@ -254,83 +588,6 @@ library(tidyr)
 library(knitr)
 library(ggplot2)
 
-# Function to run RSCA test for a specific injury threshold
-run_threshold_RSCA <- function(data, threshold, score_cols = c("S0", "S1", "S2", "S3", "S4")) {
-  # Determine which scores are above the threshold
-  above_cols <- score_cols[as.numeric(gsub("S", "", score_cols)) >= threshold]
-  
-  # Calculate affected counts (fish with injuries at or above threshold)
-  affected <- rowSums(data[, above_cols, drop = FALSE])
-  
-  # Run the RSCA test
-  result <- run_RSCA(
-    group = data$tmt,
-    replicate = data$tank,
-    affected = affected,
-    total = data$total
-  )
-  
-  # Calculate p-value
-  p_value <- 2 * (1 - pnorm(abs(result$Z)))
-  
-  # Return results
-  list(
-    threshold = threshold,
-    threshold_label = paste0("S", threshold, "+"),
-    Z = result$Z,
-    p_value = p_value,
-    interm_values = result$interm_values
-  )
-}
-
-# Run RSCA tests for all thresholds on the simulated data
-run_all_threshold_tests <- function(data, max_score = 4) {
-  # Create a list to store results
-  results_list <- list()
-  
-  # Score column names
-  score_cols <- paste0("S", 0:max_score)
-  
-  # Run tests for each threshold (1 to max_score)
-  for (threshold in 1:max_score) {
-    results_list[[threshold]] <- run_threshold_RSCA(data, threshold, score_cols)
-  }
-  
-  # Compile results into a data frame
-  results_df <- data.frame(
-    Threshold = sapply(results_list, function(x) x$threshold_label),
-    Z_statistic = sapply(results_list, function(x) x$Z),
-    P_value = sapply(results_list, function(x) x$p_value)
-  )
-  
-  # Calculate proportions for each threshold by treatment
-  prop_by_treatment <- data.frame(Treatment = unique(data$tmt))
-  
-  for (threshold in 1:max_score) {
-    above_cols <- score_cols[as.numeric(gsub("S", "", score_cols)) >= threshold]
-    
-    # Calculate proportions by treatment
-    props <- data %>%
-      group_by(tmt) %>%
-      summarize(
-        affected = sum(rowSums(across(all_of(above_cols)))),
-        total = sum(total),
-        proportion = affected / total
-      ) %>%
-      select(tmt, proportion)
-    
-    # Add to the data frame
-    prop_by_treatment[[paste0("S", threshold, "+")]] <- 
-      props$proportion[match(prop_by_treatment$Treatment, props$tmt)]
-  }
-  
-  # Return both results and proportions
-  list(
-    results = results_df,
-    proportions = prop_by_treatment,
-    detailed_results = results_list
-  )
-}
 
 # Run the analysis on the simulated data
 threshold_results <- run_all_threshold_tests(sim_data)
@@ -514,7 +771,6 @@ dat1 <- exampleHistData.Sub[,c(1:5,5+cids)]%>%tidyr::pivot_longer(-(1:5),values_
 ggplot(dat1,aes(x=Treatment,fill=factor(Response)))+geom_bar(aes(y=counts/total),stat = "identity")+scale_fill_viridis_d()+labs(title="Example Histopath Data",subtitle = "subset: F2 generation, 15 week age and female")+facet_wrap(~Endpoint,drop = T)+ scale_y_continuous(labels = percent)
 ```
 
-## Using R's prop.test instead
 
 
 
diff --git a/vignettes/drcHelper.Rmd b/vignettes/drcHelper.Rmd
index 4396297..6f20800 100644
--- a/vignettes/drcHelper.Rmd
+++ b/vignettes/drcHelper.Rmd
@@ -72,17 +72,45 @@ p
 ## Adding ECx and ECx CI's to the plots
 
 ```{r}
+p1 <- plot.modList(modList[1])
+addECxCI(p1,object=modList[[1]],EDres=NULL,trend="Decrease",endpoint="EC", respLev=c(10,20,50),
+                     textAjust.x=0.01,textAjust.y=0.3,useObsCtr=FALSE,d0=NULL,textsize = 4,lineheight = 0.5,xmin=0.012)+ ylab("Response Variable [unit]") + xlab("Concentration [µg a.s./L]")
 ## addECxCI(p)
 ```
 
+## Report ECx
+
+```{r}
+resED <- t(edResTab[1:3, c(2,4,5,6)])
+colnames(resED) <- paste("EC", c(10,20,50))
+knitr::kable(resED,caption = "Response Variable at day N",digits = 3)
+```
+**Calculate specific ECx:  **
+```{r}
+mod <-modList[[1]]
+edres <- ED.plus(mod,c(5,10,20,50),trend="Decrease")
+edres%>%knitr::kable(.,digits = 3)
+```
+
+
+## Model Output
+
+```{r}
+modsum <- summary(mod)
+knitr::kable(coef(modsum),digits = 3)
+```
+
+
+
+
 
 
 
-## Outputting
+## Additional Notes
 
 _To be written_
 
-## Dependencies
+### Dependencies
 
 ```r
 library(drcHelper)