---
title: 'Appendix S9: Codes for used statistical analyses - ULTRASOUND'
output:
  html_document:
    toc: true
    toc_depth: '2'
    df_print: paged
  html_notebook:
    toc: true
    toc_depth: 2
  pdf_document:
    toc: true
    toc_depth: '2'
---

### Appendix S10: Codes for used statistical analyses

**Title:** *A workflow to optimize spatial sampling in ecoacoustic studies*\
**Journal:** *Landscape Ecology*

**Authors:**\
Martínez-Arias, V. M.\textsuperscript{1}; Paniagua-Villada, C.\textsuperscript{1}; Daza, J. M.\textsuperscript{1}

**Author affiliations at the time the work was conducted:**\
1. Grupo Herpetológico de Antioquia, Calle 70 No. 52 - 21, Postal Code 050010, Medellín, Colombia.

**Corresponding author:**\
Victor M. Martínez-Arias\
📧 Email: [vmanuel.martinez\@udea.edu.co](mailto:vmanuel.martinez@udea.edu.co){.email}\


# 0. Exploration graphic

This graphic allow to visually compare the indices in the range 0-24khz, for all hours and localities.

```{r}
library(ggplot2)
library(readxl)
library(dplyr)

Statistic <- "mean"  # Set to "mean" or "median"
# Load data (replace path with your own)
data <- read_excel("#/DF1-1_Interpolations_variograms_and_error_for_subsamples.xlsx")
# Filter based on combined conditions
filtered_data <- data %>%
  filter(
    (Index %in% c("BI", "NDSI") & Statistic == Statistic) |
    (Index %in% c("ACItf", "ACIft", "NP") &
     Frequency_Range == "R0-24" & Statistic == "mean")
  ) %>%
  mutate(
    Subsample_Color = case_when(
      Subsample == "SLIC" ~ "SLIC",
      Subsample == "halfbasins" ~ "HalfBasins",
      Subsample == "basins" ~ "Basins",
      Subsample == "Landcover" ~ "Landcover",
      TRUE ~ "Random"
    )
  )
# Add synthetic category "ALL"
data_all <- filtered_data %>%
  mutate(Locality = "ALL")
# Merge both datasets
complete_data <- bind_rows(filtered_data, data_all)
# Create plot
plot <- ggplot(complete_data, aes(x = Sample_percent, y = Pearson_Correlation)) +
  geom_point(data = filter(complete_data, Subsample_Color == "Random"),
             aes(color = Subsample_Color),
             alpha = 0.25, size = 1) +
  geom_smooth(method = "lm", color = "blue", se = FALSE, linewidth = 0.8) + 
  geom_point(data = filter(complete_data, Subsample_Color != "Random"),
             aes(color = Subsample_Color), size = 2) +
  scale_color_manual(
    values = c(
      "SLIC" = "purple",
      "HalfBasins" = "skyblue",
      "Basins" = "navy",
      "Random" = "gray50",
      "Landcover" = "green"
    )
  ) +
  facet_grid(Locality ~ Index) +
  theme_minimal(base_size = 12) +
  labs(
    x = "Sample percent",
    y = "Subsample_performance",
    color = "Sampling method"
  ) +
  theme(
    legend.position = "right",
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12, face = "bold")
  )
# Save as PDF (set your desired path)
ggsave(
  filename = "#/0_Subsampleperformance_vs_SamplePercent_mean.pdf",
  plot = plot,
  width = 15.27,
  height = 8.15,
  units = "in"
)

print(plot)



```



# 0.1 Subsampling percent of the initial sample
This section summarizes the effective sampling effort associated with each spatial subsampling strategy, expressed as the percentage of recorder locations retained relative to the full sampling grid deployed at each locality.

We use a unified database that records, for every combination of locality, acoustic index, frequency band, and time period, the number of recorder locations selected under each sampling strategy (Basins, Halfbasins, Landcover, SLIC, and random subsamples at 20%, 40%, 60%, and 80%). Because recordings affected by rain or corrupted files were removed prior to temporal aggregation, the effective number of recorder locations may vary slightly across frequency bands and time periods. These variations are minimal and reflect quality-control procedures rather than differences in spatial design. To facilitate comparison across strategies and localities, the number of retained recorder locations is normalized by the total number of available locations per locality, yielding a percentage-based measure of sampling effort. This normalization allows direct assessment of how efficiently each strategy represents the landscape relative to the full grid.

The first output of this section is a faceted plot summarizing the distribution of sampling effort across strategies and localities. For each strategy, boxplots represent variation across indices, frequency bands, and time periods; points indicate median values, and vertical lines show the minimum–maximum range. This visualization supports direct comparison of how many recorder locations are required by each strategy to achieve comparable spatial representation.

The second output is a summary table reporting the mean and standard deviation of the percentage of recorder locations retained by each strategy in each locality. This table complements the graphical summary by providing concise numerical descriptors of sampling efficiency and stability, and is referenced in the Supplementary Information.

Graphic:
```{r}
# ===============================
# Libraries
# ===============================
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)

# ===============================
# Load data
# ===============================
df <- read_excel("~/Dropbox/Doctorate project - VMMA/Thesis-associated_papers/P001_Sampling_design/Online_Resoures-SUPPLEMENTARY_INFORMATION/Supplementary_information_S12_Associated_Databases/Sampling_count_database.xlsx")

# ===============================
# Wide → Long
# ===============================
df_long <- df %>%
  pivot_longer(
    cols = c(
      basins,
      halfbasins,
      Landcover,
      SLIC,
      ssample20,
      ssample40,
      ssample60,
      ssample80
    ),
    names_to = "Strategy",
    values_to = "Num_Points"
  )

# ===============================
# Normalize to percentage
# ===============================
df_long <- df_long %>%
  mutate(
    Percent_Used = (Num_Points / Num_Points_Total) * 100
  )

# ===============================
# Strategy order
# ===============================
df_long$Strategy <- factor(
  df_long$Strategy,
  levels = c(
    "basins",
    "halfbasins",
    "Landcover",
    "SLIC",
    "ssample20",
    "ssample40",
    "ssample60",
    "ssample80"
  )
)

# ===============================
# Plot (FACETED)
# ===============================
p <- ggplot(df_long, aes(x = Locality, y = Percent_Used)) +

  # Rango min–max (muy sutil)
  stat_summary(
    fun.min = min,
    fun.max = max,
    geom = "linerange",
    linewidth = 0.3,
    color = "grey40"
  ) +

  # Boxplot fino
  geom_boxplot(
    width = 0.18,
    color = "black",
    fill = "white",
    linewidth = 0.35,
    outlier.shape = NA
  ) +

  # Mediana
  stat_summary(
    fun = median,
    geom = "point",
    size = 1.6,
    color = "black"
  ) +

  facet_wrap(~ Strategy, ncol = 2) +

  labs(
    x = "Locality",
    y = "Recorder locations used (% of full grid)"
  ) +

  theme_classic(base_family = "Times New Roman") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 11),
    axis.title = element_text(size = 13),
    panel.spacing = unit(0.8, "lines")
  )

print(p)

ggsave(
  "Subsampling_efficiency_box_range.png",
  p,
  width = 8,
  height = 10,
  dpi = 300
)

```

Table:
```{r}
# ===============================
# Libraries
# ===============================
library(readxl)
library(dplyr)
library(tidyr)
library(openxlsx)

# ===============================
# Load data
# ===============================
df <- read_excel(
  "~/Dropbox/Doctorate project - VMMA/Thesis-associated_papers/P001_Sampling_design/Online_Resoures-SUPPLEMENTARY_INFORMATION/Supplementary_information_S12_Associated_Databases/Sampling_count_database.xlsx"
)

# ===============================
# Wide → Long
# ===============================
df_long <- df %>%
  pivot_longer(
    cols = c(
      basins,
      halfbasins,
      Landcover,
      SLIC,
      ssample20,
      ssample40,
      ssample60,
      ssample80
    ),
    names_to = "Strategy",
    values_to = "Num_Points"
  )

# ===============================
# Normalize to percentage
# ===============================
df_long <- df_long %>%
  mutate(
    Percent_Used = (Num_Points / Num_Points_Total) * 100
  )

# ===============================
# Summary table: mean & SD
# ===============================
summary_table <- df_long %>%
  group_by(Locality, Strategy) %>%
  summarise(
    Mean_Percent = mean(Percent_Used, na.rm = TRUE),
    SD_Percent   = sd(Percent_Used, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(Strategy, Locality)

# ===============================
# Export to Excel
# ===============================
output_path <- "/home/vmartinezarias/Dropbox/Doctorate project - VMMA/Thesis-associated_papers/P001_Sampling_design/Online_Resoures-SUPPLEMENTARY_INFORMATION/Supplementary_information_S12_Associated_Databases/Sampling_summary_mean_SD.xlsx"

write.xlsx(
  summary_table,
  file = output_path,
  overwrite = TRUE
)

cat("✅ Summary table exported to:", output_path, "\n")

```




# OBJECTIVE 1 SIZE DEPENDANCE

S. objective 1 - Analyze the relationship between sample size and the spatial representativeness of interpolated soundscapes using different sampling strategies in heterogeneous landscapes;

## 1.1. Generate the pearson correlation between Sampling percent and design precision.

To evaluate how sampling size influences the performance of each subsampling strategy, we calculated a second Pearson correlation coefficient, Pcoef2. This coefficient quantifies the correlation between the proportion of sampled data (i.e., Sample_percent) and the subsampling performance, previously calculated as Pcoef1 (Pearson correlation between full-sample and subsampled raster surfaces).

The analysis was performed across all combinations of localities, subsampling designs, frequency ranges, time periods, acoustic indices, and summary statistics (e.g., mean, median). For each unique combination, we filtered the data and applied ggscatterstats() from the ggstatsplot package to compute the Pearson correlation (Pcoef2), along with associated test statistics.

The resulting outputs were consolidated into a dataset named DF2, which was later used for statistical comparisons between strategies and visualized using violin and density plots.

Used filters: \| Locality \| Type \| Frequency_Range \| Time_Period \| Statistic \| Index \| Pearson_Correlation \|

Resulting table will be named "DF2", which includes the calculated Pcoef2.

