---
title: "Appendix B: Cross-Method Validation of Scale-Adaptive Interaction Test Results via GAMM and Aligned Rank Transform (ART)"
author: "Cristobal Gallardo Alba <gallardoalba@pm.me>"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    toc_depth: 2
  rmarkdown::pdf_document:
    toc: true
    fig_width: 10
    fig_height: 6.5
    toc_depth: 2
    latex_engine: xelatex
header-includes:
  - \usepackage{float}
  - \usepackage{tabu}
bibliography: TSENAT.bib
suppress-bibliography: true
vignette: |
  %\VignetteIndexEntry{3. TSENAT Appendix B: Cross-Method Validation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r global-setup-appendix-b, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 10,
    fig.height = 6.5,
    message = FALSE,
    warning = FALSE,
    dev = "png",
    dpi = 300
)
suppressPackageStartupMessages(library(kableExtra))
```

# Introduction

The primary goal of this vignette is to validate that discoveries of scale-dependent entropic index × group interactions from Scale-Adaptive Interaction Tests (SAIT) generalize to non-parametric statistical frameworks with minimal assumptions.

## Two Complementary Validation Approaches

**Method 1: Generalized Additive Mixed Models (GAMM)** with ARIMA-Ordered Measurement Structure.

- Framework: Semi-parametric model combining flexible smooth functions (thin-plate splines) with mixed-effects structure; assumes additive model Y = f₁(q) + f₂(condition) + f₃(q, condition) + subject-intercept + ε where f terms are smooth rather than linear.
- Paired design structure: Uses random intercepts by subject (~1|subject) with AR(1) correlation structure to accommodate repeated entropy measurements across q-values within each subject. This hierarchical approach preserves paired sample structure while modeling autocorrelation in q-ordered measurements.
- Distributional assumption: Uses Gaussian residuals (standard for mixed-effects models); automatically detects bounded support [0,1] entropy but applies Gaussian family due to `mgcv::gamm()` limitations with extended families in paired designs.
- Treatment of entropic structure: Treats q-values as time-like ordered measurements within each subject, applies ARIMA(1,1,0) first-differencing to remove monotone trend and ensure stationarity; AR(1) correlation structure captured through nlme::corAR1() specification.
- Heteroscedasticity handling: Automatically detects variance heterogeneity (Breusch-Pagan test); however, heteroscedasticity-based variance weights are not applied in paired designs due to `mgcv::gamm()` limitations.
- Key advantage: Captures smooth nonlinear q × condition interactions while accommodating paired sample correlation structure through random intercepts and AR(1) correlation modeling.

**Method 2: Aligned Rank Transform (ART)** — State-of-the-art non-parametric interaction testing, with Hochberg Step-Up Correction. The Conover-Iman Rank Transform is also available via the `method='rt'` option as an alternative approach.

- Framework: Pure non-parametric method via the ARTool package [@S273]. Strips main effects before ranking ("alignment") to preserve interaction structure [@S271; @S272]. Zero parametric or distributional assumptions.
- Distributional assumption: None; operates entirely on aligned ranks. Valid under any continuous distribution and correlation structure.
- Treatment of q-ordering: Two-way ART on raw entropy values, treating q-values as categorical factors; aligned ranks decomposed via ANOVA with proper F-test for interactions. Because q is treated as a factor, the continuous ordering of q-values is discarded in the ART analysis — this is a fundamental difference from GAMM, which models q as a continuous smooth predictor.
- Heteroscedasticity handling: ART inherently robust to heteroscedasticity through alignment and ranking; Hochberg step-up correction controls FWER across multiple tests.
- Key advantage: State-of-the-art for non-parametric interaction testing. Properly handles factorial interactions by removing main effects before ranking. More reliable Type I error control than classical rank-transform methods for interaction terms.


## Why Two Methods for One Question?

Statistical testing in transcriptomics faces a fundamental challenge: no single method is universally optimal. Different approaches trade-off power for robustness:

| Aspect | Semi-Parametric (GAMM) | Non-Parametric (ART) |
|--------|------------------------|----------------------|
| Assumptions | Normality, homoscedasticity | None; aligned ranks, fully non-parametric |
| Power | Highest (if assumptions hold) | Good; reduced but robust to violations |
| Interaction testing | F-test on penalized smooth terms | Proper interaction tests via alignment [@S271] |
| q-value treatment | Continuous smooth predictor (preserves ordinal structure) | Categorical factor (discards ordering) |
| Robustness | Moderate | Highest |
| Outlier sensitivity | Moderate potential for bias | Minimal; inherently resistant |

Validation strategy: Concordance between both methods confirms robustness across statistical frameworks.

***

# Setup

This section initializes the analysis environment by loading required packages, setting a random seed for reproducibility, and preparing the test dataset.

```{r load-packages, message=FALSE}
suppressPackageStartupMessages({
    library(ggplot2)
    library(SummarizedExperiment)
    library(dplyr)
    library(gridExtra)
})

