Getting Started with ecoGLMM

Overview

ecoGLMM provides a reproducible workflow for fitting and comparing ecological generalized linear mixed models. For each response and environmental predictor, the package can fit an additive model and a model containing an interaction with a temporal factor. Candidate models are compared using AICc, and the selected models can be summarized, diagnosed, plotted, and exported.

This vignette presents a complete analysis with simulated data. The lightweight model example is evaluated when the vignette is built, whereas computationally intensive diagnostics and file exports are displayed without being run. Users can copy and run every example in an interactive R session.

Install and load the package

# Install the CRAN version when available:
install.packages("ecoGLMM")

# Install the development version from GitHub:
# install.packages("remotes")
remotes::install_github("andre-fcsantos/ecoGLMM", upgrade = "never")

library(ecoGLMM)

Required data structure

The input must be a data frame with:

one or more numeric response columns;
numeric environmental predictors;
a temporal factor, such as sampling period;
a grouping variable for the random intercept, such as site.

The final analysis data must contain all these columns in the same data frame. The generic template can import either one analysis-ready table or separate acoustic and environmental tables. In separate mode, it validates the join keys and combines the tables without changing the number of acoustic observations.

Each row of the resulting data set represents one observation. The example below creates 150 observations from 10 sites and five sampling periods.

set.seed(123)

sites <- paste0("Site_", seq_len(10))
period_levels <- c("T-0", "T-1", "T-2", "T-3", "T-4")

example_data <- expand.grid(
  site = sites,
  period = period_levels,
  replicate = seq_len(3),
  stringsAsFactors = FALSE
)

n <- nrow(example_data)
example_data$forest <- runif(n, 10, 95)
example_data$urban <- runif(n, 0, 60)

site_effect <- rnorm(length(sites), 0, 0.25)
names(site_effect) <- sites

linear_predictor <-
  -0.5 +
  0.015 * example_data$forest -
  0.012 * example_data$urban +
  site_effect[example_data$site]

example_data$acoustic_index <- plogis(
  linear_predictor + rnorm(n, 0, 0.35)
)

example_data$acoustic_index <- adjust_beta(example_data$acoustic_index)

head(example_data)
#>     site period replicate   forest     urban acoustic_index
#> 1 Site_1    T-0         1 34.44409 50.847190      0.4929702
#> 2 Site_2    T-0         1 77.00594 29.851636      0.5304418
#> 3 Site_3    T-0         1 44.76304 23.274542      0.3855983
#> 4 Site_4    T-0         1 85.05648 14.786940      0.8146743
#> 5 Site_5    T-0         1 89.93972  6.665788      0.6441318
#> 6 Site_6    T-0         1 13.87230 23.399666      0.3685527

Configure the responses

The configuration object requires the columns response and family. An optional label column provides a readable name for tables and figures.

Supported family names are "beta", "gamma", "gaussian", "poisson", and "nbinom2".

response_config <- data.frame(
  response = "acoustic_index",
  family = "beta",
  label = "Simulated acoustic index"
)

For ecoacoustic analyses, prepare_acoustic_indices() can transform NDSI from the interval [-1, 1] to (0, 1) and move boundary values away from zero and one before beta regression.

# analysis_data <- prepare_acoustic_indices(
#   data = raw_data,
#   ndsi = "NDSI",
#   beta_responses = c("AEI", "ACT"),
#   ndsi_output = "NDSI_beta"
# )

Run the candidate-model analysis

Here, additive and period-interaction models are fitted separately for forest and urban. Predictors are standardized by default. All candidate models use the same complete-case data set, allowing valid AICc comparisons.

fit <- run_ecoglmm(
  data = example_data,
  config = response_config,
  predictors = c("forest", "urban"),
  period = "period",
  group = "site",
  period_levels = period_levels,
  standardize = TRUE,
  complete_cases = TRUE,
  include_additive = TRUE,
  include_interactions = TRUE,
  competitive_delta = 2
)