```{r}
# =============================================================================
# SCRIPT: GENERATE DF2 (CORRELATIONS) - ALL FREQUENCIES
# =============================================================================

# -----------------------------------------------------------------------------
# 1. SETUP AND LIBRARIES
# -----------------------------------------------------------------------------
# Ensure pacman is installed to manage libraries easily
if (!require("pacman")) install.packages("pacman")

# Load necessary libraries
pacman::p_load(ggstatsplot, dplyr, readxl, openxlsx, future, future.apply, stringr)

# Configure Parallel Processing
# We use (Total Cores - 1) to keep the system responsive
n_cores <- parallel::detectCores(logical = FALSE) - 1
plan(multisession, workers = n_cores)

message(paste("✅ Parallel processing configured with", n_cores, "cores."))

# -----------------------------------------------------------------------------
# 2. LOAD DATA
# -----------------------------------------------------------------------------
# Input path definition
input_path <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF1-1_Interpolations_variograms_error_and_performance_of_subsamples_US.xlsx"

message("📂 Loading data from sheet 'DATA'...")

# Load the specific sheet and ensure key columns are characters
data <- read_excel(input_path, sheet = "DATA") %>%
  mutate(across(c(Locality, Type, Subsample, Frequency_Range, Time_Period, Statistic, Index), as.character))

# -----------------------------------------------------------------------------
# 3. AUDIO CATEGORY CLASSIFICATION (CRITICAL FIX)
# -----------------------------------------------------------------------------
# This section ensures that both Audible and Ultrasound are identified correctly
# based on the ranges you provided.

data <- data %>%
  mutate(Audio_Category = case_when(
    # Ultrasound Ranges
    Frequency_Range %in% c("R24-96", "R12-24") ~ "Ultrasound",
    
    # Audible Ranges
    Frequency_Range %in% c("R0-2", "R2-12", "R0-24") ~ "Audible",
    
    # Fallback
    TRUE ~ "Unknown" 
  ))

# Verify that we have both categories
message("📊 Data distribution by Audio Category:")
print(table(data$Frequency_Range, data$Audio_Category))

# -----------------------------------------------------------------------------
# 4. DATA EXPANSION (GLOBAL ANALYSIS)
# -----------------------------------------------------------------------------
# Duplicate the dataset, labeling the copy as Locality = "ALL" 
# to calculate global correlations later.
data_all <- data %>% mutate(Locality = "ALL")
expanded_data <- bind_rows(data, data_all)

# Identify unique combinations to process
# We include 'Audio_Category' here to ensure it is passed to the final table
combinations <- expanded_data %>%
  distinct(Locality, Type, Audio_Category, Frequency_Range, Time_Period, Statistic, Index, Subsample)

message(paste("🔄 Total combinations to process:", nrow(combinations)))

# -----------------------------------------------------------------------------
# 5. PROCESSING FUNCTION
# -----------------------------------------------------------------------------
process_combination <- function(filter_row) {
  
  # Filter the main dataset for the current combination
  filtered_data <- expanded_data %>%
    filter(
      Locality == filter_row$Locality,
      Type == filter_row$Type,
      Frequency_Range == filter_row$Frequency_Range,
      Time_Period == filter_row$Time_Period,
      Statistic == filter_row$Statistic,
      Index == filter_row$Index,
      Subsample == filter_row$Subsample
    )
  
  # --- Data Sufficiency Checks ---
  # 1. Minimum 3 data points required for correlation
  # 2. Must have variation in X (Sample_percent)
  # 3. Must have variation in Y (Subsample_performance)
  if (nrow(filtered_data) < 3 ||
      length(unique(filtered_data$Sample_percent)) < 2 ||
      length(unique(filtered_data$Subsample_performance)) < 2) {
    return(NULL)
  }
  
  # --- Calculate Statistics (ggscatterstats) ---
  # Wrapped in tryCatch to prevent the loop from breaking on a single error
  temp_plot <- tryCatch({
    ggstatsplot::ggscatterstats(
      data = filtered_data,
      x = Sample_percent,
      y = Subsample_performance,
      type = "parametric", # Pearson Correlation
      messages = FALSE
    )
  }, error = function(e) return(NULL))
  
  if (is.null(temp_plot)) return(NULL)
  
  # Extract the statistical table from the plot object
  stats <- tryCatch(ggstatsplot::extract_stats(temp_plot), error = function(e) return(NULL))
  if (is.null(stats)) return(NULL)
  
  subtitle_data <- stats$subtitle_data
  caption_data <- stats$caption_data
  
  # Add the identification columns (Metadata) to the results
  for (col in names(filter_row)) {
    subtitle_data[[col]] <- filter_row[[col]]
    caption_data[[col]] <- filter_row[[col]]
  }
  
  # --- FIX: PREVENT OVERWRITING 'TYPE' ---
  # Previous code overwrote 'Type' (Strategy). 
  # Here we create a NEW column for the paradigm.
  subtitle_data$Stat_Paradigm <- "Frequentist"
  caption_data$Stat_Paradigm <- "Bayesian"
  
  bind_rows(subtitle_data, caption_data)
}

# -----------------------------------------------------------------------------
# 6. PARALLEL EXECUTION
# -----------------------------------------------------------------------------
message("🚀 Starting parallel processing...")

# Run the function across all combinations
results_list <- future_lapply(1:nrow(combinations), function(i) {
  filter_row <- combinations[i, ]
  process_combination(filter_row)
}, future.seed = TRUE)

# Combine the list of dataframes into one single dataframe
final_table <- bind_rows(results_list)

# -----------------------------------------------------------------------------
# 7. EXPORT RESULTS
# -----------------------------------------------------------------------------
if (nrow(final_table) > 0) {
  
  # Remove list-type columns that cause issues in Excel
  final_table_clean <- final_table %>%
    select(where(~ !is.list(.)))
  
  # Define output path
  output_path <- "/home/vmartinezarias/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF2_SampleSize_vs_Subsample_performance_Correlations(Pcoeff2)_COMPLETE.xlsx"
  
  message("💾 Saving file...")
  write.xlsx(final_table_clean, output_path)
  
  message("✅ SUCCESS! File generated at: ", output_path)
  message("   Total rows generated: ", nrow(final_table_clean))
  
} else {
  warning("⚠️ The final table is empty. Check if the input filters match the data.")
}



```

## 1.2. Visual inspection of sampling size effects

To explore the relationship between sample size and subsampling performance across different categorical variables (e.g., locality, index, strategy), we created an interactive visualization tool using a Shiny application. This tool allows the dynamic generation of violin plots showing the distribution of Pcoef2 values (correlation between sample size and subsample performance) across user-defined categories.

The application permits filtering by multiple metadata fields (e.g., frequency range, time period, statistic, or acoustic index) and dynamically updates the visualization accordingly. Violin plots were generated using the ggplot2 package, with optional overlay of jittered data points. Users can export the resulting plots as PDFs for reporting or supplementary use.

This tool is based on the table DF2, and facilitates visual comparison of the influence of sample size across grouping levels in the dataset.

This visual inspection is only for ALL dataset.