set.seed(42)

if (!requireNamespace("devtools", quietly = TRUE)) {
    stop("Please install devtools to run this vignette against the local package source.")
}
if ("package:TSENAT" %in% search()) {
    detach("package:TSENAT", unload = TRUE, character.only = TRUE)
}

pkg_root <- normalizePath("..", mustWork = FALSE)
devtools::load_all(pkg_root)

# Load preprocessed dataset
data(readcounts)
readcounts <- as.matrix(readcounts)
mode(readcounts) <- "numeric"

# Load metadata
metadata_df <- read.table(
    system.file("extdata", "metadata.tsv", package = "TSENAT"),
    header = TRUE, sep = "\t"
)

gff3_dataset <- system.file("extdata", "annotation.gff3.gz", package = "TSENAT")

# Configure analysis parameters first (best practice: fail-fast validation)
config <- TSENAT_config(
  sample_col = "sample",
  condition_col = "condition",
  subject_col = "paired_samples",
  q = seq(0, 2, by = 0.05),
  nthreads = 3,
  paired = TRUE,
  control = "normal"
)

# Build analysis object: creates SummarizedExperiment and initializes TSENATAnalysis
analysis <- build_analysis(
    readcounts = readcounts,
    tx2gene = gff3_dataset,
    metadata = metadata_df,
    tpm = tpm,
    effective_length = effective_length,
    config = config
)

