Two decades of community responses reveal a consistent increase in species occurrences in south-eastern Pacific rocky shores.

README

Description

########################################################################################
This repository contains seven R scripts for setting up, fitting and post-processing HMSC analyses of intertidal species. 

Ensure that your “working directory” (setwd) includes the following folders: data, models, and results.

The scripts must be executed in the following order: S0, S1, S2, S3, S4, S5, and S6.
########################################################################################

S0: Data input
The S0 script defines the sampling effort and timing of each survey. It includes the macroalgae and invertebrates abundance datasets, species traits, taxonomy and sea surface temperature (SST) data.

Species abundance (and METADATA)
Running S0 will load two data frames: da1_algae (macroalgae) and da1_invr (invertebrates). Both have the same number of rows and columns, with following variables in columns 1-9:

SITE: All sampled sites (22 levels).
HEIGHT: Intertidal zone (3 levels: HIGH, MID and LOW).
N_CUA: Number of quadrats per survey (50 x 50 cm).
MONTH: Month of survey.
YEAR: Year of survey.
BENCH: Platform at each site (2 levels: 1 and 2).
SOUTHLAT: South Latitude of site.
WESTLONG: West Longitude of site.
SEASON: Season of survey.

The remaining columns contain the species recorded at each survey, depending on the data frame: da1_algae (macroalgae) or da1_invr (invertebrates). 

In da1_algae, species abundances are expressed as percentage cover (%). 

In da1_invr, species abundances are recorded according to mobility described in the species traits data frame (see below):

• Sessile species: percentage cover (%).
• Mobile species: number of individuals per quadrat.

########################################################################################

Species Traits and Taxonomy (and METADATA)
Species traits and taxonomy are defined in the data frames algae_tr (macroalgae) and invr_tr (invertebrates). Both data frames share the taxonomic classification but include different sets of traits.

• algae_tr (macroalgae traits): 

LIFE_CYCLE_MOR: Life cycle (3 levels).
LIFESPAN: Life span (2 levels).
MAX_LENGTH_cm: Maximum length that the species would reach (cm).

• invr_tr (invertebrates traits): 

TIERING: Tiering (3 levels).
FEEDING: Feeding mechanisms (4 levels).
LARVAL_DEV: Larval development (3 levels).
MAX_SIZE_cm: Maximum length that the species would reach (cm).

Species taxonomy is provided in the objects my.taxonomy_algae (macroalgae) and my.taxonomy_invr (invertebrates).

########################################################################################

Sea Surface Temperature (SST) data (and METADATA)
The data frame da2 contains five environmental predictors with following variables: 

SITE: All sampled sites (22 levels).
LAT: South Latitude of site.
LONG: West Longitude of site.
YEAR: Year of survey.
ANN_ANOM: Standardized annual SST anomaly per site.
ANN_MEAN: Annual mean SST per site.
qFifth: 5^th quintile of SST per site.
qFifty: 50^th quintile (median) of SST per site.
qNinety5: 95^th quintile of SST per site.

########################################################################################

XDATA (and METADATA)
The essential data for HMSC analyses include a XDATA object, which is defined by five variables: 

year: Year of survey.
height: Intertidal zone (3 levels: HIGH, MID and LOW).
season: Season of survey.
ann_anom: Standardized annual SST anomaly per site.
ann_mean: Annual mean of SST per site.

########################################################################################

Setting up for HMSC models
By saving the objects described above in a single data frame (called “raw_data_hmsc”) in your data/ folder, you obtain the inputs required for the HMSC analyses:

1) Y_algae and Y_invr: Matrices of intertidal species abundances (macroalgae and invertebrates, respectively).
2) X_Data: Matrix of environmental predictors.
3) TrData_algae and TrData_invr: Macroalgal and invertebrate traits.
4) my.taxonomy_algae and my.taxonomy_invr: Phylogenetic-taxonomic relationships.
5) studyDesign: Spatiotemporal sampling design of the study.

###################################END#################################################

S1: Define models
To define the HMSC models, start by loading the raw_data_hmsc object:
load("data/raw_data_hmsc.RData")

This will load eight data frames:

1)    First and second: "my.taxonomy_algae" and "my.taxonomy_invr".
Taxonomic classification tables (6 columns): kingdom, phylum, family, genus, species.

2)    studyDesign.
Spatiotemporal sampling design (3 columns): site, bench, year.

3)   TrData_algae and TrData_invr.
Species traits tables:

• TrData_algae (3 columns): life_cycle, max_length", lifespan".
• TrData_invr (4 columns): tiering, feeding, larva, max_size.

 4)    XData.
Environmental predictors (4 columns): year, height, ann_anom, ann_mean.

 5)    Y_algae and Y_invr.
Species abundances matrices:

• Y_algae: 54 macroalgal species.
• Y_invr: 57 invertebrate species.

########################################################################################

Saving model definitions: Once the objects have been loaded and structured, define and save the initial HMSC models for both presence-absence and abundance of macroalgae and invertebrates. Stored these model objects in the models/ folder so they are available for the subsequent steps (S2-S6), which involve fitting, evaluating, and postprocessing.

########################################################################################

End of S1: At this stage, you have prepared the data, defined the model structures, and saved the base models. Proceed to S2, where the models are fitted.

###################################END#################################################

S2: Check MCMC convergence
Evaluates MCMC convergence of the fitted HMSC models.

Input: fitted models objects.
Output: convergence diagnostics and plots (results/MCMC_convergence.txt and results/MCMC_convergence.pdf).

###################################END#################################################

S3: Evaluate model fit (CV + WAIC)
Computes in-sample fit (MF), cross-validated fit (MFCV) via n-fold CV, and WAIC for the fitted HMSC models across increasing MCMC lengths (matching S2 runs).

Input: fitted models objects (models_thin_*_samples_*_chains_4.Rdata).
Output: MF, MFCV, and WAIC saved to models/MF_thin_*_samples_*_chains_4_nfolds_*.Rdata.
Notes: defaults to 2-fold CV; supports parallel evaluation via nParallel. 

###################################END#################################################

S4: Visualize model fit (explanatory vs. predictive)
Loads the latest model-fit objects (MF/MFCV) and plots in-sample vs. cross-validated performance for each model: Tjur R², R², AUC, and SR².

Input: models/MF_thin_*_samples_*_chains_4_nfolds_*.Rdata (highest available run auto-selected).
Output: results/model_fit_nfolds_<n>.pdf with scatterplots of explanatory vs. predictive power.
Note: defaults to 2-fold CV unless n-folds is changed.

###################################END#################################################

S5: Parameter estimates & variance partitioning
Generates summaries and plots of parameter estimates and variance partitioning, including trait effects and species associations.

Input: fitted models (models_thin_*_samples_*_chains_4.Rdata)
Output:

• parameter_estimates.pdf (variance partitioning, Beta, Gamma, Omega plots).
• parameter_estimates.txt (species order/orderings).
• CSV/XLSX files with numerical estimates (Beta, Gamma, Omega, VP).

###################################END#################################################

S6: Site-specific community responses (CWMR)
Computes community-weighted mean responses (CWMR) per site for each covariate using posterior Beta estimates, for both occurrence and abundance; also summarizes fractions of positively/negatively supported effect (≥0.95).

Input: models/models_thin_100_samples_250_chains_4.Rdata.
Output: site_specific_responses_<model>.RData containing:           

• CWMR.occurrence, CWMR.abundance.
• pos.occurrence, pos.abundance, neg.occurrence, neg.abundance.

###################################END#################################################