```{r}
library(shiny)
library(ggplot2)
library(readxl)
library(dplyr)

# ==============================================================================
# UI DEFINITION
# ==============================================================================
ui <- fluidPage(
  titlePanel("Distribution of Pearson Coefficient (Publication Ready)"),
  
  sidebarLayout(
    sidebarPanel(
      # File upload
      fileInput("file", "Upload Excel file", accept = c(".xlsx")),
      
      # Main variable
      selectInput("categorical_var", "Main categorical variable (X-axis):",
                  choices = c("Locality", "Frequency_Range", "Time_Period", 
                              "Statistic", "Index", "Subsample", "Type")),
      
      # Checkbox to toggle points
      checkboxInput("show_points", "Show individual points (Jitter)", value = TRUE),
      
      # Checkbox to split by Audio Type
      checkboxInput("split_audio", "Split by Audio Type (Audible vs Ultrasound)", value = FALSE),
      
      hr(), 
      
      h4("Dynamic Filters"),
      p("Select 'All' to include everything or filter by specific groups:"),
      uiOutput("additional_filters"),
      
      hr(),
      
      h4("Paper Formatting (PDF)"),
      # Controls for font size and dimensions
      numericInput("axis_title_size", "Axis title font size:", value = 14, min = 8, max = 30),
      numericInput("axis_text_size", "Axis text font size:", value = 12, min = 6, max = 30),
      
      splitLayout(
        numericInput("pdf_width", "Width (mm):", value = 170, min = 50, max = 300),
        numericInput("pdf_height", "Height (mm):", value = 110, min = 50, max = 500)
      ),
      helpText("Max width usually 170mm. Adjust height as needed."),
      
      # Download button
      downloadButton("download_pdf", "Download PDF for Paper")
    ),
    
    mainPanel(
      tabsetPanel(
        tabPanel("Plot Preview", plotOutput("violin_plot_output", height = "600px")),
        tabPanel("Data Table", tableOutput("filtered_table_output"))
      )
    )
  )
)

# ==============================================================================
# SERVER LOGIC
# ==============================================================================
server <- function(input, output, session) {
  
# 1. Load Data
  raw_data <- reactive({
    req(input$file)
    
    # --- MODIFICACIÓN AQUÍ ---
    # Se especifica que cargue la hoja "DATA"
    df <- read_excel(input$file$datapath, sheet = "DATA")
    
    # Check for correct column and filter
    if("Statistical_paradigm" %in% names(df)){
      df <- df %>% filter(Statistical_paradigm == "Frequentist")
    }
    
    # Create Audio_Type column based on Frequency_Range
    if("Frequency_Range" %in% names(df)){
      df <- df %>% mutate(
        Audio_Type = ifelse(grepl("24-96", Frequency_Range), "Ultrasound", "Audible")
      )
    }
    
    df %>% filter(!is.na(estimate))
  })
  
  # 2. Dynamic UI for Filters
  output$additional_filters <- renderUI({
    req(raw_data())
    all_categories <- c("Locality", "Frequency_Range", "Time_Period", 
                        "Statistic", "Index", "Subsample", "Type")
    filter_candidates <- setdiff(all_categories, input$categorical_var)
    
    lapply(filter_candidates, function(col) {
      if(col %in% names(raw_data())) {
        options_list <- c("All", sort(unique(as.character(raw_data()[[col]]))))
        selectInput(inputId = paste0("filter_", col),
                    label = paste("Filter by", col),
                    choices = options_list, selected = "All")
      }
    })
  })
  
  # 3. Filter Logic
  filtered_data <- reactive({
    req(raw_data())
    df <- raw_data()
    all_categories <- c("Locality", "Frequency_Range", "Time_Period", 
                        "Statistic", "Index", "Subsample", "Type")
    active_filters <- setdiff(all_categories, input$categorical_var)
    
    for (col in active_filters) {
      input_id <- paste0("filter_", col)
      selected_val <- input[[input_id]]
      if (!is.null(selected_val) && selected_val != "All") {
        if(col %in% names(df)){
          df <- df %>% filter(.data[[col]] == selected_val)
        }
      }
    }
    df
  })
  
  # 4. Create Plot Object
  create_plot <- reactive({
    req(filtered_data())
    df <- filtered_data()
    validate(need(nrow(df) > 0, "No data available with the current filters."))
    
    # Base Plot with Colors
    p <- ggplot(df, aes(x = .data[[input$categorical_var]], y = estimate, fill = Type)) +
      
      # VIOLIN: Added lwd=0.3 to make border thinner
      geom_violin(alpha = 0.5, color = "darkblue", trim = TRUE, 
                  position = position_dodge(width = 0.9), lwd = 0.3) +
      
      # BOXPLOT: 
      # 1. lwd = 0.3 (thinner lines)
      # 2. outlier.shape = NA (removes the "fat dots" you circled, since jitter shows them anyway)
      geom_boxplot(width = 0.1, color = "black", alpha = 0.8, 
                   position = position_dodge(width = 0.9), show.legend = FALSE, 
                   lwd = 0.3, outlier.shape = NA) +
      
      scale_fill_manual(values = c("Random" = "lightpink", "Tessells" = "skyblue"))
    
    # Conditionally add points (Jitter)
    if (input$show_points) {
      p <- p + geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.9),
                          alpha = 0.6, size = 1.0, color = "darkslategray", stroke = 0) 
                          # Reduced size slightly and removed stroke for cleaner look
    }
    
    # Split by Audio Type SIDE-BY-SIDE
    if (input$split_audio) {
      p <- p + facet_wrap(~ Audio_Type, nrow = 1) 
    }
    
    # Styling and Text
    p <- p + 
      scale_y_continuous(limits = c(-1, 1), breaks = seq(-1, 1, 0.2)) + 
      labs(
        x = input$categorical_var, 
        y = "Pearson Correlation Estimate",
        fill = "Method",
        title = NULL 
      ) +
      theme_minimal() +
      theme(
        text = element_text(family = "serif"), 
        axis.title = element_text(size = input$axis_title_size, face = "bold"),
        axis.text = element_text(size = input$axis_text_size, color = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid.minor = element_blank(),
        legend.position = "top",
        strip.text = element_text(size = input$axis_text_size, face = "bold", family = "serif")
      )
      
    p
  })
  
  # 5. Output to Screen
  output$violin_plot_output <- renderPlot({
    create_plot()
  })
  
  # 6. Output Table
  output$filtered_table_output <- renderTable({
    head(filtered_data(), 20) 
  })
  
  # 7. PDF Download Handler
  output$download_pdf <- downloadHandler(
    filename = function() {
      paste0("Pearson_Plot_", Sys.Date(), ".pdf")
    },
    content = function(file) {
      ggsave(file, 
             plot = create_plot(), 
             device = "pdf", 
             width = input$pdf_width,  
             height = input$pdf_height, 
             units = "mm") 
    }
  )
}

shinyApp(ui = ui, server = server)
```


## 1.3. Normality assessment

Before a statistical test to evaluate the size dependance, normality must be checked. To determine whether parametric tests could be applied to compare the sampling strategies, we evaluated the distribution of Pcoef2 values (i.e., correlations between sampling proportion and subsample performance) across all sampling designs. We applied the Shapiro–Wilk test for normality within each subsample category using the rstatix package. Additionally, we visually inspected the distributions by plotting histograms with density curves, faceted by subsample strategy (e.g., basins, SLIC, ssample40).

Results from the Shapiro–Wilk tests and visual assessments revealed widespread deviations from normality, justifying the use of non-parametric statistical tests in subsequent analyses. All results were saved as part of the reproducible workflow, including the test results and distribution plots.

The input is DF2_SampleSize_vs_Subsample_performance_Correlations(Pcoeff2).xlsx. Procedure must be done for mean and median

```{r}
# ==============================================================================
# SETUP AND LIBRARIES
# ==============================================================================

# Load required packages
library(readxl)
library(dplyr)
library(ggplot2)
library(rstatix)   # For normality tests and summaries
library(ggpubr)    # For visualization
library(writexl)   # For exporting Excel files

# Define input and output paths (Update these paths as needed)
input_path <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF2_SampleSize_vs_Subsample_performance_Correlations(Pcoeff2)_COMPLETE.xlsx"
output_dir <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/" # Ensure this directory exists

# Define the statistic you want to analyze ("mean", "median", etc.)
target_statistic <- "median" 

# ==============================================================================
# DATA LOADING AND PROCESSING
# ==============================================================================

# Load dataset
data <- read_excel(input_path, sheet='DATA')

# Filter data
# Note: Based on your sample, 'Statistical_paradigm' contains "Frequentist", 
# while 'Type' contains values like "Tessells" or "Random".
data_frec <- data %>%
  filter(
    Statistical_paradigm == "Frequentist",  # Corrected column name based on your data snippet
    !is.na(estimate),
    Statistic == target_statistic           # Filters dynamically based on the variable above
  )

# Check if data exists after filtering
if(nrow(data_frec) == 0) {
  stop("Error: The filtered dataset is empty. Check the spelling of 'Frequentist' or the target statistic.")
}

# ==============================================================================
# STATISTICAL ANALYSIS (SHAPIRO-WILK)
# ==============================================================================

# Perform Shapiro-Wilk normality test by Subsample
# Note: Shapiro-Wilk requires at least 3 observations per group.
normality <- data_frec %>%
  group_by(Subsample) %>%
  filter(n() >= 3) %>%  # Safety check: exclude groups with insufficient data
  shapiro_test(estimate)

print("### Shapiro-Wilk normality test by Subsample:")
print(normality)

# ==============================================================================
# VISUALIZATION
# ==============================================================================

# Create histogram and density plots
plot <- ggplot(data_frec, aes(x = estimate, fill = Subsample)) +
  # Updated '..density..' to 'after_stat(density)' for newer ggplot versions
  geom_histogram(aes(y = after_stat(density)), position = "identity", alpha = 0.5, bins = 30) +
  geom_density(alpha = 0.7) +
  facet_wrap(~Subsample) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "right",
    axis.text = element_text(size = 10),
    strip.text = element_text(size = 12, face = "bold")
  ) +
  labs(
    title = paste("Distribution of Pearson Estimates (Statistic:", target_statistic, ")"),
    x = "Pearson Estimate", 
    y = "Density"
  )

# Display plot in RStudio
print(plot)

# ==============================================================================
# EXPORT RESULTS
# ==============================================================================

# Construct dynamic filenames based on the statistic used
plot_filename <- paste0(output_dir, "R1-3_", target_statistic, "_Distribution_Estimates.pdf")
excel_filename <- paste0(output_dir, "R1-3_", target_statistic, "_Shapiro_Normality_Test.xlsx")

# Save plot to PDF
ggsave(
  filename = plot_filename,
  plot = plot,
  width = 15.27,
  height = 8.15,
  units = "in"
)

# Export normality test results to Excel
write_xlsx(normality, excel_filename)

message(paste("Files saved successfully:\n", plot_filename, "\n", excel_filename))
```

#VOY: ACABO DE HACER EL DE ARRIBA. DEBO REVISAR RESULTADOS ÁNIMO voy


## 1.4. Dependence of subsampling strategies on sample size

To evaluate whether different subsampling strategies exhibit distinct levels of sensitivity to sample size, we compared the strength of the correlation (Pcoef2) between sample size (`Sample_percent`) and subsampling performance across strategies.

We first filtered DF2 to include only frequentist estimates where the Pearson correlation was computed between sample size and design performance (`parameter1 == Sample_percent` and `parameter2 == Pearson_Correlation`) using the mean statistic.

Sampling strategies were then grouped into two categories: `"random"` (e.g., `ssample*`) and `"tessellation"` (e.g., grid-based designs such as `SLIC` or `basins`). For each group, we summarized the distribution of Pcoef2 values through descriptive statistics (mean, median, SD, IQR, CV).

Given the lack of normality (as previously assessed), we employed a non-parametric Wilcoxon rank-sum test to compare Pcoef2 distributions between the two strategy types.

Due to the lack of normality, a **Wilcoxon rank-sum test** was used to compare the distributions of Pcoef2 between groups. We also applied a
Fligner–Killeen test to evaluate differences in dispersion between strategies.

All results were visualized using violin and boxplots. Summary tables and test outputs were exported for reproducibility.

