Package {SppTrend}

Title:	Analyzing Linear Trends in Species Occurrence Data
Version:	0.5
Description:	Provides a comparative framework to detect species-specific spatial and thermal responses to climate change using opportunistic occurrence data. Species temporal trends in geographic position (via Earth-Centred Earth-Fixed vector analysis) and environmental variables (temperature and elevation) are contrasted against the overall trend of the complete dataset, allowing classification of species into ecologically interpretable response categories. Approach described in Lobo et al. (2023) <doi:10.1002/ece3.10674>.
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.3
Suggests:	testthat (≥ 3.0.0), devtools
Imports:	data.table, dplyr, ggplot2, patchwork, sf, stats, stringr, terra
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2026-06-29 11:14:14 UTC; mario
Author:	Mario Mingarro [aut, cre], Emilio García-Roselló [aut], Jorge M. Lobo [aut]
Maintainer:	Mario Mingarro <mario_mingarro@mncn.csic.es>
Repository:	CRAN
Date/Publication:	2026-06-29 11:50:07 UTC

SppTrend: Analyzing Linear Trends in Species Occurrence Data

Description

Provides a comparative framework to detect species-specific spatial and thermal responses to climate change using opportunistic occurrence data. Species temporal trends in geographic position (via Earth-Centred Earth-Fixed vector analysis) and environmental variables (temperature and elevation) are contrasted against the overall trend of the complete dataset, allowing classification of species into ecologically interpretable response categories. Approach described in Lobo et al. (2023) doi:10.1002/ece3.10674.

Details

Methodology

SppTrend evaluates species-specific responses to climate change by contrasting individual temporal trends against the overall trend estimated from the complete dataset. This comparative approach explicitly accounts for sampling bias shared across species within a taxonomic group.

Species responses are decomposed into two complementary dimensions:

• Spatial: Temporal shifts in geographic position (latitude and longitude), analysed jointly via Earth-Centred Earth-Fixed (ECEF) vector analysis.

• Environmental: Temporal changes in temperature and elevation conditions associated with species occurrences.

Response categories

Spatial (from spp_trend_spatial):

• SA (Spatial Adaptation): movement faster than or directionally different from the global pool, in a poleward direction.

• SD (Spatial Discordance): movement slower than the global pool, or in a non-poleward direction.

• SC (Spatial Conformance): speed and direction indistinguishable from the global pattern.

Environmental (from spp_trend_environmental):

• TT (Thermal Tolerance): significant positive trend in temperature.

• TA (Thermal Adjustment): significant negative trend in temperature.

• TC (Thermal Conformance): temperature trend not significantly different from the overall trend.

• SA/SD/SC: equivalent classes applied to elevation trends.

Workflow

SppTrend provides a structured workflow for analysing these trends:

1. Rapid diagnostic: Quick visual check of spatial distribution and temperature trends using get_fast_info.

2. Environmental data integration (optional): Add ERA5 temperature values with get_era5_tme or elevation with get_elevation.

3. Spatial trend analysis: Estimate temporal changes in species geographic position using ECEF vector analysis with spp_trend_spatial.

4. Environmental trend analysis: Estimate temporal trends in temperature and elevation for each species using spp_trend_environmental.

More details

Source code: https://github.com/MarioMingarro/SppTrend

Author(s)

Maintainer: Mario Mingarro mario_mingarro@mncn.csic.es (ORCID)

Authors:

• Emilio García-Roselló (ORCID)

• Jorge M. Lobo (ORCID)

Extract elevation from DEM

Description

This function retrieves elevation values from a Digital Elevation Model (DEM) based on their geographic coordinates (lon/lat).

Usage

get_elevation(data, dem_file)

Arguments

Value

The input data frame data with a new column (ele) containing the extracted elevation values.

Extract temperature data from ERA5 NetCDF file

Description

This function retrieves mean monthly air temperature values associated with species occurrence records based on their geographic coordinates (lon/lat) and sampling date (year/month).

Usage

get_era5_tme(data, nc_file)

Arguments

Value

The input data frame data with a new column named (tme), containing the temperature values.

Quick visual diagnostic of the input data

Description