fit
#> ecoGLMM analysis
#> Responses: 1 
#> Observations: 150 
#> 
#>        Response      Model Predictor     Type      AICc Akaike_weight
#>  acoustic_index forest_Add    forest Additive -250.0743     0.9403152
summary(fit)
#> $overview
#>         Response      Model Predictor     Type      AICc Akaike_weight
#> 1 acoustic_index forest_Add    forest Additive -250.0743     0.9403152
#> 
#> $coefficients
#>         Response      Model     Type Predictor        Term     Estimate
#> 1 acoustic_index forest_Add Additive    forest (Intercept)  0.043224004
#> 2 acoustic_index forest_Add Additive    forest      forest  0.363635345
#> 3 acoustic_index forest_Add Additive    forest   periodT-1 -0.048800242
#> 4 acoustic_index forest_Add Additive    forest   periodT-2 -0.105173540
#> 5 acoustic_index forest_Add Additive    forest   periodT-3 -0.002203573
#> 6 acoustic_index forest_Add Additive    forest   periodT-4 -0.086877852
#>           SE   Statistic      p_value     CI_low    CI_high R2_marginal
#> 1 0.09111297  0.47440014 6.352146e-01 -0.1353541 0.22180214   0.6738737
#> 2 0.03457496 10.51730269 7.190255e-26  0.2958697 0.43140103   0.6738737
#> 3 0.10258905 -0.47568667 6.342976e-01 -0.2498711 0.15227060   0.6738737
#> 4 0.10260697 -1.02501359 3.053568e-01 -0.3062795 0.09593243   0.6738737
#> 5 0.10240943 -0.02151729 9.828330e-01 -0.2029224 0.19851522   0.6738737
#> 6 0.10406477 -0.83484401 4.038055e-01 -0.2908411 0.11708536   0.6738737
#>   R2_conditional R2_marginal_performance R2_conditional_performance      AICc
#> 1      0.8207683               0.6738737                  0.8207683 -250.0743
#> 2      0.8207683               0.6738737                  0.8207683 -250.0743
#> 3      0.8207683               0.6738737                  0.8207683 -250.0743
#> 4      0.8207683               0.6738737                  0.8207683 -250.0743
#> 5      0.8207683               0.6738737                  0.8207683 -250.0743
#> 6      0.8207683               0.6738737                  0.8207683 -250.0743
#>   Delta_AICc Akaike_weight
#> 1          0     0.9403152
#> 2          0     0.9403152
#> 3          0     0.9403152
#> 4          0     0.9403152
#> 5          0     0.9403152
#> 6          0     0.9403152
#> 
#> $interaction_tests
#> $interaction_tests$acoustic_index
#>   Predictor    Chisq df   p_value
#> 1    forest 3.741820  4 0.4420726
#> 2     urban 3.519027  4 0.4749912
#> 
#> 
#> $preparation
#> $preparation$removed_rows
#> [1] 0
#> 
#> $preparation$original_n
#> [1] 150
#> 
#> $preparation$analysis_n
#> [1] 150
#> 
#> $preparation$standardize
#> [1] TRUE
#> 
#> $preparation$scaling
#> $preparation$scaling$forest
#>   center    scale 
#> 52.85439 24.40825 
#> 
#> $preparation$scaling$urban
#>   center    scale 
#> 29.70327 16.48889

Inspect selection and coefficients