```{r}
# ==============================================================================
# SETUP AND LIBRARIES
# ==============================================================================

# Load required packages
library(readxl)
library(dplyr)
library(psych)
library(writexl)
library(ggplot2)
library(scales) # For squish logic in plots
library(broom)  # For tidying statistical test results

# Define input path (Adjust as needed)
input_path <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF2_SampleSize_vs_Subsample_performance_Correlations(Pcoeff2)_COMPLETE.xlsx"

# Define output path/prefix for saved files
output_dir <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/"

# Variable to control which statistic to analyze ("median" or "mean")
target_stat <- "median"

# ==============================================================================
# DATA LOADING AND PROCESSING
# ==============================================================================

# Load dataset
data <- read_excel(input_path)

# Filter valid Pearson correlations
data_frec <- data %>%
  filter(
    parameter1 == "Sample_percent",
    parameter2 == "Subsample_performance",
    effectsize == "Pearson correlation",
    Statistical_paradigm == "Frequentist", 
    Statistic == target_stat,
    !is.na(estimate)
  )

# 1. CLASSIFY AUDIO TYPE AND SAMPLING STRATEGY
data_frec <- data_frec %>%
  mutate(
    # Create Audio_Type based on Frequency_Range
    Audio_Type = ifelse(grepl("24-96", Frequency_Range), "Ultrasound", "Audible"),
    # Create Type_group (Random vs Tessellation)
    Type_group = case_when(
      grepl("^ssample", Subsample) ~ "Random",
      TRUE ~ "Tessellation"
    )
  )

# ==============================================================================
# STATISTICAL ANALYSIS (SPLIT BY AUDIO TYPE)
# ==============================================================================

# 1. Descriptive statistics (Grouped by Audio_Type AND Type_group)
summary_dependency <- data_frec %>%
  group_by(Audio_Type, Type_group) %>%
  summarise(
    n = n(),
    mean_Pcoef2 = mean(estimate, na.rm = TRUE),
    median_Pcoef2 = median(estimate, na.rm = TRUE),
    sd_Pcoef2 = sd(estimate, na.rm = TRUE),
    iqr_Pcoef2 = IQR(estimate, na.rm = TRUE),
    cv_Pcoef2 = sd(estimate, na.rm = TRUE) / abs(mean(estimate, na.rm = TRUE)),
    .groups = "drop"
  )

print("### Summary of sample size dependency (Split by Audio Type):")
print(summary_dependency)

# Function to run tests for a specific subset of data
run_tests <- function(df, label) {
  # Wilcoxon (Location)
  w_test <- wilcox.test(estimate ~ Type_group, data = df)
  w_tidy <- tidy(w_test) %>% mutate(Audio_Type = label, Test_Type = "Wilcoxon")
  
  # Fligner-Killeen (Homogeneity of Variance)
  fl_test <- fligner.test(estimate ~ Type_group, data = df)
  fl_tidy <- tidy(fl_test) %>% mutate(Audio_Type = label, Test_Type = "Fligner-Killeen")
  
  # F-test (Ratio of Variances)
  f_test <- var.test(estimate ~ Type_group, data = df)
  f_tidy <- data.frame(
    Audio_Type = label,
    Test_Type = "F-test (Variance)",
    F_statistic = f_test$statistic,
    df1 = f_test$parameter[1],
    df2 = f_test$parameter[2],
    p_value = f_test$p.value,
    CI_lower = f_test$conf.int[1],
    CI_upper = f_test$conf.int[2],
    Ratio_of_variances = f_test$estimate
  )
  
  return(list(wilcox = w_tidy, fligner = fl_tidy, ftest = f_tidy))
}

# Run tests for Audible
audible_data <- data_frec %>% filter(Audio_Type == "Audible")
res_audible <- run_tests(audible_data, "Audible")

# Run tests for Ultrasound
ultrasound_data <- data_frec %>% filter(Audio_Type == "Ultrasound")
res_ultrasound <- run_tests(ultrasound_data, "Ultrasound")

# Combine results into unified dataframes
wilcox_results <- bind_rows(res_audible$wilcox, res_ultrasound$wilcox)
fligner_results <- bind_rows(res_audible$fligner, res_ultrasound$fligner)
ftest_results <- bind_rows(res_audible$ftest, res_ultrasound$ftest)

print("### F-test Results (Split):")
print(ftest_results)

# ==============================================================================
# PLOTTING (FACETED BY AUDIO TYPE)
# ==============================================================================

final_plot <- ggplot(data_frec, aes(x = Type_group, y = estimate, fill = Type_group)) +
  # Violin
  geom_violin(alpha = 0.5, trim = FALSE, color = "darkblue", lwd = 0.3) +
  # Boxplot
  geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0.8, color = "black", lwd = 0.3) +
  # Colors
  scale_fill_manual(values = c("Random" = "#F8766D", "Tessellation" = "#00BFC4")) +
  # Facet Wrap to split the plot visually
  facet_wrap(~ Audio_Type) + 
  # Labels
  labs(title = paste("Sample size dependency (Pcoef2) -", target_stat),
       subtitle = "Split by Audio Type",
       y = "Pcoef2 (correlation)",
       x = "Sampling strategy") +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "none",
    strip.text = element_text(face = "bold", size = 12) # Style for "Audible"/"Ultrasound" labels
  ) +
  scale_y_continuous(limits = c(-1, 1), oob = scales::squish)

# Display plot
print(final_plot)

# ==============================================================================
# EXPORTING
# ==============================================================================

# Dynamic filenames
xlsx_name <- paste0(output_dir, "R1-4_", toupper(target_stat), "_STRATEGY_SPLIT_AUDIO.xlsx")
pdf_name  <- paste0(output_dir, "R1-4_", toupper(target_stat), "_Pcoef2_SPLIT_AUDIO_PLOT.pdf")

# Export to Excel (All results in structured sheets)
write_xlsx(
  list(
    Data_Filtered = data_frec,
    Summary_By_Audio = summary_dependency,
    F_Test_Results = ftest_results,
    Wilcoxon_Results = wilcox_results,
    Fligner_Results = fligner_results
  ),
  path = xlsx_name
)

# Export Plot
ggsave(
  filename = pdf_name,
  plot = final_plot,
  width = 10,  # Slightly wider to accommodate side-by-side plots
  height = 6,
  units = "in"
)

message(paste("Process complete. Files saved:\n", xlsx_name, "\n", pdf_name))

```


# 1.5 Sample size and subsample performance
```{r}
# ==============================================================================
# SETUP AND LIBRARIES
# ==============================================================================
library(readxl)
library(dplyr)
library(ggplot2)
library(scales) # Needed for squish

# Define paths (Adjust to your actual paths)
input_file <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF2_SampleSize_vs_Subsample_performance_Correlations(Pcoeff2)_COMPLETE.xlsx"
output_file <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/Figure4_Split.pdf" # Replace '#' with your actual output directory

# Choose statistic (consistent with your text/table, e.g., "median" or "mean")
target_stat <- "median"

# ==============================================================================
# DATA LOADING AND PROCESSING
# ==============================================================================

# 1. Load data
df <- read_excel(input_file, sheet= 'DATA')

# 2. Filter data
df_pcoef2 <- df %>%
  filter(
    parameter1 == "Sample_percent",
    parameter2 == "Subsample_performance", 
    effectsize == "Pearson correlation",
    Statistical_paradigm == "Frequentist", 
    Statistic == target_stat,              
    !is.na(estimate)
  )

# 3. Classify Strategy AND Audio Type
df_pcoef2 <- df_pcoef2 %>%
  mutate(
    # Strategy: Random vs Tessellation
    Strategy = ifelse(grepl("^ssample", Subsample), "Random", "Tessellation"),
    # Audio Type: Audible vs Ultrasound (based on "24-96")
    Audio_Type = ifelse(grepl("24-96", Frequency_Range), "Ultrasound", "Audible")
  )

# ==============================================================================
# PLOTTING
# ==============================================================================

# 4. Create Plot (Paper Style with Facets)
p <- ggplot(df_pcoef2, aes(x = Strategy, y = estimate, fill = Strategy)) +
  # Violin: Thinner border
  geom_violin(trim = FALSE, alpha = 0.6, color = "black", linewidth = 0.3) +
  # Boxplot: Nested, no outliers
  geom_boxplot(width = 0.18, outlier.shape = NA, alpha = 0.85, linewidth = 0.35) +
  # Split by Audio Type
  facet_wrap(~ Audio_Type) +
  # Manual Colors
  scale_fill_manual(
    values = c("Random" = "#F8766D", "Tessellation" = "#00BFC4")
  ) +
  # Limits -1 to 1 
  scale_y_continuous(limits = c(-1, 1), oob = scales::squish) +
  # Labels
  labs(
    x = NULL,
    y = "Pearson's correlation"
  ) +
  # Theme: Classic with Times New Roman
  theme_classic(base_family = "Times New Roman", base_size = 12) +
  theme(
    legend.position = "none",
    axis.text.x = element_text(size = 12, face = "bold", color = "black"),
    axis.text.y = element_text(size = 11, color = "black"),
    axis.title.y = element_text(size = 12, color = "black"),
    axis.line = element_line(linewidth = 0.4),
    # Style for the facet labels (Audible/Ultrasound)
    strip.background = element_blank(), # Clean background
    strip.text = element_text(size = 12, face = "bold", family = "Times New Roman", color = "black")
  )

# Display plot in RStudio
print(p)

# ==============================================================================
# EXPORT
# ==============================================================================

# 5. Save with EXACT dimensions
ggsave(
  filename = output_file,
  plot = p,
  width = 164.747,
  height = 107.244,
  units = "mm",
  device = cairo_pdf 
)

message(paste("Figure 4 saved to:", output_file))
```



# OBJECTIVE 2: EXPLORE EFFECTIVENESS OF STRATIFIED SAMPLING DESIGNS

Compare the performance of different sampling strategies, including random and stratified approaches, genarated by segmentation methods.

To assess the effectiveness of tessellation-based sampling strategies, we conducted a multi-step statistical analysis using R. First, we calculated two metrics for each subsample: (1) subsample performance in the form of Pearson correlation between the subsample and its full reference dataset; and (2) sampling efficiency, computed as the ratio between Pearson correlation and the sampling effort (Subsample performance/Sample_percent).

We then grouped data by summary statistic type (mean or median) and computed descriptive statistics for each tessellation strategy, including sample size, central tendency (mean, median), dispersion (standard deviation), and coefficient of variation.

Normality of Pearson correlations was assessed using the Shapiro–Wilk test within each subsampling group, using shapiro_test() from the rstatix package.

Given widespread non-normality, we applied Kruskal–Wallis rank-sum tests to compare both raw performance and efficiency across subsampling strategies. Where differences were detected, we used Dunn’s post-hoc tests with Benjamini–Hochberg correction to evaluate pairwise contrasts.

Violin and boxplots were generated to visually compare the distribution of both performance and efficiency across subsampling methods. All procedures were run separately for mean and median statistics, and results (summary tables and plots) were exported to PDF and Excel formats.

Packages we used: dplyr, ggplot2, ggpubr, rstatix, FSA, writexl.

The entry database is: DF1-1_Interpolations_variograms_and error_for_subsamples.xlsx