This function provides a quick visual diagnostic of the input data. It generates a map showing the spatial distribution of occurrence records together with a time-series plot derived from a NetCDF environmental dataste, including a linear trend analysis. Using the geographic coordinates of the occurrence records, the function extracts the complete climate time-series (from the earliest to the latest year represented in the data) for the corresponding occupied cells. All temperature values from occupied cells are then added annually to estimate and visualise the overall temperature trend (including slope and associated p-value). This diagnostic step allows users to quickly assess the climate trajectory of the regions where the species have been recorded and to evaluate whether sufficient temporal and environmental variation is present for subsequent analyses.

Usage

get_fast_info(data, nc_file)

Arguments

Value

Invisibly returns a composite plot. Displays a composite plot showing the geographic distribution and the thermal trend with its corresponding global slope and p-value.

Analyse Long-Term Environmental Niche Trends Across Species

Description

spp_trend_environmental() models temporal trends in environmental variables - Temperature and/or Elevation - for each species in a community dataset and evaluates whether species-specific shifts differ from the pooled community trend.

For every requested response variable the function fits:

1. A global ordinary least-squares (OLS) regression pooled across all valid species, capturing the average community-level temporal trend.

2. An independent OLS regression per species, yielding species-specific temporal slopes together with their standard errors, t-statistics, and 95 \

Species-specific slopes are then compared against the global slope using Welch-Satterthwaite t-tests. Raw p-values from both the individual species regressions and the slope-difference tests are corrected for multiple testing with the Benjamini-Hochberg False Discovery Rate (FDR) procedure, applied separately for each response variable. An optional sample-size correction penalises the significance threshold for well-sampled species (N > 100).

Usage

spp_trend_environmental(
  data,
  spp = NULL,
  responses = c("Temperature", "Elevation"),
  min_records = 20,
  min_years = 5,
  alpha_ref = 0.05,
  use_Ncorrected_alpha = TRUE,
  individual_lm_min_thresholds = c(Temperature = 0.01, Elevation = 50),
  individual_lm_threshold_type = "speciesslope",
  verbose = TRUE
)

Arguments

Value

A named list with the following elements:

params

Named list storing the values of all function arguments (except data and spp) as supplied by the caller.Useful for reproducibility and logging.

species_filter

Data frame with one row per species x response combination, summarising the filtering step.Columns:

Species

Character.Species name.

response

Character.Response variable ("Temperature" or "Elevation").

n_records

Integer.Total number of non-missing records for this species and response.

n_years

Integer.Number of distinct calendar years with at least one non-missing record.

retained

Logical.TRUE if the species meets both min_records and min_years thresholds and was therefore included in the modelling step.

environmental_global_lm

Data frame with one row per response variable containing statistics from the pooled (community-level) OLS regression.Columns:

response

Character.Response variable.

model

Character.Always "global_lm".

global_lm_slope

Numeric.Estimated temporal slope of the global model (units per year: deg C yr^{-1} or m yr^{-1}).

global_lm_se

Numeric.Standard error of the global slope.

global_lm_t

Numeric.t-statistic for the global slope (H0: slope = 0).

global_lm_p_value

Numeric.Two-tailed p-value for the global slope test.

global_lm_conf_low

Numeric.Lower bound of the 95 \ confidence interval for the global slope.

global_lm_conf_high

Numeric.Upper bound of the 95 \ confidence interval for the global slope.

n_records

Integer.Total number of records used to fit the global model (all valid records across retained species).

n_years

Integer.Number of distinct calendar years represented in the global model dataset.

environmental_individual_lm

Data frame with one row per species x response combination for all species that passed the filtering step. Contains species-level OLS statistics plus comparisons against the global model.Columns:

Species

Character.Species name.

response

Character.Response variable.

n_records

Integer.Number of records used to fit this species model.

n_years

Integer.Number of distinct years in the species data for this response.

individual_slope

Numeric.Estimated temporal slope for this species (deg C yr^{-1} or m yr^{-1}).

individual_slope_se

Numeric.Standard error of the species-specific slope.

individual_t

Numeric.t-statistic for the species slope (H0: slope = 0).