fit$selection$acoustic_index
#>        Model Predictor        Type   n  K      AICc Delta_AICc Akaike_weight
#> 1 forest_Add    forest    Additive 150  8 -250.0743   0.000000  9.403152e-01
#> 2 forest_Int    forest Interaction 150 12 -244.5600   5.514275  5.968481e-02
#> 3  urban_Add     urban    Additive 150  8 -200.0745  49.999795  1.306038e-11
#> 4  urban_Int     urban Interaction 150 12 -194.3374  55.736864  7.415957e-13
#>   Evidence_ratio Competitive
#> 1   1.000000e+00        TRUE
#> 2   1.575468e+01       FALSE
#> 3   7.199752e+10       FALSE
#> 4   1.267962e+12       FALSE
fit$interaction_tests$acoustic_index
#>   Predictor    Chisq df   p_value
#> 1    forest 3.741820  4 0.4420726
#> 2     urban 3.519027  4 0.4749912
fit$coefficients
#>         Response      Model     Type Predictor        Term     Estimate
#> 1 acoustic_index forest_Add Additive    forest (Intercept)  0.043224004
#> 2 acoustic_index forest_Add Additive    forest      forest  0.363635345
#> 3 acoustic_index forest_Add Additive    forest   periodT-1 -0.048800242
#> 4 acoustic_index forest_Add Additive    forest   periodT-2 -0.105173540
#> 5 acoustic_index forest_Add Additive    forest   periodT-3 -0.002203573
#> 6 acoustic_index forest_Add Additive    forest   periodT-4 -0.086877852
#>           SE   Statistic      p_value     CI_low    CI_high R2_marginal
#> 1 0.09111297  0.47440014 6.352146e-01 -0.1353541 0.22180214   0.6738737
#> 2 0.03457496 10.51730269 7.190255e-26  0.2958697 0.43140103   0.6738737
#> 3 0.10258905 -0.47568667 6.342976e-01 -0.2498711 0.15227060   0.6738737
#> 4 0.10260697 -1.02501359 3.053568e-01 -0.3062795 0.09593243   0.6738737
#> 5 0.10240943 -0.02151729 9.828330e-01 -0.2029224 0.19851522   0.6738737
#> 6 0.10406477 -0.83484401 4.038055e-01 -0.2908411 0.11708536   0.6738737
#>   R2_conditional R2_marginal_performance R2_conditional_performance      AICc
#> 1      0.8207683               0.6738737                  0.8207683 -250.0743
#> 2      0.8207683               0.6738737                  0.8207683 -250.0743
#> 3      0.8207683               0.6738737                  0.8207683 -250.0743
#> 4      0.8207683               0.6738737                  0.8207683 -250.0743
#> 5      0.8207683               0.6738737                  0.8207683 -250.0743
#> 6      0.8207683               0.6738737                  0.8207683 -250.0743
#>   Delta_AICc Akaike_weight
#> 1          0     0.9403152
#> 2          0     0.9403152
#> 3          0     0.9403152
#> 4          0     0.9403152
#> 5          0     0.9403152
#> 6          0     0.9403152
fit$preparation
#> $removed_rows
#> [1] 0
#> 
#> $original_n
#> [1] 150
#> 
#> $analysis_n
#> [1] 150
#> 
#> $standardize
#> [1] TRUE
#> 
#> $scaling
#> $scaling$forest
#>   center    scale 
#> 52.85439 24.40825 
#> 
#> $scaling$urban
#>   center    scale 
#> 29.70327 16.48889

A model with the smallest AICc is ranked first. Models with Delta_AICc <= 2 are marked as competitive. Akaike weights express relative support within the candidate set. The likelihood-ratio tests compare the nested additive and interaction models for each predictor.

The columns R2_marginal and R2_conditional are calculated with MuMIn::r.squaredGLMM(). Alternative estimates from performance::r2_nakagawa() are provided in columns ending in _performance.

Run residual diagnostics

diagnostics <- diagnose_models(
  fit,
  nsim = 1000,
  seed = 123
)

diagnostics$summary

# Inspect the DHARMa residual plot:
plot(diagnostics$residuals$acoustic_index)

The summary reports uniformity, dispersion, and outlier tests. A "Passed" status means that all three p-values exceed 0.05; it does not replace graphical inspection or ecological judgment.

Create figures

plot_effect(fit, "acoustic_index")
#> You are calculating adjusted predictions on the population-level (i.e.
#>   `type = "fixed"`) for a *generalized* linear mixed model.
#>   This may produce biased estimates due to Jensen's inequality. Consider
#>   setting `bias_correction = TRUE` to correct for this bias.
#>   See also the documentation of the `bias_correction` argument.

plot_selection(fit, "acoustic_index", metric = "delta")

