Movement phases shape wolf habitat selection in human-dominated Central European landscapes

Code Repository: Movement phases shape wolf habitat selection in human-dominated Central European landscapes Overview This repository contains the R and shell scripts used to run the integrated Step Selection Analysis (iSSA) and the subsequent post-hoc evaluations of animal movement distributions presented in the manuscript, together with the two fitted model objects. The code documents the pipeline from model fitting (conditional logistic regression approximated via a Poisson formulation with stratum-specific intercepts, following Muff et al. 2020, J. Anim. Ecol. 89, 80-92) to the extraction of updated movement parameters. All three movement phases are estimated within a single joint model, so that every comparison between phases is obtained as a linear contrast within one variance-covariance matrix of the fixed effects, rather than by comparing coefficients across separately fitted models. The upstream preparation of the input dataset from raw GPS fixes and environmental layers is described in the Methods and is not included here. Data Availability The GPS tracking data and the annotated input dataset derived from them (new_dat_issa_n50.RData, observed steps with 50 generated random steps each) are not included. They contain precise locations of individuals of a strictly protected species, and their public release would expose the animals to a risk of persecution. They are withheld for conservation reasons, as explained in the Data accessibility statement of the manuscript. The fitted model objects contain no locations and are included: n50_iSSA_general_idonly.RData (about 14 GB): output of script 01, the general model. n50_iSSA_inter_idonly.RData (about 15 GB): output of script 02, the joint phase model underlying all comparisons between phases reported in the manuscript. Each file holds a fitted glmmTMB object named model_result. Loading them requires several hundred gigabytes of memory. Scripts 03 and 04 run on these files alone. Script 05 additionally needs the input dataset to fit the tentative movement distributions, and script 06 needs the output of script 05. Software R 4.3.3, glmmTMB 1.1.9, TMB 1.9.17, multcomp 1.4.28, circular 0.5.1, tidyverse 2.0.0. The models were fitted on an HPC cluster with 8 CPUs and up to 600 GB of memory. Directory Structure and File Descriptions The scripts are numbered sequentially to represent the analytical pipeline. Each .R script has a matching .sh batch script (PBS) for submitting the job to a high-performance computing (HPC) cluster, where all steps were run. The batch scripts expect the environment variables ISSA_PROJECT and ISSA_CONDA_ENV to point to the project directory and to the conda environment holding R. 01_model_general.R: Fits the general iSSA to the annotated GPS data (observed steps and 50 random steps per observed step), with no term interacted with the movement phase, so the estimates describe habitat selection generalised across individuals. Same data and random-effect structure as the joint model (script 02). Writes results/n50_iSSA_general_idonly.RData. 02_model_phase.R: The core modelling script. It fits the joint iSSA with 66 fixed effects and 18 uncorrelated random slopes specified at the individual level, (0 + covariate | id), one for each of the nine environmental covariates and one for each of their nine interactions with diel period, giving 18 estimated variance parameters. The random-effect structure is shared across movement phases. The variance of the stratum-specific intercepts is fixed at a large value, as in Muff et al. 2020, to implement the Poisson approximation. The model uses a cell-means parametrisation, in which every covariate is interacted with the movement phase, so that each coefficient is a direct estimate for the given phase rather than a difference from a reference phase. Writes results/n50_iSSA_inter_idonly.RData. 03_contrasts_general.R: Derives day, night, and night-against-day estimates from the general model produced by script 01, using linear hypothesis tests (multcomp package) with Benjamini-Hochberg correction applied within scenario. Writes results/general_daynight_idonly.csv. 04_contrasts_phase.R: Derives all comparisons from the joint model produced by script 02 as linear contrasts, in three families, each corrected separately with the Benjamini-Hochberg method: effects within each phase (day, night and night-against-day), differences between phases by day, and differences between phases by night. Writes results/posthoc_within_strategy_idonly.csv and results/posthoc_between_strategy_idonly.csv. 05_movement_distributions.R: Updates the tentative step-length (gamma) and turn-angle (von Mises) distributions, i.e. those fitted to observed steps before the habitat-selection update, using the fitted joint model coefficients, separately for each phase and diel period, and exports the parameters required to draw the movement figures. Writes distr_tentative_idonly.RData, distr_updated_idonly.RData and results/movement_plot_data.RData. 06_displacement_turnangle.R: Calculates the expected displacement per movement step, using the delta method to propagate variance from both the tentative distributions and the model, and tests displacement and turn-angle concentration between phases. Requires distr_tentative_idonly.RData, written by script 05. Writes results/displacement_idonly_results.RData and tables in Tables/. Notes Phase labels in file and column names use the term strategy for historical reasons; in the manuscript these are referred to as movement phases. Turn-angle concentration is reported in the manuscript as the absolute value |κ|. A negative concentration parameter denotes concentration around ±π (reversals).