individual_p_value

Numeric.Two-tailed p-value for the species slope test (raw, before FDR correction).

individual_conf_low

Numeric.Lower bound of the 95 \ confidence interval for the species slope.

individual_conf_high

Numeric.Upper bound of the 95 \ confidence interval for the species slope.

global_lm_slope

Numeric.Global slope for this response variable (copied from environmental_global_lm for convenience).

slope_diff_signed

Numeric.Signed difference between the species slope and the global slope (individual_slope - global_lm_slope).Positive values indicate the species trend is steeper than the community average.

slope_diff_direction

Character.Direction of slope_diff_signed: "positive", "negative", or "zero".

slope_diff

Numeric.Absolute value of slope_diff_signed.

slope_diff_t

Numeric.Welch-Satterthwaite t-statistic for the test of equality between the species slope and the global slope.

slope_diff_df

Numeric.Welch-Satterthwaite degrees of freedom for slope_diff_t.

slope_diff_p_value

Numeric.Two-tailed p-value for the slope-difference test (raw, before FDR correction).

individual_direction

Character.Direction of the species slope itself: "positive", "negative", or "zero".

individual_p_adj_fdr

Numeric.FDR-adjusted (BH method) p-value for the species slope test, computed within each response variable separately.

slope_diff_p_adj_fdr

Numeric.FDR-adjusted (BH method) p-value for the slope-difference test.

alpha

Numeric.Effective significance threshold applied to this species.Equals alpha_ref for species with <= 100 records, or the sample-size-corrected value when use_Ncorrected_alpha = TRUE and N > 100.

individual_significant

Logical.TRUE if individual_p_adj_fdr < alpha for this species.

slope_diff_significant

Logical.TRUE if slope_diff_p_adj_fdr < alpha for this species.

environmental_comparison

Data frame extending environmental_individual_lm with the classification columns added by classify_lm_individual_slope().Contains all the columns described for environmental_individual_lm plus:

individual_slope_above_threshold

Logical.TRUE if the relevant quantity (species slope or slope difference, depending on individual_lm_threshold_type) exceeds the minimum ecological threshold defined in individual_lm_min_thresholds. Always TRUE when no threshold is specified.Retained for downstream quality audits.

individual_response_class

Character.Final environmental response classification for this species and response variable. Possible values: "TT" (Thermal Tolerance), "TA" (Thermal Adjustment), "TC" (Thermal Conformance) for Temperature; "SA" (Spatial Adaptation), "SD" (Spatial Discordance), "SC" (Spatial Conformance) for Elevation.A non-conformance class (TT, TA, SA, SD) is assigned only when both individual_significant and slope_diff_significant are TRUE and individual_slope_above_threshold is TRUE. Otherwise the species defaults to the conformance class (TC or SC).

Classification system

Each species is assigned an environmental response class reflecting the direction and significance of its temporal trend (or its deviation from the global trend, depending on individual_lm_threshold_type):

"TT" - Thermal Tolerance

Temperature response.The species is moving toward warmer environments over time: its individual slope (or its deviation from the global slope, under "slopediff") is positive, statistically significant, and exceeds the minimum ecological threshold.

"TA" - Thermal Adjustment

Temperature response.The species is shifting toward cooler environments: slope (or deviation) is negative, significant, and exceeds the threshold.

"TC" - Thermal Conformance

Temperature response.Default class assigned when the individual trend is non-significant or not distinguishable from the global trend.

"SA" - Spatial Adaptation

Elevation response.The species is tracking higher elevations over time: slope (or deviation) is positive, significant, and exceeds the threshold.

"SD" - Spatial Discordance

Elevation response.The species is shifting to lower elevations: slope (or deviation) is negative, significant, and exceeds the threshold.

"SC" - Spatial Conformance

Elevation response.Default class when the elevation trend is non-significant or indistinguishable from the global trend.

Time variable

Internally the function constructs a continuous, mean-centred time predictor (time_cont_c) from Year and Month as Year + (Month - 1) / 12 - mean(...). Slopes are therefore expressed in environmental units per year (deg C yr^{-1} for Temperature, m yr^{-1} for Elevation).

Examples