plot_selection(fit, "acoustic_index", metric = "weight")

plot_coefficients(fit)

Export the analysis

exported_files <- export_ecoglmm(
  object = fit,
  path = file.path(tempdir(), "ecoGLMM_results"),
  diagnostics = diagnostics,
  figures = TRUE
)

exported_files

The output directory contains an Excel workbook with the overview, coefficients, model-selection results, interaction tests, and diagnostics. When figures = TRUE, publication-quality PNG files are also created.

Use the complete generic analysis template

The installed package includes a complete generic template that accepts either one analysis-ready table or separate acoustic and environmental tables in CSV or Excel format. When separate files are used, the template joins them by a user-defined site key and checks for duplicated or unmatched values. The example uses fictional column names and does not contain data or settings from any unpublished study.

Locate, copy, and open an editable version of the template with:

template_path <- system.file(
  "examples",
  "ecoGLMM_template.R",
  package = "ecoGLMM"
)

file.copy(
  from = template_path,
  to = file.path(tempdir(), "ecoGLMM_template.R"),
  overwrite = FALSE
)

file.edit(file.path(tempdir(), "ecoGLMM_template.R"))

For a persistent copy, choose your own writable destination instead of tempdir(). Edit only Section 01 of the copied file before running the complete script. The full template is reproduced below directly from the file distributed with the package, ensuring that the vignette and script remain synchronized.

###############################################################
## ecoGLMM — GENERIC ANALYSIS TEMPLATE
##
## Edit only Section 01 before running the complete script.
## Input: either one analysis-ready table or separate acoustic
## and environmental tables in CSV or Excel format.
## Each row of the final data set must represent one observation.
###############################################################

###############################################################
## 01. USER SETTINGS — EDIT THIS SECTION
###############################################################

# Choose "single" for one analysis-ready table or "separate" to
# join acoustic and environmental tables using the configured keys.
INPUT_MODE <- "single"

SINGLE_DATA_FILE <- NULL
ACOUSTIC_DATA_FILE <- NULL
ENVIRONMENTAL_DATA_FILE <- NULL

SINGLE_EXCEL_SHEET <- 1
ACOUSTIC_EXCEL_SHEET <- 1
ENVIRONMENTAL_EXCEL_SHEET <- 1

# These may be identical. Separate names are supported when the
# site identifier differs between the two input tables.
ACOUSTIC_JOIN_COLUMN <- "site"
ENVIRONMENTAL_JOIN_COLUMN <- "site"

GROUP_COLUMN <- "site"
PERIOD_COLUMN <- "period"
PERIOD_LEVELS <- c("T-0", "T-1", "T-2", "T-3", "T-4")
PREDICTORS <- c("forest", "urban")

# Set to TRUE when NDSI must be converted from [-1, 1] to (0, 1)
# and other beta-distributed responses must be moved away from
# exact zero and one.
PREPARE_ACOUSTIC_INDICES <- FALSE
NDSI_COLUMN <- "NDSI"
NDSI_OUTPUT_COLUMN <- "NDSI_beta"
BETA_RESPONSE_COLUMNS <- c("AEI", "ACT")
BETA_EPSILON <- 0.001

# Optional cleaning of non-analytical columns after import.
DROP_EMPTY_COLUMNS <- TRUE
DROP_COLUMNS <- c("-")

# One row per response. Supported families:
# "beta", "gamma", "gaussian", "poisson", and "nbinom2".
RESPONSE_CONFIG <- data.frame(
  response = c("acoustic_index"),
  family = c("beta"),
  label = c("Acoustic index"),
  stringsAsFactors = FALSE
)

STANDARDIZE_PREDICTORS <- TRUE
COMPLETE_CASES <- TRUE
COMPETITIVE_DELTA <- 2
RUN_DIAGNOSTICS <- TRUE
DHARMA_SIMULATIONS <- 1000
RANDOM_SEED <- 123
# Set an explicit output directory for permanent results. The default
# uses a temporary directory and will not persist after this R session.
OUTPUT_DIRECTORY <- file.path(tempdir(), "ecoGLMM_results")