# Apply filtering for quality control
analysis <- filter_analysis(
    analysis,
    stringency = "medium"
)
```

## Computing Tsallis Entropy with Pseudocount Regularization

Pseudocounts are critical for statistical robustness when applying non-parametric tests to sparse RNA-seq count data. RNA-seq experiments typically contain many zero or near-zero counts, which creates two problems for rank-based methods:

1. **Ties in rankings:** Sparse counts produce many tied values (especially zeros), which reduces the discriminatory power of rank tests. Rank-based statistics depend on unique orderings; when many observations are identical, the test statistic contains less information, reducing statistical power.

2. **Library size artifacts:** Small count differences due to sequencing depth variations overshadow true biological differences. Normalized library sizes (accounting for sequencing depth) ensure that entropy estimates reflect true biological diversity rather than technical artifacts [@S195; @S069].

By adding small pseudocounts proportional to library size, we achieve two benefits:

- **Regularization** breaks ties and prevents zero-inflation bias that compromises rank-based inference.
- **Normalization** makes entropy estimates comparable across samples with different sequencing depths. This approach is standard in RNA-seq analysis and is particularly important for entropy-based diversity metrics, which require positive values for logarithmic transforms.

The `calculate_diversity()` function implements automatic pseudocount selection, scaling pseudocount magnitude and diversity estimates by median library size to ensure robust statistical inference across the range of sequencing depths in your dataset.

```{r compute-tsallis-entropy, message=FALSE, results='hide'}
# Compute diversity using S4 wrapper with bootstrap confidence intervals
analysis <- calculate_diversity(
  analysis, 
  norm = TRUE,
  pseudocount = "auto"
)
```

# Aligned Rank Transform (ART)

TSENAT uses the **Aligned Rank Transform (ART)** as the default non-parametric method for testing Q×Condition interactions. ART is the state-of-the-art for non-parametric factorial analysis, implemented via the ARTool R package [@S273; @S272]. The classical Conover-Iman Rank Transform is also available via `method='rt'` as an alternative approach [@S270].

1. **Robustness to Distribution Violations**

Tsallis entropy values exhibit inherent distributional properties that often violate parametric assumptions:

- Non-normality: Entropy commonly exhibits bounded support (e.g., normalized entropy in [0,1]), producing skewed or bimodal distributions rather than normal distributions.
- Heteroscedasticity: Variance in entropy estimates varies systematically across q-values. Low q-values (emphasizing rare transcripts) produce volatile entropy estimates with high variance; high q-values (emphasizing abundant transcripts) produce stable estimates with lower variance.
- Outlier sensitivity: Rare isoforms and count variability can produce extreme entropy values that heavily influence parametric tests.

ART avoids these issues through alignment and ranking — stripping main effects before ranking preserves interaction structure while eliminating distributional dependence. The aligned ranks are uniformly distributed and contain no outliers, making the method valid for ANY continuous distribution [@S271; @S272]. Critically, rank ordering is unaffected by normalization choice; whether entropy is normalized, log-transformed, or left raw, the test produces identical results [@S270].

2. **Paired Design for Ordered q-Value Structure**

Tsallis entropy has a fundamental sequential property: entropy curves are smooth functions of q. As q increases from 0 (rare-species-emphasizing) to ∞ (common-species-emphasizing), entropy values change systematically. This ordered structure is critical to interpreting q-dependent patterns and exhibits **AR(1) autocorrelation**: consecutive q-values produce correlated entropy estimates.

ART handles this structure through its alignment step: within-subject ranking after stripping main effects preserves the pairing structure while removing the influence of AR(1) dependence [@S271; @S272]. Because aligned ranks are constructed from residuals after removing main effects, autocorrelation in the original entropy values does not propagate into the test statistic — rank-based tests are invariant to monotonic transformations of the data, and alignment renders the residuals approximately exchangeable under the null hypothesis of no interaction [@S019; @C018]. However, this robustness is not absolute: strong autocorrelation can reduce effective sample size and affect power even when the test remains valid.

3. **Interaction Testing: Testing q × Condition Effects**

The implementation uses ART (alignment + ranking + ANOVA), the modern state-of-the-art approach:

$$F_{q \times \text{condition}} = \frac{MS_{\text{interaction}}}{MS_{\text{residual}}}$$

where:

- Main effects are stripped via alignment: residuals = original − main effect estimates
- Aligned data are then ranked: $R = \text{rank}(\text{aligned residuals})$ (within-subject for paired designs)
- Two-way ANOVA decomposes aligned ranked data: $R = \mu + \alpha_q + \beta_{\text{condition}} + \gamma_{q \times \text{condition}} + \epsilon$
- $MS_{\text{interaction}}$ = sum of squares for $q \times$ condition interaction / degrees of freedom
- $MS_{\text{residual}}$ = residual sum of squares / residual degrees of freedom
- $F$-statistic follows $F$-distribution under null hypothesis of no interaction

This directly tests the core biological question: "Does entropy q-dependence differ between groups?" The aligned rank ANOVA is the modern state-of-the-art approach for non-parametric interaction testing, validated across balanced and unbalanced designs [@S272].

## Assumption Validation

To ensure valid statistical inference, we verify key data assumptions across multiple dimensions.

```{r validate-rank-transform-assumptions, message=FALSE, results='hide'}
# Validate statistical assumptions (including GAMM diagnostics)
analysis <- calculate_assumptions(
    analysis,
    checks = "all"  # Include core checks + GAMM diagnostics
)