```{r}
# ==============================================================================
# SETUP AND LIBRARIES
# ==============================================================================
library(readxl)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(rstatix)
library(writexl)
library(scales) 
library(patchwork) # INDISPENSABLE para unir gráficas (arriba/abajo)

# 1. Define Paths
input_path <- "/home/vmartinezarias/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF1-1_Interpolations_variograms_error_and_performance_of_subsamples_US.xlsx"
output_dir <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/"

# ==============================================================================
# DATA LOADING & PREP
# ==============================================================================

df <- read_excel(input_path) %>%
  rename(Pearson_Correlation = Subsample_performance) %>% 
  mutate(across(c(Subsample, Type, Index, Statistic), as.character)) %>%
  mutate(
    Pearson_rel = Pearson_Correlation / Sample_percent,
    Audio_Type = ifelse(grepl("24-96", Frequency_Range), "Ultrasound", "Audible")
  )

# ==============================================================================
# LOOP: MEAN & MEDIAN
# ==============================================================================

for (stat in c("mean", "median")) {
  
  # Filter: Only Tessellations, current statistic
  df_tess <- df %>%
    filter(Type == "Tessells", Statistic == stat, !is.na(Pearson_Correlation))
  
  # --------------------------------------------------------------------------
  # 1. STATISTICS (Keep your existing calculation logic)
  # --------------------------------------------------------------------------
  
  # Summary Stats
  summary_precision <- df_tess %>%
    group_by(Audio_Type, Subsample) %>%
    summarise(
      n = n(),
      mean_Pearson = mean(Pearson_Correlation, na.rm = TRUE),
      median_Pearson = median(Pearson_Correlation, na.rm = TRUE),
      sd_Pearson = sd(Pearson_Correlation, na.rm = TRUE),
      mean_eff = mean(Pearson_rel, na.rm = TRUE),
      median_eff = median(Pearson_rel, na.rm = TRUE),
      sd_eff = sd(Pearson_rel, na.rm = TRUE),
      cv_eff = sd_eff / abs(mean_eff),
      .groups = "drop"
    ) %>%
    arrange(Audio_Type, desc(median_eff))
  
  # Statistical Tests Function
  run_stat_tests <- function(sub_data, label) {
    kw_p <- kruskal_test(Pearson_Correlation ~ Subsample, data = sub_data) %>% mutate(Metric="Precision", Audio=label)
    ph_p <- dunn_test(Pearson_Correlation ~ Subsample, data = sub_data, p.adjust.method = "BH") %>% mutate(Metric="Precision", Audio=label)
    kw_e <- kruskal_test(Pearson_rel ~ Subsample, data = sub_data) %>% mutate(Metric="Efficiency", Audio=label)
    ph_e <- dunn_test(Pearson_rel ~ Subsample, data = sub_data, p.adjust.method = "BH") %>% mutate(Metric="Efficiency", Audio=label)
    return(list(kw = bind_rows(kw_p, kw_e), posthoc = bind_rows(ph_p, ph_e)))
  }
  
  res_aud <- run_stat_tests(df_tess %>% filter(Audio_Type == "Audible"), "Audible")
  res_us <- run_stat_tests(df_tess %>% filter(Audio_Type == "Ultrasound"), "Ultrasound")
  
  # Export Stats to Excel
  write_xlsx(
    list(
      Summary_Stats = summary_precision,
      KW_Tests = bind_rows(res_aud$kw, res_us$kw),
      PostHoc = bind_rows(res_aud$posthoc, res_us$posthoc)
    ),
    path = paste0(output_dir, "KW_Tessellation_Stats_", stat, "_Split.xlsx")
  )
  
  # --------------------------------------------------------------------------
  # 2. PLOTTING (The Combined Figure 5)
  # --------------------------------------------------------------------------
  
  # Define a consistent theme for the paper
  paper_theme <- theme_bw(base_family = "Times New Roman", base_size = 11) +
    theme(
      legend.position = "none",
      axis.text = element_text(color = "black"),
      axis.title = element_text(face = "bold"),
      strip.background = element_blank(),
      strip.text = element_text(face = "bold", size = 12),
      panel.grid.major.x = element_blank(),
      panel.grid.minor = element_blank()
    )
  
  # --- PLOT A: EFFICIENCY (TOP) ---
  p_eff <- ggplot(df_tess, aes(x = Subsample, y = Pearson_rel, fill = Subsample)) +
    geom_violin(trim = FALSE, alpha = 0.8, size = 0.2) +
    geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5, size = 0.2) +
    facet_wrap(~ Audio_Type) + # Split columns
    labs(y = "Sampling efficiency\n(Pearson / % sample)", x = NULL) +
    paper_theme +
    theme(
      axis.text.x = element_blank(),  # Hide X labels (redundant with bottom plot)
      axis.ticks.x = element_blank(),
      axis.title.x = element_blank()
    )
  
  # --- PLOT B: PERFORMANCE/PRECISION (BOTTOM) ---
  p_acc <- ggplot(df_tess, aes(x = Subsample, y = Pearson_Correlation, fill = Subsample)) +
    geom_violin(trim = FALSE, alpha = 0.8, size = 0.2) +
    geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5, size = 0.2) +
    facet_wrap(~ Audio_Type) + # Split columns
    scale_y_continuous(limits = c(-1, 1), oob = scales::squish) +
    labs(y = "Subsample performance\n(Pearson's correlation)", x = NULL) +
    paper_theme +
    theme(
      strip.text = element_blank(), # Hide facet titles (already on top plot)
      axis.text.x = element_text(angle = 0, hjust = 0.5, size = 10, face = "bold") 
    )
  
  # --- COMBINE PLOTS WITH PATCHWORK ---
  final_plot <- p_eff / p_acc + 
    plot_layout(heights = c(1, 1)) # Equal height
  
  # --------------------------------------------------------------------------
  # 3. EXPORT FIGURE (Exact dimensions)
  # --------------------------------------------------------------------------
  ggsave(
    filename = paste0(output_dir, "Figure5_Tessellation_Composite_", stat, ".pdf"),
    plot = final_plot,
    width = 170,   # Max width for paper
    height = 160,  # Adjusted height for two stacked panels
    units = "mm",
    device = cairo_pdf
  )
  
  message(paste("Generated stats and Figure 5 for:", stat))
}

message("Process complete. Check output folder.")
```




# OBJECTIVE 3:

Specific objective 3: To explore how the interplay of categorical (e.g., landscape, time, index) and continuous variables (e.g., sample size, performance) explains the multivariate structure that differentiates stratified subsampling strategies.

## 3.1 FAMD analyses

A Factorial Analysis of Mixed Data (FAMD) is conducted to explore the multivariate structure of stratified subsampling strategies by integrating both categorical and continuous variables. The dataset was filtered to retain only relevant observations (e.g., median statistics from tessellation-based designs) and included key landscape, temporal, acoustic, and performance variables. The variable Subsample, which defines the sampling strategy under evaluation, was treated as a supplementary qualitative variable to avoid influencing the dimensional construction.

FAMD was performed using the FactoMineR package, extracting the coordinates of individuals, the contributions of quantitative and qualitative variables, and the proportion of variance explained by each dimension. Results were organized into separate Excel sheets for subsequent visualization and statistical analysis.

```{r}
# ==============================================================================
# SCRIPT: DUAL FAMD ANALYSIS (AUDIBLE & ULTRASOUND)
# ==============================================================================

# 1. Load required libraries
library(readxl)
library(FactoMineR)
library(factoextra)
library(dplyr)
library(openxlsx)

# 2. Load dataset
# Replace the path with your actual file path
file_path <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF1-1_Interpolations_variograms_error_and_performance_of_subsamples_US.xlsx"
data <- read_excel(file_path, sheet = "DATA")

# 3. Pre-process: Define Audio Categories
# We create a split column based on your Frequency Ranges
data <- data %>%
  mutate(Audio_Category = case_when(
    Frequency_Range %in% c("R24-96") ~ "Ultrasound",
    Frequency_Range %in% c("R0-2", "R2-12", "R0-24", "R12-24") ~ "Audible",
    TRUE ~ "Unknown"
  ))

# 4. Define a function to perform FAMD
# This allows us to run the exact same logic twice (for Audible and Ultra)
perform_famd_and_extract <- function(input_data, prefix) {
  
  message(paste("Processing:", prefix, "| Rows:", nrow(input_data)))
  
  # --- A. Filter and Select ---
  # Using the logic from your snippet:
  # 1. Filter valid values (Subsample_performance represents Pearson_Correlation here)
  # 2. Filter Statistic and Type (Adjust "median" to "mean" if your data uses mean)
  data_famd <- input_data %>%
    filter(
      !is.na(Subsample_performance), 
      Statistic == "median",  #<-- WARNING: Check if your Excel says "median" or "mean"
      Type == "Tessells"      #<-- WARNING: Check if "Tessells" is the exact string
    ) %>%
    select(Locality, Frequency_Range, Time_Period, Index, Subsample, Sample_percent, Subsample_performance)
  
  # Check if data exists after filtering
  if(nrow(data_famd) < 5) {
    warning(paste("Not enough data for", prefix, "- Check your filters (Statistic/Type)."))
    return(NULL)
  }

  # --- B. Convert to Factors ---
  data_famd <- data_famd %>%
    mutate(across(c(Locality, Frequency_Range, Time_Period, Index, Subsample), as.factor)) %>%
    mutate(
      Sample_percent = as.numeric(Sample_percent),
      Subsample_performance = as.numeric(Subsample_performance)
    )

  # --- C. Clean Constant Variables ---
  # FAMD crashes if a factor has only 1 level. We remove them.
  data_famd <- data_famd %>% select(where(~ is.numeric(.) || n_distinct(.) > 1))

  # --- D. Identify Supplementary Variable ---
  if ("Subsample" %in% names(data_famd)) {
    supplementary_var_index <- which(colnames(data_famd) == "Subsample")
  } else {
    supplementary_var_index <- NULL
  }

  # --- E. Run FAMD ---
  famd_result <- FAMD(data_famd, 
                      ncp = 5, 
                      sup.var = supplementary_var_index, 
                      graph = FALSE)

  # --- F. Extract Results ---
  # 1. Individuals
  coord_ind <- as.data.frame(famd_result$ind$coord)
  coord_ind$ID <- rownames(coord_ind)

  # 2. Quantitative Variables
  coord_quanti <- as.data.frame(famd_result$quanti.var$coord)
  coord_quanti$Variable <- rownames(coord_quanti)

  # 3. Qualitative Variables
  coord_quali <- as.data.frame(famd_result$quali.var$coord)
  coord_quali$Category <- rownames(coord_quali)

  # 4. Eigenvalues (Variance)
  eig <- as.data.frame(famd_result$eig)
  colnames(eig) <- c("Eigenvalue", "Variance_Percent", "Cumulative_Variance_Percent")
  eig$Dimension <- rownames(eig)

  # Return a list containing all dataframes
  return(list(
    ind = coord_ind, 
    quanti = coord_quanti, 
    quali = coord_quali, 
    eig = eig
  ))
}

# 5. Execute Analysis for both groups
# Split data
data_audible <- data %>% filter(Audio_Category == "Audible")
data_ultra   <- data %>% filter(Audio_Category == "Ultrasound")

# Run FAMD
results_audible <- perform_famd_and_extract(data_audible, "Audible")
results_ultra   <- perform_famd_and_extract(data_ultra, "Ultrasound")

# 6. Save to Excel
wb <- createWorkbook()

# Helper function to add sheets
add_results_to_wb <- function(wb, res_list, prefix) {
  if (is.null(res_list)) return()
  
  addWorksheet(wb, paste0(prefix, "_Ind_Coord"))
  writeData(wb, paste0(prefix, "_Ind_Coord"), res_list$ind)
  
  addWorksheet(wb, paste0(prefix, "_Quanti_Var"))
  writeData(wb, paste0(prefix, "_Quanti_Var"), res_list$quanti)
  
  addWorksheet(wb, paste0(prefix, "_Quali_Var"))
  writeData(wb, paste0(prefix, "_Quali_Var"), res_list$quali)
  
  addWorksheet(wb, paste0(prefix, "_Explained_Var"))
  writeData(wb, paste0(prefix, "_Explained_Var"), res_list$eig)
}

# Add Audible sheets
if (!is.null(results_audible)) {
  add_results_to_wb(wb, results_audible, "AUD")
}

# Add Ultrasound sheets
if (!is.null(results_ultra)) {
  add_results_to_wb(wb, results_ultra, "ULTRA")
}

# Export
output_file <- paste0(dirname(file_path), "/FAMD_results_Audible_and_Ultrasound.xlsx")
saveWorkbook(wb, output_file, overwrite = TRUE)

message("SUCCESS: FAMD analysis completed for both groups. Results saved.")
```