###############################################################
## 02. PACKAGE AND DATA IMPORT
###############################################################

if (!requireNamespace("ecoGLMM", quietly = TRUE)) {
  stop(
    "Install ecoGLMM before running this script. ",
    "See https://github.com/andre-fcsantos/ecoGLMM",
    call. = FALSE
  )
}

read_input_table <- function(path, excel_sheet = 1) {
  extension <- tolower(tools::file_ext(path))

  if (extension %in% c("xlsx", "xls")) {
    if (!requireNamespace("readxl", quietly = TRUE)) {
      stop("Install the readxl package to import Excel files.", call. = FALSE)
    }
    output <- readxl::read_excel(path, sheet = excel_sheet)
  } else if (extension == "csv") {
    output <- utils::read.csv(
      path,
      stringsAsFactors = FALSE,
      check.names = FALSE
    )
  } else {
    stop("Use an .xlsx, .xls, or .csv input file.", call. = FALSE)
  }

  as.data.frame(output)
}

INPUT_MODE <- match.arg(tolower(INPUT_MODE), c("single", "separate"))

if (INPUT_MODE == "single") {
  if (is.null(SINGLE_DATA_FILE)) {
    message("Select the analysis-ready data file.")
    SINGLE_DATA_FILE <- file.choose()
  }

  analysis_data <- read_input_table(
    SINGLE_DATA_FILE,
    excel_sheet = SINGLE_EXCEL_SHEET
  )

  cat(
    "Imported one analysis-ready table: ",
    basename(SINGLE_DATA_FILE),
    "\n",
    sep = ""
  )
} else {
  if (is.null(ACOUSTIC_DATA_FILE)) {
    message("Select the acoustic data file.")
    ACOUSTIC_DATA_FILE <- file.choose()
  }
  if (is.null(ENVIRONMENTAL_DATA_FILE)) {
    message("Select the environmental data file.")
    ENVIRONMENTAL_DATA_FILE <- file.choose()
  }

  acoustic_data <- read_input_table(
    ACOUSTIC_DATA_FILE,
    excel_sheet = ACOUSTIC_EXCEL_SHEET
  )
  environmental_data <- read_input_table(
    ENVIRONMENTAL_DATA_FILE,
    excel_sheet = ENVIRONMENTAL_EXCEL_SHEET
  )

  cat(
    "Imported acoustic table: ", basename(ACOUSTIC_DATA_FILE), "\n",
    "Imported environmental table: ", basename(ENVIRONMENTAL_DATA_FILE), "\n",
    sep = ""
  )

  if (!ACOUSTIC_JOIN_COLUMN %in% names(acoustic_data)) {
    stop(
      "The acoustic join column was not found: ",
      ACOUSTIC_JOIN_COLUMN,
      call. = FALSE
    )
  }
  if (!ENVIRONMENTAL_JOIN_COLUMN %in% names(environmental_data)) {
    stop(
      "The environmental join column was not found: ",
      ENVIRONMENTAL_JOIN_COLUMN,
      call. = FALSE
    )
  }
  if (anyNA(environmental_data[[ENVIRONMENTAL_JOIN_COLUMN]])) {
    stop("The environmental join column contains missing values.", call. = FALSE)
  }
  if (anyDuplicated(environmental_data[[ENVIRONMENTAL_JOIN_COLUMN]])) {
    stop(
      "Each join value must occur only once in the environmental table.",
      call. = FALSE
    )
  }

  overlapping_columns <- intersect(
    setdiff(names(acoustic_data), ACOUSTIC_JOIN_COLUMN),
    setdiff(names(environmental_data), ENVIRONMENTAL_JOIN_COLUMN)
  )

  if (length(overlapping_columns) > 0) {
    stop(
      "Non-key columns occur in both input tables: ",
      paste(overlapping_columns, collapse = ", "),
      ". Rename or remove them before joining.",
      call. = FALSE
    )
  }

  unmatched_values <- setdiff(
    unique(acoustic_data[[ACOUSTIC_JOIN_COLUMN]]),
    unique(environmental_data[[ENVIRONMENTAL_JOIN_COLUMN]])
  )

  if (length(unmatched_values) > 0) {
    warning(
      "Environmental data were not found for: ",
      paste(unmatched_values, collapse = ", "),
      call. = FALSE
    )
  }

  acoustic_n <- nrow(acoustic_data)
  row_order_column <- ".ecoGLMM_input_row"

  if (row_order_column %in% c(
    names(acoustic_data),
    names(environmental_data)
  )) {
    stop(
      "Reserved internal column found in an input table: ",
      row_order_column,
      call. = FALSE
    )
  }

  acoustic_data[[row_order_column]] <- seq_len(acoustic_n)

  analysis_data <- merge(
    acoustic_data,
    environmental_data,
    by.x = ACOUSTIC_JOIN_COLUMN,
    by.y = ENVIRONMENTAL_JOIN_COLUMN,
    all.x = TRUE,
    sort = FALSE
  )

  if (nrow(analysis_data) != acoustic_n) {
    stop(
      "The join changed the number of acoustic observations. ",
      "Check the configured join columns.",
      call. = FALSE
    )
  }

  analysis_data <- analysis_data[
    order(analysis_data[[row_order_column]]),
    ,
    drop = FALSE
  ]
  analysis_data[[row_order_column]] <- NULL

  cat(
    "Joined ", acoustic_n, " acoustic observations to ",
    nrow(environmental_data), " environmental records.\n",
    sep = ""
  )
}