# Get assumptions
assumptions_text <- results(analysis, type = "assumptions")
print(assumptions_text)
```

```{r display-assumption-results, message=FALSE, echo=FALSE}
# Extract assumption check results as structured list
check_results <- results(analysis, type = "assumptions", format = "list")

# Display assumptions table as formatted kable
if (!is.null(check_results) && !is.null(check_results$assumptions_table)) {
  assumptions_df <- check_results$assumptions_table
  
  knitr::kable(
    assumptions_df,
    caption = "**Supplementary Table 14 | Assumption Checks for Aligned Rank Transform (ART).** Evaluates homogeneity of variance and normality for rank-based interaction test validity."
  )
}
```


**Interpretation of Results:** The exchangeability test detects strong serial correlation (p ≈ 0), indicating that consecutive samples are more correlated than expected by chance. This violation is **not problematic** for the Aligned Rank Transform because ART's alignment step removes main effects before ranking, after which the aligned residuals are approximately exchangeable under the null hypothesis of no interaction [@S270; @S271]. While strong autocorrelation can reduce effective sample size and affect power, the test itself remains valid. The three reasons outlined above (distributional robustness, paired ordered structure handling, and proper interaction testing via alignment) together make ART ideally suited for multi-q entropy validation without requiring strong distributional assumptions.

## Aligned Rank Transform (ART): Testing q × Condition Interactions

Now we use the **Aligned Rank Transform (ART, default)** — the state-of-the-art non-parametric method for interaction testing — combined with Hochberg step-up correction for multiple testing across genes. ART properly handles factorial interactions by stripping main effects before ranking, providing more reliable Type I error control than classical rank-transform methods.

```{r rank-transform-test-analysis, message=FALSE, results='hide'}
# Run Aligned Rank Transform (ART, default method)
# This preserves the GAMM results in 'analysis' for later concordance comparison

# Run ART for q × condition interaction
analysis_rank <- calculate_rank_transform(
    analysis,
    multicorr = "hochberg"
)

# Extract ALL results (no filtering) for summary statistics
rank_transform_results_all <- results(analysis_rank, type = "rank_test", rankBy = "pvalue")

# Extract top significant genes for display
rank_transform_results <- results(analysis_rank, type = "rank_test", rankBy = "pvalue", n = 20, filterFDR = 0.05)
print(head(rank_transform_results, n = 10))
```

```{r rank-transform-results-summary, message=FALSE, echo=FALSE}
# Create summary statistics table with additional metrics using ALL tested genes
# Handle NA and NaN values properly
mean_effect <- mean(rank_transform_results_all$effect_size_eta2, na.rm = TRUE)
median_effect <- median(rank_transform_results_all$effect_size_eta2, na.rm = TRUE)

rank_summary_stats <- data.frame(
  Metric = c(
    "Genes tested", 
    "Significant (p < 0.05)", 
    "Significant (adj_p < 0.05, FWER-controlled)", 
    "NAs",
    "Mean effect size ($\\eta^2$)",
    "Median effect size ($\\eta^2$)",
    "Strong effect genes ($\\eta^2$ > 10%)"
  ),
  Value = c(
    nrow(rank_transform_results_all),
    sum(rank_transform_results_all$p_value < 0.05, na.rm = TRUE),
    sum(rank_transform_results_all$adj_p_value < 0.05, na.rm = TRUE),
    sum(is.na(rank_transform_results_all$p_value)),
    if (is.na(mean_effect) || !is.finite(mean_effect)) "Not available" else sprintf("%.1f%%", mean_effect * 100),
    if (is.na(median_effect) || !is.finite(median_effect)) "Not available" else sprintf("%.1f%%", median_effect * 100),
    sum(rank_transform_results_all$p_value < 0.05 & rank_transform_results_all$effect_size_eta2 > 0.10, na.rm = TRUE)
  )
)

