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
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.3685527Configure 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.
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.48889Inspect 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.48889A 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.
Export the analysis
exported_files <- export_ecoglmm(
object = fit,
path = file.path(tempdir(), "ecoGLMM_results"),
diagnostics = diagnostics,
figures = TRUE
)
exported_filesThe 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. When a diagnostics object is supplied,
DHARMa_diagnostics.pdf is generated with the residual
panels and a formal-test summary for every selected response model.
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"
)
if (!nzchar(template_path) || !file.exists(template_path)) {
stop("The generic analysis template was not found in the installed package.")
}
template_copy <- file.path(tempdir(), "ecoGLMM_template.R")
copied <- file.copy(
from = template_path,
to = template_copy,
overwrite = TRUE
)
if (!copied || !file.exists(template_copy) || file.size(template_copy) == 0) {
stop("The generic analysis template could not be copied correctly.")
}
if (interactive()) {
file.edit(template_copy)
}
template_copyFor 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")
# Because ecoGLMM fits one environmental predictor per candidate model,
# classical VIF is not required for the current model structure. The optional
# correlation screen below identifies potentially redundant ecological
# hypotheses without removing variables automatically.
CHECK_PREDICTOR_CORRELATIONS <- TRUE
CORRELATION_METHOD <- "spearman"
CORRELATION_THRESHOLD <- 0.70
# 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. DATA-QUALITY AND PREDICTOR-REDUNDANCY CHECKS
###############################################################
required_config_columns <- c("response", "family")
missing_config_columns <- setdiff(
required_config_columns,
names(RESPONSE_CONFIG)
)
if (length(missing_config_columns) > 0) {
stop(
"RESPONSE_CONFIG is missing required columns: ",
paste(missing_config_columns, collapse = ", "),
call. = FALSE
)
}
RESPONSE_CONFIG$response <- as.character(RESPONSE_CONFIG$response)
RESPONSE_CONFIG$family <- tolower(as.character(RESPONSE_CONFIG$family))
if (
anyNA(RESPONSE_CONFIG$response) ||
any(!nzchar(RESPONSE_CONFIG$response))
) {
stop("RESPONSE_CONFIG contains an empty response name.", call. = FALSE)
}
if (anyDuplicated(RESPONSE_CONFIG$response)) {
stop("Each response must occur only once in RESPONSE_CONFIG.", call. = FALSE)
}
supported_families <- c(
"beta", "gamma", "gaussian", "poisson", "nbinom2"
)
unsupported_families <- setdiff(
unique(RESPONSE_CONFIG$family),
supported_families
)
if (length(unsupported_families) > 0) {
stop(
"Unsupported families in RESPONSE_CONFIG: ",
paste(unsupported_families, collapse = ", "),
call. = FALSE
)
}
if (anyDuplicated(names(analysis_data))) {
duplicated_names <- unique(names(analysis_data)[duplicated(names(analysis_data))])
stop(
"Duplicated column names were found: ",
paste(duplicated_names, collapse = ", "),
call. = FALSE
)
}
numeric_columns <- unique(c(
PREDICTORS,
RESPONSE_CONFIG$response
))
for (column_name in numeric_columns) {
if (!is.numeric(analysis_data[[column_name]])) {
original_values <- analysis_data[[column_name]]
converted_values <- suppressWarnings(as.numeric(gsub(
",",
".",
as.character(original_values),
fixed = TRUE
)))
unsafe_values <- is.na(converted_values) & !is.na(original_values)
if (any(unsafe_values)) {
stop(
"Column could not be safely converted to numeric: ",
column_name,
call. = FALSE
)
}
analysis_data[[column_name]] <- converted_values
}
non_finite <- !is.na(analysis_data[[column_name]]) &
!is.finite(analysis_data[[column_name]])
if (any(non_finite)) {
stop(
"Column contains infinite or non-finite values: ",
column_name,
call. = FALSE
)
}
}
constant_predictors <- PREDICTORS[vapply(
analysis_data[PREDICTORS],
function(column) length(unique(stats::na.omit(column))) < 2,
logical(1)
)]
if (length(constant_predictors) > 0) {
stop(
"Predictors with fewer than two observed values: ",
paste(constant_predictors, collapse = ", "),
call. = FALSE
)
}
observed_groups <- unique(stats::na.omit(analysis_data[[GROUP_COLUMN]]))
observed_periods <- unique(stats::na.omit(analysis_data[[PERIOD_COLUMN]]))
if (length(observed_groups) < 2) {
stop("At least two grouping levels are required.", call. = FALSE)
}
if (length(observed_periods) < 2) {
stop("At least two observed period levels are required.", call. = FALSE)
}
if (!is.null(PERIOD_LEVELS)) {
unknown_periods <- setdiff(observed_periods, PERIOD_LEVELS)
if (length(unknown_periods) > 0) {
stop(
"Observed period values are absent from PERIOD_LEVELS: ",
paste(unknown_periods, collapse = ", "),
call. = FALSE
)
}
}
for (row_index in seq_len(nrow(RESPONSE_CONFIG))) {
response_name <- RESPONSE_CONFIG$response[row_index]
response_family <- RESPONSE_CONFIG$family[row_index]
response_values <- stats::na.omit(analysis_data[[response_name]])
if (!length(response_values)) {
stop("Response contains no observed values: ", response_name, call. = FALSE)
}
invalid_values <- switch(
response_family,
beta = response_values <= 0 | response_values >= 1,
gamma = response_values <= 0,
poisson = response_values < 0 | response_values != floor(response_values),
nbinom2 = response_values < 0 | response_values != floor(response_values),
gaussian = rep(FALSE, length(response_values))
)
if (any(invalid_values)) {
family_requirement <- switch(
response_family,
beta = "values strictly between zero and one",
gamma = "strictly positive values",
poisson = "non-negative integer counts",
nbinom2 = "non-negative integer counts"
)
stop(
"Response `", response_name, "` is configured as ",
response_family, " but contains invalid values. It requires ",
family_requirement, ".",
call. = FALSE
)
}
}
analysis_columns <- unique(c(
RESPONSE_CONFIG$response,
PREDICTORS,
PERIOD_COLUMN,
GROUP_COLUMN
))
column_role <- function(column_name) {
roles <- character(0)
if (column_name %in% RESPONSE_CONFIG$response) roles <- c(roles, "response")
if (column_name %in% PREDICTORS) roles <- c(roles, "predictor")
if (identical(column_name, PERIOD_COLUMN)) roles <- c(roles, "period")
if (identical(column_name, GROUP_COLUMN)) roles <- c(roles, "group")
paste(roles, collapse = "; ")
}
data_quality <- data.frame(
Variable = analysis_columns,
Role = vapply(analysis_columns, column_role, character(1)),
Class = vapply(
analysis_data[analysis_columns],
function(column) paste(class(column), collapse = "/"),
character(1)
),
Non_missing = vapply(
analysis_data[analysis_columns],
function(column) sum(!is.na(column)),
integer(1)
),
Missing = vapply(
analysis_data[analysis_columns],
function(column) sum(is.na(column)),
integer(1)
),
Unique = vapply(
analysis_data[analysis_columns],
function(column) length(unique(stats::na.omit(column))),
integer(1)
),
stringsAsFactors = FALSE,
row.names = NULL
)
complete_analysis_rows <- stats::complete.cases(
analysis_data[analysis_columns]
)
if (!any(complete_analysis_rows)) {
stop(
"No complete observations remain across the configured analysis columns.",
call. = FALSE
)
}
cat("\nDATA-QUALITY SUMMARY\n")
print(data_quality, row.names = FALSE)
cat(
"Complete observations available to the shared candidate set: ",
sum(complete_analysis_rows), " of ", nrow(analysis_data), "\n",
"Grouping levels: ", length(observed_groups), "\n",
"Observed period levels: ", length(observed_periods), "\n",
sep = ""
)
predictor_correlations <- data.frame(
Predictor_1 = character(0),
Predictor_2 = character(0),
Correlation = numeric(0),
Absolute_correlation = numeric(0),
Above_threshold = logical(0),
stringsAsFactors = FALSE
)
if (CHECK_PREDICTOR_CORRELATIONS && length(PREDICTORS) > 1) {
CORRELATION_METHOD <- match.arg(
tolower(CORRELATION_METHOD),
c("pearson", "spearman", "kendall")
)
if (
length(CORRELATION_THRESHOLD) != 1 ||
!is.finite(CORRELATION_THRESHOLD) ||
CORRELATION_THRESHOLD < 0 ||
CORRELATION_THRESHOLD > 1
) {
stop("CORRELATION_THRESHOLD must be between zero and one.", call. = FALSE)
}
correlation_matrix <- stats::cor(
analysis_data[PREDICTORS],
use = "pairwise.complete.obs",
method = CORRELATION_METHOD
)
pair_indices <- which(upper.tri(correlation_matrix), arr.ind = TRUE)
predictor_correlations <- data.frame(
Predictor_1 = rownames(correlation_matrix)[pair_indices[, "row"]],
Predictor_2 = colnames(correlation_matrix)[pair_indices[, "col"]],
Correlation = correlation_matrix[pair_indices],
stringsAsFactors = FALSE
)
predictor_correlations$Absolute_correlation <- abs(
predictor_correlations$Correlation
)
predictor_correlations$Above_threshold <- !is.na(
predictor_correlations$Absolute_correlation
) & predictor_correlations$Absolute_correlation >= CORRELATION_THRESHOLD
predictor_correlations <- predictor_correlations[
order(-predictor_correlations$Absolute_correlation, na.last = TRUE),
,
drop = FALSE
]
rownames(predictor_correlations) <- NULL
cat(
"\nPREDICTOR CORRELATIONS — ",
toupper(CORRELATION_METHOD), "\n",
sep = ""
)
print(predictor_correlations, row.names = FALSE)
flagged_pairs <- predictor_correlations[
predictor_correlations$Above_threshold,
,
drop = FALSE
]
if (nrow(flagged_pairs) > 0) {
warning(
nrow(flagged_pairs),
" predictor pair(s) reached |correlation| >= ",
CORRELATION_THRESHOLD,
". ecoGLMM fits one environmental predictor per candidate model, so ",
"classical VIF is not required. Review whether the flagged predictors ",
"represent redundant ecological hypotheses; no variables were removed.",
call. = FALSE
)
}
}
preflight_checks <- list(
data_quality = data_quality,
predictor_correlations = predictor_correlations,
complete_observations = sum(complete_analysis_rows),
total_observations = nrow(analysis_data),
grouping_levels = length(observed_groups),
observed_period_levels = length(observed_periods),
correlation_method = if (CHECK_PREDICTOR_CORRELATIONS) {
CORRELATION_METHOD
} else {
NA_character_
},
correlation_threshold = CORRELATION_THRESHOLD
)
###############################################################
## 05. 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]])
}
###############################################################
## 06. 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]]])
}
###############################################################
## 07. 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)
data_quality_file <- file.path(OUTPUT_DIRECTORY, "data_quality.csv")
correlation_file <- file.path(
OUTPUT_DIRECTORY,
"predictor_correlations.csv"
)
preflight_file <- file.path(OUTPUT_DIRECTORY, "preflight_checks.rds")
utils::write.csv(data_quality, data_quality_file, row.names = FALSE)
utils::write.csv(
predictor_correlations,
correlation_file,
row.names = FALSE
)
saveRDS(preflight_checks, preflight_file)
generated_files <- ecoGLMM::export_ecoglmm(
object = results,
path = OUTPUT_DIRECTORY,
diagnostics = diagnostics,
figures = TRUE,
diagnostic_pdf = TRUE
)
generated_files <- c(
generated_files,
data_quality = data_quality_file,
predictor_correlations = correlation_file,
preflight_checks = preflight_file
)
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.
Preflight checks and predictor redundancy
Before fitting any model, the generic template checks that configured response and predictor columns are numeric and finite, that predictors are not constant, and that group and period variables contain at least two observed levels. It also verifies that response values satisfy the selected distribution: beta responses must lie strictly between zero and one, Gamma responses must be positive, and Poisson or negative-binomial responses must be non-negative integers.
When CHECK_PREDICTOR_CORRELATIONS <- TRUE, the
template calculates all pairwise correlations among the environmental
predictors. The method and warning threshold are controlled in Section
01:
This screen is descriptive. It warns when two predictors may represent redundant ecological hypotheses, but it never removes variables or changes the candidate set. Classical variance inflation factors (VIF) are not calculated because ecoGLMM currently fits one environmental predictor per candidate model; VIF diagnoses collinearity among multiple predictors included in the same model. If multivariable candidate models are added in a future release, VIF or a comparable conditional-collinearity diagnostic should be assessed for those models.
The template saves data_quality.csv,
predictor_correlations.csv, and
preflight_checks.rds beside the regular model outputs so
the checks can be reported and reproduced.
Reproducibility recommendations
- Preserve the original data and analyse a separate working copy.
- Record the response distribution selected for each ecological index.
- Keep
complete_cases = TRUEwhen 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.
- Overview
- Install and load the package
- Required data structure
- Configure the responses
- Run the candidate-model analysis
- Inspect selection and coefficients
- Run residual diagnostics
- Create figures
- Export the analysis
- Use the complete generic analysis template
- Preflight checks and predictor redundancy
- Reproducibility recommendations