if (DROP_EMPTY_COLUMNS) {
  empty_columns <- vapply(
    analysis_data,
    function(column) all(is.na(column)),
    logical(1)
  )
  analysis_data <- analysis_data[, !empty_columns, drop = FALSE]
}

columns_to_drop <- intersect(DROP_COLUMNS, names(analysis_data))
if (length(columns_to_drop) > 0) {
  analysis_data[columns_to_drop] <- NULL
}

structure_columns <- unique(c(
  PREDICTORS,
  PERIOD_COLUMN,
  GROUP_COLUMN
))

missing_structure_columns <- setdiff(
  structure_columns,
  names(analysis_data)
)

if (length(missing_structure_columns) > 0) {
  stop(
    "Columns not found in the imported data: ",
    paste(missing_structure_columns, collapse = ", "),
    call. = FALSE
  )
}

cat(
  "Prepared ", nrow(analysis_data), " rows and ",
  ncol(analysis_data), " columns.\n",
  "Available columns: ",
  paste(names(analysis_data), collapse = ", "),
  "\n",
  sep = ""
)

###############################################################
## 03. OPTIONAL RESPONSE TRANSFORMATIONS
###############################################################

# Configure this operation in Section 01. When enabled, NDSI is
# converted from [-1, 1] to (0, 1), and configured beta responses
# are moved away from exact boundary values.
if (PREPARE_ACOUSTIC_INDICES) {
  analysis_data <- ecoGLMM::prepare_acoustic_indices(
    data = analysis_data,
    ndsi = NDSI_COLUMN,
    beta_responses = BETA_RESPONSE_COLUMNS,
    ndsi_output = NDSI_OUTPUT_COLUMN,
    epsilon = BETA_EPSILON
  )
}

# Response validation occurs after optional transformations so newly
# created columns such as NDSI_beta are recognized correctly.
missing_responses <- setdiff(
  RESPONSE_CONFIG$response,
  names(analysis_data)
)

if (length(missing_responses) > 0) {
  transformation_hint <- ""

  if (
    NDSI_OUTPUT_COLUMN %in% missing_responses &&
    !PREPARE_ACOUSTIC_INDICES
  ) {
    transformation_hint <- paste0(
      " To create ", NDSI_OUTPUT_COLUMN,
      ", set PREPARE_ACOUSTIC_INDICES <- TRUE in Section 01."
    )
  }

  stop(
    "Response columns not found after data preparation: ",
    paste(missing_responses, collapse = ", "),
    ".",
    transformation_hint,
    call. = FALSE
  )
}