knitr::kable(rank_summary_stats, caption = "**Supplementary Table 15 | Aligned Rank Transform (ART) Results Summary.** Non-parametric test of *q* × condition interactions via ART. Compare with GAMM results (main vignette Table 3) to assess parametric vs. non-parametric method concordance.", label = "rank-transform-test-summary", align = c("l", "r"))

# Show top genes results
rank_sig <- rank_transform_results_all[!is.na(rank_transform_results_all$p_value), ]
if (nrow(rank_sig) > 0) {
  rank_sig <- rank_sig[order(rank_sig$adj_p_value), ]
  rank_sig <- rank_sig[1:min(10, nrow(rank_sig)), ]
  
  # Select columns: gene, p-values, f-statistic, effect size, interaction class
  cols_to_show <- intersect(c("gene", "p_value", "adj_p_value", "f_statistic", "effect_size_eta2", 
                              "interaction_class"), 
                            colnames(rank_sig))
  if (length(cols_to_show) > 0) {
    rank_display <- rank_sig[, cols_to_show, drop = FALSE]
    
    # Format p-values with scientific notation for very small values
    if ("p_value" %in% colnames(rank_display)) {
      rank_display$p_value <- sapply(rank_display$p_value, function(x) {
        if (is.na(x)) "NA" else if (x >= 0.9999) "1" else if (x < 0.0001) format(signif(x, 2), scientific=TRUE) else sprintf("%.6f", x)
      })
    }
    
    if ("adj_p_value" %in% colnames(rank_display)) {
      rank_display$adj_p_value <- sapply(rank_display$adj_p_value, function(x) {
        if (is.na(x)) "NA" else if (x >= 0.9999) "1" else if (x < 0.0001) format(signif(x, 2), scientific=TRUE) else sprintf("%.6f", x)
      })
    }
    
    # Replace test method name for clarity (removed - column no longer included)
    # if ("test_method" %in% colnames(rank_display)) {
    #   rank_display$test_method <- sapply(rank_display$test_method, function(x) {
    #     if (is.na(x)) "NA" else if (x == "rt_paired") "Conover-Iman Rank Transform (paired)" else x
    #   })
    # }
    
    # Format other numeric columns
    if ("f_statistic" %in% colnames(rank_display)) {
      rank_display$f_statistic <- sapply(rank_display$f_statistic, function(x) {
        if (is.na(x)) "NA" else sprintf("%.4f", x)
      })
    }
    
    if ("effect_size_eta2" %in% colnames(rank_display)) {
      rank_display$effect_size_eta2 <- sapply(rank_display$effect_size_eta2, function(x) {
        if (is.na(x)) "NA" else sprintf("%.4f", x)
      })
    }
    
    knitr::kable(rank_display, 
                 caption = "**Supplementary Table 16 | Top genes identified by Aligned Rank Transform (ART).** Ranked by adjusted p-value; comparison with GAMM (main vignette Table 3) reveals method-specific sensitivities.",
                 col.names = c("Gene", "P-value", "Adj. P-value", "F-Statistic", 
                               "Effect Size (η²)", "Interaction Class")[1:ncol(rank_display)],
                 row.names = FALSE)
  }
}
```

Visualize q-curves for top genes:

```{r extended-fig-1-rank-transform-q-curves, fig.width=12, fig.height=8, message=FALSE, out.width="98%", fig.pos="H", fig.cap="**Supplementary Figure 1 | Aligned Rank Transform (ART) results for scale-dependent isoform switching.** Q-curves for top genes identified by ART; compare visual patterns with GAMM results (main vignette Figure 2) to assess method agreement."}
# Plot top genes from Aligned Rank Transform (ART) (using all genes, not just significant ones)
# This ensures the plot displays n_top genes regardless of significance threshold
top_genes_plot <- plot_diversity_spectrum(analysis, sait_res = rank_transform_results_all, n_top = 4)
print(top_genes_plot)
```

***

# Generalized Additive Mixed Model Approach (GAMM)

This section loads precomputed GAMM results generated in the main vignette ([TSENAT.Rmd](TSENAT.Rmd)).

```{r generalized-additive-model-analysis, message=FALSE, results='hide'}
# Load precomputed GAMM analysis object from RDS file
analysis_sait <- readRDS(
    system.file("extdata", "analysis_sait.rds", package = "TSENAT")
)