## 3.2 Graphical abstract plots for FAMD

This script generates a visual summary of the Factorial Analysis of Mixed Data (FAMD) results. It includes a scree plot showing the variance explained by each dimension, along with bar plots displaying the top contributing variables to Dimensions 1 through 4. All plots are styled consistently using a minimal theme and arranged using the patchwork package into a single combined figure. The final output is saved as a high-resolution PDF for publication or reporting purposes.

```{r}
# ==============================================================================
# 3.2 GRAPHICAL ABSTRACT PLOTS (AUDIBLE & ULTRASOUND)
# ==============================================================================

# Load required libraries
library(FactoMineR)
library(factoextra)
library(patchwork) # For combining ggplot-based charts
library(ggplot2)
library(dplyr)

# ------------------------------------------------------------------------------
# FUNCTION: GENERATE AND SAVE FAMD PLOTS
# ------------------------------------------------------------------------------
# This function re-runs the FAMD calculation (needed for plotting) 
# and generates the PDF with your specific layout.

create_famd_abstract_plot <- function(input_data, category_name) {
  
  message(paste("Generating Graphical Abstract for:", category_name))
  
  # --- 1. DATA PREPARATION (Must match analysis logic) ---
  # Filter valid values and select columns
  df_plot <- input_data %>%
    filter(!is.na(Subsample_performance)) %>%
    select(Locality, Frequency_Range, Time_Period, Index, Subsample, Sample_percent, Subsample_performance) %>%
    na.omit()
  
  # Remove constant variables (FAMD crashes if variance is 0)
  df_plot <- df_plot %>% select(where(~ is.numeric(.) || n_distinct(.) > 1))
  
  # Safety check
  if(nrow(df_plot) < 5) {
    warning(paste("Skipping", category_name, "- Not enough data."))
    return(NULL)
  }

  # Identify supplementary variable
  sup_var_idx <- if ("Subsample" %in% names(df_plot)) which(colnames(df_plot) == "Subsample") else NULL

  # --- 2. RUN FAMD (Required for fviz functions) ---
  famd_result <- FAMD(df_plot, ncp = 5, sup.var = sup_var_idx, graph = FALSE)

  # --- 3. CREATE INDIVIDUAL PLOTS ---
  
  # Scree plot
  plot_scree <- fviz_screeplot(famd_result, addlabels = TRUE, 
                               ggtheme = theme_minimal()) +
    ggtitle(paste("Scree Plot -", category_name)) +
    coord_cartesian(ylim = c(0, 50)) # Adjusted Y-limit to see variance clearly

  # Contribution plots (Dims 1 to 4)
  plot_contrib_dim1 <- fviz_contrib(famd_result, choice = "var", axes = 1, top = 10,
                                    ggtheme = theme_minimal()) +
    ggtitle("Contribution to Dim.1 (Top 10)") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  plot_contrib_dim2 <- fviz_contrib(famd_result, choice = "var", axes = 2, top = 10,
                                    ggtheme = theme_minimal()) +
    ggtitle("Contribution to Dim.2 (Top 10)") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  plot_contrib_dim3 <- fviz_contrib(famd_result, choice = "var", axes = 3, top = 10,
                                    ggtheme = theme_minimal()) +
    ggtitle("Contribution to Dim.3 (Top 10)") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  plot_contrib_dim4 <- fviz_contrib(famd_result, choice = "var", axes = 4, top = 10,
                                    ggtheme = theme_minimal()) +
    ggtitle("Contribution to Dim.4 (Top 10)") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  # --- 4. COMBINE PLOTS (PATCHWORK) ---
  combined_plot <- plot_scree /
    (plot_contrib_dim1 + plot_contrib_dim2) /
    (plot_contrib_dim3 + plot_contrib_dim4)

  # Add global title
  combined_plot <- combined_plot + 
    plot_annotation(
      title = paste("      FAMD -", category_name),
      theme = theme(plot.title = element_text(size = 16, face = "bold"))
    )

  # --- 5. SAVE PDF ---
  # Using the directory of the original file path defined in previous chunk
  # If 'file_path' is lost, hardcode your path here.
  out_path <- paste0(dirname(file_path), "/FAMD_summary_plot_", category_name, ".pdf")
  
  ggsave(
    filename = out_path,
    plot = combined_plot,
    device = "pdf",
    width = 8,
    height = 10,
    units = "in"
  )
  
  message(paste("PDF saved:", out_path))
}

# ------------------------------------------------------------------------------
# EXECUTE FOR BOTH GROUPS
# ------------------------------------------------------------------------------
# Assumes 'data_audible' and 'data_ultra' exist from the previous script.
# If not, reload your data before running this.

if(exists("data_audible")) {
  create_famd_abstract_plot(data_audible, "AUDIBLE")
} else {
  warning("data_audible object not found.")
}

if(exists("data_ultra")) {
  create_famd_abstract_plot(data_ultra, "ULTRASOUND")
} else {
  warning("data_ultra object not found.")
}
```



## 3.3 Variables projection and individuals plot

```{r}
# ==============================================================================
# 3.3 VARIABLES PROJECTION AND INDIVIDUALS PLOT (AUDIBLE & ULTRASOUND)
# ==============================================================================

# 1. Load required libraries
library(readxl)
library(FactoMineR)
library(factoextra)
library(dplyr)
library(ggplot2)
library(patchwork)

# ------------------------------------------------------------------------------
# FUNCTION: GENERATE PROJECTION PLOTS
# ------------------------------------------------------------------------------
create_projection_plots <- function(input_data, category_name, file_path_ref) {
  
  message(paste("Generating Projection Plots for:", category_name))
  
  # --- A. Data Preparation ---
  # Filter based on your specific criteria (Median, Tessells, etc.)
  # Note: I replaced 'Pearson_Correlation' with 'Subsample_performance' based on previous context.
  df_famd <- input_data %>%
    filter(
      !is.na(Subsample_performance),
      # Uncomment and adjust these lines if you need strict filtering:
      # Statistic == "median", 
      # Type == "Tessells"
    ) %>%
    select(Locality, Frequency_Range, Time_Period, Index, Subsample, Sample_percent, Subsample_performance) %>%
    na.omit()
  
  # Remove constant variables to prevent FAMD errors
  df_famd <- df_famd %>% select(where(~ is.numeric(.) || n_distinct(.) > 1))
  
  # Ensure factors are correct
  df_famd <- df_famd %>% 
    mutate(across(where(is.character), as.factor)) %>%
    mutate(across(where(is.factor), droplevels))
  
  # Check data sufficiency
  if(nrow(df_famd) < 5) {
    warning(paste("Skipping", category_name, "- Not enough data points."))
    return(NULL)
  }

  # --- B. Run FAMD ---
  # Identify supplementary variable index if it exists
  sup_idx <- if ("Subsample" %in% names(df_famd)) which(colnames(df_famd) == "Subsample") else NULL
  
  famd_result <- FAMD(df_famd, ncp = 5, sup.var = sup_idx, graph = FALSE)
  
  # --- C. Variable Projection Plots ---
  # Dim 1 vs 2
  plot_var_12 <- fviz_famd_var(famd_result, axes = c(1, 2), repel = TRUE) + 
    ggtitle("Variables: Dim 1 vs Dim 2") +
    theme_minimal()
  
  # Dim 2 vs 3
  plot_var_23 <- fviz_famd_var(famd_result, axes = c(2, 3), repel = TRUE) + 
    ggtitle("Variables: Dim 2 vs Dim 3") +
    theme_minimal()
  
  # Dim 1 vs 3
  plot_var_13 <- fviz_famd_var(famd_result, axes = c(1, 3), repel = TRUE) + 
    ggtitle("Variables: Dim 1 vs Dim 3") +
    theme_minimal()

  # --- D. Individual Projection Plots ---
  # Note: habillage colors points by 'Subsample'. addEllipses draws confidence ovals.
  
  # Dim 1 vs 2
  plot_ind_12 <- fviz_famd_ind(famd_result, axes = c(1, 2), habillage = "Subsample", 
                               addEllipses = TRUE, label = "none") + 
    ggtitle("Individuals: Dim 1 vs Dim 2") +
    theme_minimal()
  
  # Dim 2 vs 3
  plot_ind_23 <- fviz_famd_ind(famd_result, axes = c(2, 3), habillage = "Subsample", 
                               addEllipses = TRUE, label = "none") + 
    ggtitle("Individuals: Dim 2 vs Dim 3") +
    theme_minimal()
  
  # Dim 1 vs 3
  plot_ind_13 <- fviz_famd_ind(famd_result, axes = c(1, 3), habillage = "Subsample", 
                               addEllipses = TRUE, label = "none") + 
    ggtitle("Individuals: Dim 1 vs Dim 3") +
    theme_minimal()

  # --- E. Combine Plots (Patchwork) ---
  # Layout: Left column = Variables, Right column = Individuals
  combined_plot <- (plot_var_12 | plot_ind_12) /
                   (plot_var_23 | plot_ind_23) /
                   (plot_var_13 | plot_ind_13) +
                   plot_annotation(
                     title = paste("FAMD Projections:", category_name),
                     subtitle = "Left: Variables | Right: Individuals (colored by Subsample)",
                     theme = theme(plot.title = element_text(size = 16, face = "bold"))
                   )

  # --- F. Save to PDF ---
  # Construct output path based on the input file location
  out_path <- paste0(dirname(file_path_ref), "/FAMD_variables_individuals_clean_", category_name, ".pdf")
  
  ggsave(
    filename = out_path,
    plot = combined_plot,
    device = "pdf",
    width = 12,
    height = 14,
    units = "in"
  )
  
  message(paste("Saved plot to:", out_path))
}

# ------------------------------------------------------------------------------
# 2. EXECUTION
# ------------------------------------------------------------------------------
# Define the file path (used for saving outputs)
# Update this path if necessary
input_path_ref <- "~/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF1-1_Interpolations_variograms_error_and_performance_of_subsamples_US.xlsx"

# NOTE: Assumes 'data_audible' and 'data_ultra' exist from the previous chunks.
# If they don't, please run the "Split Data" chunk first.

if(exists("data_audible")) {
  create_projection_plots(data_audible, "AUDIBLE", input_path_ref)
} else {
  warning("data_audible object not found. Run previous split logic.")
}

if(exists("data_ultra")) {
  create_projection_plots(data_ultra, "ULTRASOUND", input_path_ref)
} else {
  warning("data_ultra object not found. Run previous split logic.")
}
```