###############################################################
## 04. FIT AND COMPARE CANDIDATE MODELS
###############################################################

set.seed(RANDOM_SEED)

results <- ecoGLMM::run_ecoglmm(
  data = analysis_data,
  config = RESPONSE_CONFIG,
  predictors = PREDICTORS,
  period = PERIOD_COLUMN,
  group = GROUP_COLUMN,
  period_levels = PERIOD_LEVELS,
  standardize = STANDARDIZE_PREDICTORS,
  complete_cases = COMPLETE_CASES,
  include_additive = TRUE,
  include_interactions = TRUE,
  competitive_delta = COMPETITIVE_DELTA
)

print(results)
analysis_summary <- summary(results)

for (response in names(results$selection)) {
  cat("\nMODEL SELECTION — ", response, "\n", sep = "")
  print(results$selection[[response]])

  cat("\nINTERACTION TEST — ", response, "\n", sep = "")
  print(results$interaction_tests[[response]])
}

###############################################################
## 05. DIAGNOSTICS
###############################################################

diagnostics <- NULL

if (RUN_DIAGNOSTICS) {
  diagnostics <- ecoGLMM::diagnose_models(
    object = results,
    nsim = DHARMA_SIMULATIONS,
    seed = RANDOM_SEED
  )
  print(diagnostics$summary)

  # Example for graphical inspection:
  # plot(diagnostics$residuals[[RESPONSE_CONFIG$response[1]]])
}

###############################################################
## 06. FIGURES AND EXPORT
###############################################################

# Examples for the RStudio Plots pane:
# ecoGLMM::plot_effect(results, RESPONSE_CONFIG$response[1])
# ecoGLMM::plot_selection(
#   results,
#   RESPONSE_CONFIG$response[1],
#   metric = "delta"
# )
# ecoGLMM::plot_coefficients(results)

dir.create(OUTPUT_DIRECTORY, recursive = TRUE, showWarnings = FALSE)

generated_files <- ecoGLMM::export_ecoglmm(
  object = results,
  path = OUTPUT_DIRECTORY,
  diagnostics = diagnostics,
  figures = TRUE
)

saveRDS(
  results,
  file.path(OUTPUT_DIRECTORY, "ecoGLMM_results.rds")
)

if (!is.null(diagnostics)) {
  saveRDS(
    diagnostics,
    file.path(OUTPUT_DIRECTORY, "DHARMa_diagnostics.rds")
  )
}

capture.output(
  sessionInfo(),
  file = file.path(OUTPUT_DIRECTORY, "sessionInfo.txt")
)

cat(
  "\nAnalysis completed. Results were saved to:\n",
  normalizePath(OUTPUT_DIRECTORY, winslash = "/", mustWork = FALSE),
  "\n",
  sep = ""
)

Users should adapt file paths, column mappings, responses, distributions, predictors, temporal levels, output folders, and analysis options only in their own copy. After changing any setting, run the script again from Section 01.

For two input tables, set INPUT_MODE <- "separate" and provide the exact join-column name used by each table. The template verifies uniqueness in the environmental table, reports unmatched keys, preserves the original acoustic-row order, removes configured non-analytical columns, and confirms the final column names.

For ecoacoustic data requiring NDSI and beta-response preparation, set PREPARE_ACOUSTIC_INDICES <- TRUE in Section 01. The corresponding column names and boundary adjustment can be configured without editing the processing code. If the option is left disabled while the configured NDSI output is requested, the template returns a targeted instruction instead of a generic missing-column error.

Reproducibility recommendations

Preserve the original data and analyse a separate working copy.
Record the response distribution selected for each ecological index.
Keep complete_cases = TRUE when comparing models by AICc.
Inspect convergence messages and DHARMa plots, not only p-values.
Report the candidate set, AICc, delta AICc, Akaike weights, coefficients, confidence intervals, and marginal and conditional R-squared.
Save sessionInfo() with the exported results.