# Extract GAMM results for inspection using accessor function
sait_results <- results(analysis_sait, type = "sait", rankBy = "pvalue")
print(sait_results)

```

```{r gam-results-summary-table, message=FALSE, echo=FALSE}
# Create summary statistics table with additional metrics
# Handle NA and NaN values properly
mean_effect <- mean(sait_results$effect_size, na.rm = TRUE)
median_effect <- median(sait_results$effect_size, na.rm = TRUE)

summary_stats <- data.frame(
  Metric = c(
    "Genes tested", 
    "Significant (p < 0.05)", 
    "Concordant (p < 0.05 AND adj_p < 0.05)",
    "NAs",
    "Mean effect size",
    "Median effect size",
    "Strong effect genes (effect_size > 30%)",
    "Model convergence rate"
  ),
  Value = c(
    nrow(sait_results),
    sum(sait_results$p_interaction < 0.05, na.rm = TRUE),
    sum(sait_results$p_interaction < 0.05 & sait_results$adj_p_interaction < 0.05, na.rm = TRUE),
    sum(is.na(sait_results$p_interaction)),
    if (is.na(mean_effect) || !is.finite(mean_effect)) "Not available" else sprintf("%.1f%%", mean_effect * 100),
    if (is.na(median_effect) || !is.finite(median_effect)) "Not available" else sprintf("%.1f%%", median_effect * 100),
    sum(sait_results$p_interaction < 0.05 & sait_results$effect_size > 0.30, na.rm = TRUE),
    sprintf("%.1f%%", mean(sait_results$model_converged, na.rm = TRUE) * 100)
  )
)

knitr::kable(summary_stats, caption = "**Supplementary Table 17 | GAMM Results Summary.** Overview of generalized additive mixed model statistics for q × condition interaction tests across full q-spectrum with paired sample structure.", label = "gam-summary-stats", align = c("l", "r"))