# Minimal example with simulated data
set.seed(42)
n <- 600
sim_data <- data.frame(
Species = sample(paste0("Sp", 1:5), n, replace = TRUE),
Year = sample(2000:2020, n, replace = TRUE),
Month = sample(1:12, n, replace = TRUE),
Temperature = rnorm(n, mean = 15, sd = 3),
Elevation = rnorm(n, mean = 500, sd = 100))

results <- spp_trend_environmental(
data = sim_data,
responses = c("Temperature", "Elevation"),
min_records = 20,
min_years = 5,
alpha_ref = 0.05,
use_Ncorrected_alpha = FALSE,
individual_lm_min_thresholds = c(Temperature = 0.01, Elevation = 50),
individual_lm_threshold_type = "speciesslope",
verbose = FALSE
)

# Inspect the classification table
head(results$environmental_comparison)

# Species retained per response
results$species_filter[results$species_filter$retained, ]

# Global community trends
results$environmental_global_lm

Analyse Temporal Changes in Species Geographical Position Using ECEF Vector Analysis

Description

Estimates temporal trends in the spatial distribution of species by modelling the annual displacement of each species' geographical centroid in three-dimensional Earth-Centred Earth-Fixed (ECEF) Cartesian coordinates derived from the WGS84 reference ellipsoid.

For each species (and for a pooled global reference), the function computes monthly ECEF centroids from occurrence records, fits a multivariate linear regression of centroid position against a continuous time variable (\mathrm{time\_cont} = \mathrm{Year} + (\mathrm{Month} - 1) / 12), and projects the resulting three-dimensional slope vector onto the local tangent plane at the species centroid to obtain ecologically interpretable quantities: surface displacement speed (km yr^{-1}) and absolute geographical bearing.

Each species is compared against a single global reference vector estimated from all retained records pooled together. Uncertainty is propagated via covariance-based multivariate simulation, and each species is assigned a spatial class:

SA (Spatial Adaptation)

The species moves faster than the global pool, or moves in a significantly different direction that is poleward, consistent with a climate-driven distributional shift.

SD (Spatial Discordance)

The species moves slower than the global pool, or moves in a significantly different direction that is non-poleward, inconsistent with strict climatic control of distribution.

SC (Spatial Conformance)

Neither speed nor direction differs significantly from the global spatial pattern.

The ECEF coordinate system used here places the X-axis towards the intersection of the Equator and the Greenwich meridian, the Y-axis towards 90 deg E longitude, and the Z-axis towards the North Pole.

Usage

spp_trend_spatial(
  data,
  spp = NULL,
  min_records = 20,
  min_years = 5,
  spatial_simulation_n = 1000,
  spatial_probability_threshold = 0.9,
  direction_angle_threshold_deg = 68,
  random_seed = NULL,
  verbose = TRUE
)

Arguments

Details

Coordinate transformation

Latitude and longitude (decimal degrees, WGS84) are converted to ECEF Cartesian coordinates (kilometres) using the standard geodetic transformation with WGS84 parameters: semi-major axis a = 6378.137 km and flattening f = 1/298.257223563. The prime vertical radius of curvature N(\phi) = a / \sqrt{1 - e^2 \sin^2\phi} accounts for the ellipsoidal shape of the Earth.

Temporal centroids

Records are grouped by species and year-month combination. Within each group the arithmetic mean of the ECEF coordinates is computed to produce a temporal centroid. Groups with fewer records than \lfloor \texttt{min\_records} / \max(1, \texttt{min\_years}) \rfloor are discarded. Global-pool centroids use all retained species combined.

ECEF slope model

A multivariate ordinary least-squares regression is fitted:

\mathbf{C}(t) = \mathbf{\alpha} + \mathbf{\beta}\, t + \mathbf{\varepsilon}

where \mathbf{C}(t) is the 3 \times 1 vector of ECEF coordinates at time t, and \mathbf{\beta} is the slope vector (km yr^{-1}) of interest. The residual covariance matrix \hat{\mathbf{\Sigma}} is used together with (X^\top X)^{-1} to construct the 3 \times 3 covariance matrix of \hat{\mathbf{\beta}}, which propagates uncertainty into the Monte Carlo simulations.