## 3.4 KRUSKALL TEST

This graph presents the results of a statistical analysis (Kruskal-Wallis + post-hoc Dunn) applied to the coordinates of dimension 3 (Dim.3) of the FAMD analysis, with the aim of comparing the stratified sampling methods (Subsample: basins, halfbasins, SLIC).

```{r}
# ==============================================================================
# SCRIPT: FAMD & KRUSKAL-WALLIS ANALYSIS (AUDIBLE VS ULTRASOUND)
# ==============================================================================
# Objective: Compare subsampling strategies using FAMD coordinates.
# Logic: 'Subsample' is treated as an ACTIVE variable in FAMD.
# Outputs: 
#   1. Excel: Detailed Statistics (Chi2, p, df, n, letters) AND Pairwise comparisons.
#   2. PDF 1: Heatmap Comparison (Summary View).
#   3. PDF 2: Detailed Boxplots (BIGGER FONT, 240x140mm).
# ==============================================================================

# 1. Load Required Libraries
library(readxl)
library(FactoMineR)
library(dplyr)
library(rstatix)
library(multcompView)
library(ggplot2)
library(ggpubr)    
library(patchwork) 
library(tidyr)
library(stringr)
library(viridis)
library(openxlsx) 

# 2. Define Data Path
# ------------------------------------------------------------------------------
file_path <- "/home/vmartinezarias/Dropbox/Doctorado_Personal_Folder/ThesisAnalyses/PROVISIONALES_PAPERS/PAPER01 SAMPLING DESIGN/CORRECTIONS/DB_for_corrections/DF1-1_Interpolations_variograms_error_and_performance_of_subsamples_US.xlsx"
data_raw <- read_excel(file_path, sheet = "DATA")

# 3. Master Analysis Function
# ------------------------------------------------------------------------------
analyze_active_subsample <- function(data_input, group_name) {
  
  # --- A. Data Cleaning and Filtering ---
  df_work <- data_input %>%
    filter(
      !is.na(Subsample_performance), 
      Statistic == "median", 
      Type == "Tessells"
    ) %>%
    select(Locality, Frequency_Range, Time_Period, Index, Subsample, Sample_percent, Subsample_performance) %>%
    na.omit() %>%
    mutate(across(where(is.character), as.factor)) %>%
    mutate(across(where(is.factor), droplevels))
  
  n_samples <- nrow(df_work)
  cat(paste("\n>>> Processing Group:", group_name, "| Sample Size (n):", n_samples, "\n"))
  
  if(n_samples < 10) {
    warning(paste("Insufficient data for group:", group_name))
    return(NULL)
  }
  
  # --- B. FAMD Calculation ---
  famd_res <- FAMD(df_work, ncp = 5, graph = FALSE)
  coords <- as.data.frame(famd_res$ind$coord[, 1:4])
  coords$Subsample <- df_work$Subsample
  
  # --- C. Statistical Loop ---
  plot_data <- data.frame() 
  stats_data <- data.frame() 
  pairwise_data <- data.frame() 
  
  for (dim in c("Dim.1", "Dim.2", "Dim.3", "Dim.4")) {
    
    # 1. Kruskal-Wallis Test
    kw <- kruskal_test(as.formula(paste(dim, "~ Subsample")), data = coords)
    
    final_letters <- rep("a", length(levels(coords$Subsample)))
    names(final_letters) <- levels(coords$Subsample)
    
    # 2. Post-Hoc (Dunn)
    if (kw$p < 0.05) {
      dunn <- dunn_test(as.formula(paste(dim, "~ Subsample")), data = coords, p.adjust.method = "BH")
      
      # Save Pairwise
      dunn_clean <- dunn %>%
        mutate(Dimension = dim, Category = group_name) %>%
        select(Category, Dimension, group1, group2, n1, n2, statistic, p, p.adj, p.adj.signif)
      pairwise_data <- bind_rows(pairwise_data, dunn_clean)
      
      grps <- levels(coords$Subsample)
      pmat <- matrix(1, nrow=length(grps), ncol=length(grps), dimnames=list(grps, grps))
      diag(pmat) <- 0
      
      for(i in 1:nrow(dunn)){
        pmat[dunn$group1[i], dunn$group2[i]] <- dunn$p.adj[i]
        pmat[dunn$group2[i], dunn$group1[i]] <- dunn$p.adj[i]
      }
      
      try({
        lets <- multcompLetters(pmat < 0.05)$Letters
        final_letters <- lets[grps]
      }, silent = TRUE)
    }
    
    # Data for Heatmap
    tmp_plot <- data.frame(
      Subsample = names(final_letters),
      Significance_Letter = as.character(final_letters),
      Dimension = dim,
      Category = group_name
    )
    plot_data <- rbind(plot_data, tmp_plot)
    
    # Data for Excel
    tmp_stats <- data.frame(
      Dimension = dim,
      KW_Statistic = round(kw$statistic, 3),
      KW_df = kw$df,
      KW_n = kw$n,
      KW_p_value = format.pval(kw$p, digits = 3, eps = 0.001)
    )
    for(g in names(final_letters)) {
      tmp_stats[[paste0("Letter_", g)]] <- final_letters[g]
    }
    stats_data <- bind_rows(stats_data, tmp_stats)
  }
  
  return(list(
    coords = coords, 
    heatmap_data = plot_data, 
    stats = stats_data, 
    pairwise = pairwise_data 
  ))
}

# 4. Helper Function: Boxplots (FONT SIZE INCREASED + BOLD AXIS)
# ------------------------------------------------------------------------------
create_detailed_boxplots <- function(res_obj, group_name) {
  coords <- res_obj$coords
  stats_df <- res_obj$stats
  plot_list <- list()
  font_family <- "serif" # Times New Roman
  
  for(dim in c("Dim.1", "Dim.2", "Dim.3", "Dim.4")) {
    row_stats <- stats_df %>% filter(Dimension == dim)
    letters_named <- row_stats %>% select(starts_with("Letter_"))
    
    anno_df <- data.frame(
      Subsample = names(letters_named) %>% str_remove("Letter_"),
      Letter = as.character(letters_named[1,])
    )
    
    y_max_vals <- coords %>% 
      group_by(Subsample) %>% 
      summarise(max_val = max(.data[[dim]], na.rm = TRUE))
    anno_df <- left_join(anno_df, y_max_vals, by = "Subsample")
    
    # Format p-value string properly
    p_val_str <- row_stats$KW_p_value
    if(grepl("<", p_val_str)) {
      p_label <- paste("p", p_val_str)
    } else {
      p_label <- paste("p =", p_val_str)
    }
    
    subtitle_text <- paste0("KW Chi2 = ", row_stats$KW_Statistic, ", ", p_label, ", n = ", row_stats$KW_n)

    p <- ggplot(coords, aes(x = Subsample, y = .data[[dim]], fill = Subsample)) +
      stat_boxplot(geom = "errorbar", width = 0.2, lwd = 0.6) +
      geom_boxplot(outlier.shape = 19, outlier.size = 2, alpha = 0.8, lwd = 0.6) +
      
      # Letters above bars (Bold, larger size)
      geom_text(data = anno_df, aes(x = Subsample, y = max_val, label = Letter), 
                vjust = -0.8, size = 6, family = font_family, fontface = "bold", inherit.aes = FALSE) +
      
      scale_fill_brewer(palette = "Set2") + 
      
      labs(y = paste("FAMD Coordinate -", dim), x = NULL, 
           title = paste(dim),
           subtitle = subtitle_text) +
      
      theme_minimal(base_family = font_family, base_size = 14) + # Increased Base Size to 14
      theme(
        legend.position = "none",
        plot.title = element_text(face = "bold", size = 16),
        plot.subtitle = element_text(size = 12, face = "italic"),
        axis.title.y = element_text(face = "bold", size = 14), # Bold Axis Title
        axis.text.x = element_text(face = "bold", size = 12, color = "black"), # Bold Axis Text
        axis.text.y = element_text(face = "bold", size = 12, color = "black"),
        panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
        panel.grid.minor = element_blank()
      )
    
    p <- p + scale_y_continuous(expand = expansion(mult = c(0.1, 0.2)))
    plot_list[[dim]] <- p
  }
  
  final_grid <- (plot_list[["Dim.1"]] | plot_list[["Dim.2"]]) / 
                (plot_list[["Dim.3"]] | plot_list[["Dim.4"]]) +
                plot_annotation(
                  title = paste("Detailed Distribution -", group_name),
                  theme = theme(
                    text = element_text(family = font_family),
                    plot.title = element_text(face = "bold", size = 18, hjust = 0.5)
                  )
                )
  return(final_grid)
}

# 5. Execute Analysis
# ------------------------------------------------------------------------------
data_audible <- data_raw %>% filter(Frequency_Range %in% c("R0-2", "R2-12", "R0-24", "R12-24"))
data_ultra   <- data_raw %>% filter(Frequency_Range %in% c("R24-96"))

res_aud <- analyze_active_subsample(data_audible, "AUDIBLE")
res_us  <- analyze_active_subsample(data_ultra, "ULTRASOUND")

# 6. Output 1: Excel Table (Detailed Stats + Pairwise)
# ------------------------------------------------------------------------------
wb <- createWorkbook()

if(!is.null(res_aud)) {
  addWorksheet(wb, "Audible_Summary")
  writeData(wb, "Audible_Summary", res_aud$stats)
  addStyle(wb, "Audible_Summary", createStyle(textDecoration = "bold"), rows = 1, cols = 1:ncol(res_aud$stats))
  
  addWorksheet(wb, "Audible_Pairwise")
  writeData(wb, "Audible_Pairwise", res_aud$pairwise)
  addStyle(wb, "Audible_Pairwise", createStyle(textDecoration = "bold"), rows = 1, cols = 1:ncol(res_aud$pairwise))
}

if(!is.null(res_us)) {
  addWorksheet(wb, "Ultrasound_Summary")
  writeData(wb, "Ultrasound_Summary", res_us$stats)
  addStyle(wb, "Ultrasound_Summary", createStyle(textDecoration = "bold"), rows = 1, cols = 1:ncol(res_us$stats))
  
  addWorksheet(wb, "Ultrasound_Pairwise")
  writeData(wb, "Ultrasound_Pairwise", res_us$pairwise)
  addStyle(wb, "Ultrasound_Pairwise", createStyle(textDecoration = "bold"), rows = 1, cols = 1:ncol(res_us$pairwise))
}

excel_path <- paste0(dirname(file_path), "/Statistical_Results_Detailed_With_Pairwise.xlsx")
saveWorkbook(wb, excel_path, overwrite = TRUE)
message(paste("1. Excel saved with Pairwise comparisons to:", excel_path))

# 7. Output 2: Heatmap Comparison (PDF 170mm)
# ------------------------------------------------------------------------------
df_plot_combined <- bind_rows(res_aud$heatmap_data, res_us$heatmap_data)

if(nrow(df_plot_combined) > 0) {
  df_plot_combined$Dimension <- factor(df_plot_combined$Dimension, levels = c("Dim.1", "Dim.2", "Dim.3", "Dim.4"))
  # Note: Adjust labels here if you want full text
  df_plot_combined$Category <- factor(df_plot_combined$Category, levels = c("AUDIBLE", "ULTRASOUND"))
  
  font_family <- "serif"
  
  p_heat <- ggplot(df_plot_combined, aes(x = Dimension, y = Subsample, fill = Significance_Letter)) +
    geom_tile(color = "white", lwd = 0.5) +
    geom_text(aes(label = Significance_Letter), 
              family = font_family, fontface = "bold", size = 5, color = "black") +
    facet_wrap(~Category, ncol = 2) +
    scale_fill_viridis_d(option = "plasma", direction = -1, begin = 0.1, end = 0.9, name = "Sig. Group") +
    labs(x = "FAMD Dimension", y = "Subsampling Strategy") +
    theme_minimal(base_family = font_family, base_size = 12) +
    theme(
      strip.text = element_text(face = "bold", size = 11),
      axis.text = element_text(color = "black"),
      axis.title = element_text(face = "bold"),
      panel.grid = element_blank(),
      panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
      legend.position = "right"
    )
  
  heatmap_path <- paste0(dirname(file_path), "/Heatmap_Comparison_170mm.pdf")
  ggsave(heatmap_path, p_heat, device = cairo_pdf, width = 170, height = 90, units = "mm")
  message(paste("2. Heatmap PDF saved to:", heatmap_path))
}

# 8. Output 3: Detailed Boxplots (UPDATED SIZE & FONT)
# ------------------------------------------------------------------------------
# Width = 240mm, Height = 140mm
# Font: Times New Roman, Size 12+ (handled in function via base_size=14)

if(!is.null(res_aud)) {
  box_aud <- create_detailed_boxplots(res_aud, "AUDIBLE")
  box_aud_path <- paste0(dirname(file_path), "/Detailed_Boxplots_Audible_240mm.pdf")
  ggsave(box_aud_path, box_aud, device = cairo_pdf, width = 240, height = 140, units = "mm")
  message(paste("3a. Boxplots Audible saved to:", box_aud_path))
}

if(!is.null(res_us)) {
  box_us <- create_detailed_boxplots(res_us, "ULTRASOUND")
  box_us_path <- paste0(dirname(file_path), "/Detailed_Boxplots_Ultrasound_240mm.pdf")
  ggsave(box_us_path, box_us, device = cairo_pdf, width = 240, height = 140, units = "mm")
  message(paste("3b. Boxplots Ultrasound saved to:", box_us_path))
}
```