# Show top genes with additional information
if ("p_interaction" %in% colnames(sait_results)) {
  gam_sig <- sait_results[!is.na(sait_results$p_interaction), ]
  
  if (nrow(gam_sig) > 0) {
    gam_sig <- gam_sig[order(gam_sig$p_interaction), ]
    gam_sig <- gam_sig[1:min(10, nrow(gam_sig)), ]
    
    # Select columns for display: gene, p-values, effect size, test statistic
    cols_to_show <- intersect(c("gene", "p_interaction", "adj_p_interaction", "effect_size", 
                                "test_statistic"), 
                              colnames(gam_sig))
    if (length(cols_to_show) > 0) {
      gam_display <- gam_sig[, cols_to_show, drop = FALSE]
      
      # Format numeric columns (but NOT p-value columns - preserve full precision for later formatting)
      numeric_cols <- sapply(gam_display, is.numeric)
      # Only round non-p-value columns to 3 decimal places
      pval_cols <- c("p_interaction", "adj_p_interaction")
      cols_to_round <- setdiff(names(numeric_cols[numeric_cols]), pval_cols)
      if (length(cols_to_round) > 0) {
        gam_display[, cols_to_round] <- round(gam_display[, cols_to_round], 3)
      }
      
      # Format p-value columns with scientific notation for very small values
      if ("p_interaction" %in% colnames(gam_display)) {
        gam_display$p_interaction <- sapply(gam_display$p_interaction, function(x) {
          if (is.na(x)) "NA" else if (x < 0.001) sprintf("%.2e", x) else sprintf("%.6f", x)
        })
      }
      
      if ("adj_p_interaction" %in% colnames(gam_display)) {
        gam_display$adj_p_interaction <- sapply(gam_display$adj_p_interaction, function(x) {
          if (is.na(x)) "NA" else if (x < 0.001) sprintf("%.2e", x) else sprintf("%.6f", x)
        })
      }
      
      # Format effect size as percentage
      if ("effect_size" %in% colnames(gam_display)) {
        gam_display$effect_size <- sapply(gam_display$effect_size, function(x) {
          if (is.na(x)) "NA" else sprintf("%.1f%%", x * 100)
        })
      }
      
      # Format test statistic
      if ("test_statistic" %in% colnames(gam_display)) {
        gam_display$test_statistic <- sapply(gam_display$test_statistic, function(x) {
          if (is.na(x)) "NA" else sprintf("%.2f", x)
        })
      }
      
      # Rename columns for better readability — dynamic to match actual columns present
      col_name_map <- c(
        "gene" = "Gene",
        "p_interaction" = "p (interaction)",
        "adj_p_interaction" = "Adjusted p",
        "effect_size" = "Effect size",
        "test_statistic" = "Test statistic"
      )
      present_cols <- intersect(names(col_name_map), colnames(gam_display))
      gam_display <- gam_display[, present_cols, drop = FALSE]
      colnames(gam_display) <- unname(col_name_map[present_cols])
      n_cols <- ncol(gam_display)
      
      knitr::kable(gam_display, 
                   caption = "**Supplementary Table 18 | Top 10 Genes by GAMM p-value (with Effect Size and Test Statistics).** Ranked by adjusted p-value; compare with ART results (Supplementary Table 16) for method agreement assessment.",
                   label = "gam-top-genes",
                   align = c("l", rep("r", max(0, n_cols - 1))),
                   row.names = FALSE)
    }
  }
} else {
  cat("WARNING: p_interaction column not found in results\n")
  cat("Available columns:", paste(colnames(sait_results), collapse=", "), "\n")
}
```

***

# Method Concordance Analysis

A critical validation step is to compare whether the Aligned Rank Transform (ART) and GAMM produce consistent results. High concordance between methods (i.e., genes significant in both approaches) provides strong evidence that discoveries are robust to methodological choice. Discordant genes — those significant in only one method — warrant closer inspection: they may represent genuine biological signals revealed by one method's specific advantages, or artifacts of that method's assumptions. This section quantifies agreement between approaches and identifies high-confidence genes significant across both statistical frameworks.

```{r compare-rank-transform-gam-methods, message=FALSE}
# Compute concordance analysis using new two-object API
rank_results <- results(analysis_rank, type = "rank_test")
if (is.null(rank_results) || length(rank_results) == 0) {
    stop("Cannot run calculate_concordance(): rank_test results are empty. Ensure calculate_rank_transform() completed successfully.")
}

analysis_rank <- calculate_concordance(
    analysis_sait = analysis_sait,
    analysis_rank = analysis_rank,
    verbose = TRUE
)