Tangent-plane projection

The ECEF slope vector is projected onto the local east-north-up frame at the species centroid. The east and north components are used to compute surface speed and bearing; the up component is reported as radial_km_year.

Spatial classification hierarchy

1. Direction test. If prob_direction_different_from_global >= spatial_probability_threshold and angle_to_global_3d_deg >= direction_angle_threshold_deg, the direction is considered significantly different. Under this condition, poleward movement -> SA; non-poleward movement -> SD.

2. Speed test (direction not significant). If prob_speed_greater_than_global >= spatial_probability_threshold, poleward movement -> SA; non-poleward -> SD. If prob_speed_lower_than_global >= spatial_probability_threshold -> SD.

3. Default. Neither speed nor direction differs significantly -> SC.

Poleward movement is defined as a positive northward component in the Northern Hemisphere or a negative northward component in the Southern Hemisphere.

Value

A named list with three elements:

spatial

A data.frame with one row per entity analysed. The first row always corresponds to the global pooled reference (Species == "__GLOBAL_POOL__"); subsequent rows correspond to individual species. Columns:

Species

Character. Species name, or "__GLOBAL_POOL__" for the global reference row.

analysis

Character. Identifies the type of analysis: "spatial_ecef_centroid_vector_global_pool" for the global row or "spatial_ecef_centroid_vector_species" for individual species.

n_records

Integer. Total number of valid occurrence records used. For the global pool this is the sum across all temporal centroids.

n_time_steps

Integer. Number of year-month combinations (temporal centroids) used to fit the ECEF slope model.

centroid_Longitude

Numeric. Decimal degrees longitude of the mean ECEF centroid, converted back to geodetic coordinates.

centroid_Latitude

Numeric. Decimal degrees latitude of the mean ECEF centroid, converted back to geodetic coordinates.

hemisphere

Character. "Northern" if centroid_Latitude >= 0, otherwise "Southern".

poleward_status

Character. "poleward" when the northward component of the displacement vector is positive in the Northern Hemisphere or negative in the Southern Hemisphere; "non_poleward" otherwise. NA for the global pool row.

slope_ecef_x_km_year

Numeric. Annual rate of change (km yr^{-1}) of the ECEF X-coordinate centroid, estimated from the multivariate linear regression.

slope_ecef_y_km_year

Numeric. Annual rate of change (km yr^{-1}) of the ECEF Y-coordinate centroid.

slope_ecef_z_km_year

Numeric. Annual rate of change (km yr^{-1}) of the ECEF Z-coordinate centroid.

speed_chord_km_year

Numeric. Magnitude of the 3-D ECEF slope vector (km yr^{-1}), representing the straight-line chord displacement speed through the Earth's interior. Equal to \sqrt{\dot{X}^2 + \dot{Y}^2 + \dot{Z}^2}.

speed_surface_km_year

Numeric. Horizontal surface displacement speed (km yr^{-1}) obtained by projecting the ECEF slope vector onto the local tangent plane at the centroid and taking the magnitude of the east and north components. This is the ecologically preferred speed metric.

direction_bearing_deg

Numeric. Absolute geographical bearing (0-360 deg, clockwise from North) of the displacement vector on the local tangent plane. NA when the surface speed is zero or not finite.

direction_cardinal

Character. Eight-point compass direction corresponding to direction_bearing_deg: one of "N", "NE", "E", "SE", "S", "SW", "W", or "NW".

east_km_year

Numeric. Eastward component (km yr^{-1}) of the tangent-plane displacement vector. Positive values indicate eastward movement.

north_km_year

Numeric. Northward component (km yr^{-1}) of the tangent-plane displacement vector. Positive values indicate northward movement in geographic coordinates.

radial_km_year

Numeric. Radial (vertical) component (km yr^{-1}) of the ECEF slope vector projected onto the local up direction. Non-zero values indicate apparent elevation change and should be close to zero for purely horizontal distributional shifts.

global_speed_surface_km_year_full

Numeric. Surface displacement speed (km yr^{-1}) of the global reference vector (estimated from all retained records). Provided as a fixed reference for comparison. NA for the global pool row itself.