# 4 OBJECTIVE 4: LANDSCAPE REPRESENTATION

Specific objective 4: Determine the influence of landscape heterogeneity on the performance and representativeness of soundscape information obtained using tessellation-based sampling strategies

```{r}
library(tidyverse)
library(readxl)
library(ggpubr)
library(openxlsx)  # For exporting to Excel

# Load Excel file
df <- read_excel("#/DF3_Subsampling_Heterogeneity_Comparison.xlsx")

# Filter valid values
df_filtered <- df %>% filter(!is.na(L_proxy_mean))

# Extract unique localities
localities <- unique(df_filtered$Locality)

# Set PDF output path
output_pdf <- "#/Comparative_Plots_Lproxy_mean_POLYPOLES_ALLTOGETHER.pdf"
pdf_opened <- FALSE

# Initialize results table
statistical_results <- data.frame(
  Locality = character(),
  Method = character(),
  Procedence = character(),
  KS_statistic = numeric(),
  KS_p_value = numeric(),
  Wasserstein = numeric(),
  stringsAsFactors = FALSE
)

for (loc in localities) {
  cat("\n\n📍 Locality:", loc, "\n")

  df_loc <- df_filtered %>% filter(Locality == loc)
  full_vals <- df_loc %>% filter(Tessell == "complete") %>% pull(L_proxy_mean)

  for (method in unique(df_loc$Tessell)) {
    if (method == "complete") next

    df_sub_all <- df_loc %>% filter(Tessell == method & Procedence %in% c("Examinated", "Poles"))
    if (nrow(df_sub_all) == 0) next

    # Compute statistics by procedence
    for (proc in c("Examinated", "Poles")) {
      df_sub <- df_loc %>% filter(Tessell == method, Procedence == proc)
      if (nrow(df_sub) == 0) next

      sub_vals <- df_sub$L_proxy_mean
      n <- min(length(full_vals), length(sub_vals))
      ks <- ks.test(full_vals, sub_vals)
      wass <- mean(abs(sort(full_vals)[1:n] - sort(sub_vals)[1:n]))

      cat("🔎 Method:", method, "| Procedence:", proc, "\n")
      cat("   KS p-value:", ks$p.value, "\n")
      cat("   KS statistic:", ks$statistic, "\n")
      cat("   Approx. Wasserstein distance:", round(wass, 4), "\n")

      statistical_results <- rbind(statistical_results, data.frame(
        Locality = loc,
        Method = method,
        Procedence = proc,
        KS_statistic = as.numeric(ks$statistic),
        KS_p_value = as.numeric(ks$p.value),
        Wasserstein = round(wass, 4),
        stringsAsFactors = FALSE
      ))
    }

    # Generate density and Q–Q plots
    if (!pdf_opened) {
      pdf(output_pdf, width = 14, height = 6)
      pdf_opened <- TRUE
    }

    # Density plot
    p1 <- ggplot(
      df_loc %>% filter((Tessell == method & Procedence %in% c("Examinated", "Poles")) | Tessell == "complete"),
      aes(x = L_proxy_mean, fill = interaction(Tessell, Procedence, drop = TRUE))
    ) +
      geom_density(alpha = 0.4) +
      theme_minimal() +
      labs(title = "Density Plot",
           x = "L_proxy_mean", y = "Density", fill = "Tessellation / Procedence")

    # Q–Q plot
    qq_dfs <- list()
    for (proc in c("Examinated", "Poles")) {
      df_sub <- df_loc %>% filter(Tessell == method, Procedence == proc)
      if (nrow(df_sub) == 0) next

      sub_vals <- df_sub$L_proxy_mean
      n <- min(length(full_vals), length(sub_vals))

      qq_dfs[[proc]] <- data.frame(
        complete = sort(full_vals)[1:n],
        sub = sort(sub_vals)[1:n],
        Procedence = proc
      )
    }
    qq_df_all <- bind_rows(qq_dfs)

    p2 <- ggplot(qq_df_all, aes(x = complete, y = sub, color = Procedence)) +
      geom_point(alpha = 0.6) +
      geom_abline(slope = 1, intercept = 0, color = "black", linetype = "dashed") +
      theme_minimal() +
      labs(title = "Q–Q Plot", x = "Complete Distribution", y = method)

    combined <- ggarrange(p2, p1, ncol = 2, widths = c(1, 1),
                          labels = c("A", "B"),
                          common.legend = TRUE, legend = "right")

    combined <- annotate_figure(combined,
                                 top = text_grob(paste0(loc, " / ", method, " / Comparison with Examinated & Poles"),
                                                 face = "bold", size = 14))

    print(combined)
  }
}

if (pdf_opened) {
  dev.off()
  cat("\n✅ PDF successfully generated with all comparisons and statistics.\n")
} else {
  cat("\n⚠️ No plots were generated. Please check available data.\n")
}

# Save statistical results to Excel
output_excel <- "#/Comparative_Stats_Lproxy_mean_POLYPOLES_ALLTOGETHER.xlsx"
write.xlsx(statistical_results, output_excel)

```