# Extract concordance results with list format for component display
concordance_results <- results(analysis_rank, type = "concordance", format = "list")
```

```{r display-concordance-metrics, message=FALSE, echo=FALSE, results='asis'}
# Display concordance summary - results() with format="list" returns structured components
if (!is.null(concordance_results) && is.list(concordance_results)) {
  # Display summary table
  if (!is.null(concordance_results$summary_table) && nrow(concordance_results$summary_table) > 0) {
    print(knitr::kable(concordance_results$summary_table, 
                 caption = "**Supplementary Table 19 | Global Concordance Metrics: GAMM vs Aligned Rank Transform (ART).** Quantifies agreement between parametric (GAMM) and non-parametric (ART) approaches for detecting q × condition interactions.",
                 align = c("l", "r"),
                 row.names = FALSE))
  }
  
  # Display agreement distribution
  if (!is.null(concordance_results$agreement_dist) && nrow(concordance_results$agreement_dist) > 0) {
    cat("\n")
    print(knitr::kable(concordance_results$agreement_dist, 
                 caption = "**Supplementary Table 20 | Agreement Distribution: Concordance by Gene Category.** Distribution of genes across categories: both methods significant, GAMM only, ART only, and neither method.",
                 align = c("l", "r", "r"),
                 row.names = FALSE))
  }
  
  # Display high-confidence genes (concordant in both methods)
  if (!is.null(concordance_results$high_conf_table) && nrow(concordance_results$high_conf_table) > 0) {
    cat("\n")
    n_hc <- nrow(concordance_results$high_conf_table)
    print(knitr::kable(concordance_results$high_conf_table, 
                 caption = paste0("**Supplementary Table 21 | Robust Entropic Order Index Interactions: High-Confidence Genes Detected by Both Methods (n=", n_hc, ", ranked by statistical significance).** Cross-method validation: genes significant in both GAMM and ART provide strongest evidence for q × condition effects."),
                 align = c("l", "r", "r", "r", "r"),
                 row.names = FALSE))
  }
} else {
  cat("Concordance results not available\n")
}
```

## Statistical Power vs. Robustness: Understanding Method Discordance

The substantial discordance between GAMM and ART results — where GAMM identifies many more significant genes — reflects a fundamental trade-off in statistical methodology: **parametric methods maximize power when assumptions hold, while non-parametric methods sacrifice power for robustness to assumption violations**. The exact percentages depend on the dataset and analysis parameters; the pattern of GAMM detecting more interactions than ART is the consistent finding.

GAMM's superior power derives from three factors:

1. Distributional assumptions: GAMM assumes approximately normal residuals and homogeneous variance, which are reasonable after appropriate transformation for entropy data. Under these assumptions, parametric methods are theoretically optimal, achieving maximum power for a given type I error rate. ART is assumption-free but necessarily discards quantitative information through alignment and ranking, which reduces statistical power when the underlying data are approximately normal.

2. Model flexibility with penalty: GAMM uses thin-plate splines with smoothness penalties that simultaneously fit nonlinear patterns while controlling degrees of freedom. This flexibility allows GAMM to detect subtle entropic index × group interactions across all q-values.

3. Continuous vs. categorical q treatment: GAMM models q as a continuous smooth predictor across all 41 q-values, preserving the ordinal structure and borrowing strength across adjacent q-levels. ART treats q as a categorical factor — while this enables proper interaction testing via alignment, it discards the continuous ordering information, reducing sensitivity to graded q-dependent trends. This discretization is a major contributor to ART's lower power relative to GAMM, independent of distributional assumptions.

This complementary approach—combining high-power parametric tests with robust nonparametric alternatives—provides confidence that discoveries are methodologically sound rather than artifacts of statistical assumptions.

## Visualization: Method Comparison

Visual comparison of concordance patterns provides intuitive assessment of which genes align between methods and which are method-specific. The following plots display significance calls, effect sizes, and agreement metrics in a format that facilitates interpretation of both robust discoveries and method-specific findings.

```{r visualize-method-concordance, fig.width=10, fig.height=8, out.width="98%", fig.pos="H", fig.cap="**Supplementary Figure 2 | Method concordance visualization comparing Aligned Rank Transform (ART) and GAMM approaches for detecting q × condition interactions.** Scatter plot showing overlap in genes identified as significant by each method; reveals parametric-non-parametric agreement patterns."}
# Create comparison plot using S4 wrapper
# Concordance analysis results are stored in analysis metadata after calculate_concordance()
p <- plot_concordance(analysis_rank, verbose = TRUE)
grid::grid.draw(p)
```

***

# Session Information

```{r session-info-appendix-b}
sessionInfo()
```