speed_diff_surface_km_year

Numeric. Difference between the species surface speed and the global surface speed (km yr^{-1}). Positive values indicate the species is moving faster than the global pattern. NA for the global pool row.

speed_ratio_to_global

Numeric. Ratio of the species surface speed to the global surface speed. A value greater than 1 indicates faster movement than the global reference. NA when the global speed is zero or for the global pool row.

angle_to_global_3d_deg

Numeric. Three-dimensional angle (degrees, 0-180) between the species ECEF slope vector and the global ECEF slope vector. Values near 0 deg indicate nearly parallel movement; values near 90 deg indicate orthogonal movement; values near 180 deg indicate opposite movement. NA for the global pool row.

prob_speed_greater_than_global

Numeric (0-1). Estimated probability, derived from spatial_simulation_n Monte Carlo draws, that the species surface speed exceeds the global surface speed. NA for the global pool row.

prob_speed_lower_than_global

Numeric (0-1). Estimated probability that the species surface speed is lower than the global surface speed. NA for the global pool row.

prob_direction_different_from_global

Numeric (0-1). Estimated probability that the three-dimensional angle between the species slope vector and the global slope vector is at least direction_angle_threshold_deg degrees, indicating directional divergence from the global pattern. NA for the global pool row.

direction_angle_threshold_deg

Numeric. The value of the direction_angle_threshold_deg argument used in the analysis, stored here for traceability.

spatial_class_Latitudelon

Character. Joint spatial classification assigned to each species: "SA" (Spatial Adaptation), "SD" (Spatial Discordance), or "SC" (Spatial Conformance). See Details for the full decision hierarchy. NA for the global pool row.

species_filter

A data.frame with one row per species found in data (after optional spp subsetting). Records whether each species passed the data-sufficiency filters. Columns:

Species

Character. Species name.

n_records

Integer. Total number of valid records for this species after cleaning.

n_years

Integer. Number of distinct calendar years in which the species was recorded.

retained

Logical. TRUE if the species met both min_records and min_years thresholds and was included in the spatial analysis; FALSE otherwise.

metadata

A named list documenting analysis settings and methodological details. Elements:

analysis

Character. Fixed label "spatial_ecef_wgs84_vector_only".

time_variable

Character. Description of the continuous time variable: "time_cont = Year + (Month - 1) / 12".

coordinate_system

Character. Coordinate system and units: "WGS84 ECEF, kilometres".

model

Character. Formula of the multivariate regression: "cbind(ecef_x_km, ecef_y_km, ecef_z_km) ~ time_cont".

joint_spatial_classification

Character. Human-readable description of the hierarchical decision rule used to assign SA, SD, or SC, including the threshold values applied.

spatial_simulation_n

Integer. Value of the spatial_simulation_n argument used.

spatial_probability_threshold

Numeric. Value of the spatial_probability_threshold argument used.

direction_angle_threshold_deg

Numeric. Value of the direction_angle_threshold_deg argument used.

min_records

Integer. Value of the min_records argument used.

min_years

Integer. Value of the min_years argument used.

See Also

lm for the underlying regression engine.

Examples

# Minimal example with simulated data
set.seed(42)
n <- 600
sim_data <- data.frame(Species = sample(c("Sp_A", "Sp_B", "Sp_C"), n, replace = TRUE),
Year = sample(2000:2020, n, replace = TRUE),
Month = sample(1:12, n, replace = TRUE),
Latitude = runif(n, 30, 65),
Longitude = runif(n, -10, 30)
)

result <- spp_trend_spatial(
data= sim_data,
min_records = 20,
min_years = 5,
spatial_simulation_n = 100,
spatial_probability_threshold = 0.90,
direction_angle_threshold_deg = 68,
random_seed = 1,
verbose = FALSE
)

# Main results table
head(result$spatial)

# Species filtering summary
result$species_filter

# Analysis metadata
result$metadata

# Extract classified species only (excluding global pool row)
classified <- result$spatial[result$spatial$Species != "__GLOBAL_POOL__", ]
table(classified$spatial_class_Latitudelon)