Data and R code from: Spatiotemporal risk factors predict landscape-scale survivorship for a northern ungulate

These data and computer code (written in R, https://www.r-project.org) were created to statistically evaluate a suite of spatiotemporal covariates that could potentially explain pronghorn (Antilocapra americana) mortality risk in the Northern Sagebrush Steppe (NSS) ecosystem (50.0757o N, −108.7526o W). Known-fate data were collected from 170 adult female pronghorn monitored with GPS collars from 2003-2011, which were used to construct a time-to-event (TTE) dataset with a daily timescale and an annual recurrent origin of 11 November. Seasonal risk periods (winter, spring, summer, autumn) were defined by median migration dates of collared pronghorn. We linked this TTE dataset with spatiotemporal covariates that were extracted and collated from pronghorn seasonal activity areas (estimated using 95% minimum convex polygons) to form a final dataset. Specifically, average fence and road densities (km/km2), average snow water equivalent (SWE; kg/m2), and maximum decadal normalized difference vegetation index (NDVI) were considered as predictors. We tested for these main effects of spatiotemporal risk covariates as well as the hypotheses that pronghorn mortality risk from roads or fences could be intensified during severe winter weather (i.e., interactions: SWE*road density and SWE*fence density). We also compare an analogous frequentist implementation to estimate model-averaged risk coefficients. Ultimately, the study aimed to develop the first broad-scale, spatially explicit map of predicted annual pronghorn survivorship based on anthropogenic features and environmental gradients to identify areas for conservation and habitat restoration efforts.

Methods
We combined relocations from GPS-collared adult female pronghorn (n = 170) with raster data that described potentially important spatiotemporal risk covariates. We first collated relocation and time-to-event data to remove individual pronghorn from the analysis that had no spatial data available. We then constructed seasonal risk periods based on the median migration dates determined from a previous analysis; thus, we defined 4 seasonal periods as winter (11 November–21 March), spring (22 March–10 April), summer (11 April–30 October), and autumn (31 October–10 November). We used the package 'amt' in Program R to rarify relocation data to a common 4-hr interval using a 30-min tolerance. We used the package 'adehabitatHR' in Program R to estimate seasonal activity areas using 95% minimum convex polygon. We constructed annual- and seasonal-specific risk covariates by averaging values within individual activity areas. We specifically extracted values for linear features (road and fence densities), a proxy for snow depth (SWE), and a measure of forage productivity (NDVI). We resampled all raster data to a common resolution of 1 km2. Given that fence density models characterized regional-scale variation in fence density (i.e., 1.5 km2), this resolution seemed appropriate for our risk analysis. 

We fit Bayesian proportional hazards (PH) models using a time-to-event approach to model the effects of spatiotemporal covariates on pronghorn mortality risk. We aimed to develop a model to understand the relative effects of risk covariates for pronghorn in the NSS. The effect of fence or road densities may depend on SWE such that the variables interact in affecting mortality risk. Thus, our full candidate model included four main effects and two interaction terms. We used reversible-jump Markov Chain Monte Carlo (RJMCMC) to determine relative support for a nested set of Bayesian PH models. This allowed us to conduct Bayesian model selection and averaging in one step by using two custom samplers provided for the R package 'nimble'. 

For brevity, we provide the final time-to-event dataset and analysis code rather than include all of the code, GIS, etc. used to estimate seasonal activity areas and extract and collate spatial risk covariates for each individual. Rather we provide the data and all code to reproduce the risk regression results presented in the manuscript.

本数据集与配套计算机代码（基于R语言编写，https://www.r-project.org）旨在通过统计方法评估一系列时空协变量，以解释北灌丛草原（Northern Sagebrush Steppe, NSS）生态系统（坐标：50.0757°N，−108.7526°W）内叉角羚（Antilocapra americana）的死亡风险。研究于2003-2011年间收集了170只佩戴GPS项圈的成年雌性叉角羚的已知命运数据，以此构建了以每日为时间步长、以每年11月11日为重复起始点的生存时间（Time-to-Event, TTE）数据集。基于佩戴项圈的叉角羚的中位迁徙日期，划定了四季风险期：冬季、春季、夏季与秋季。本研究将该TTE数据集与从叉角羚季节性活动区（通过95%最小凸多边形估算得到）提取并整理得到的时空协变量进行关联，最终形成完整数据集。具体而言，本研究选取平均围栏与道路密度（km/km²）、平均雪水当量（Snow Water Equivalent, SWE；单位：kg/m²）以及十年最大归一化植被指数（Normalized Difference Vegetation Index, NDVI）作为预测变量。我们检验了这些时空风险协变量的主效应，同时验证了以下假说：在极端冬季天气条件下，道路或围栏对叉角羚死亡风险的影响会被加剧（即交互项：SWE×道路密度与SWE×围栏密度）。此外，我们还对比了用于估算模型平均风险系数的同类频率学派方法。本研究最终旨在基于人为活动特征与环境梯度，绘制首个大范围空间显式的叉角羚年度存活率预测地图，以识别适合开展保护与栖息地修复工作的区域。

方法
我们将佩戴GPS项圈的成年雌性叉角羚的定位点数据（n=170）与描述潜在重要时空风险协变量的栅格数据进行整合。首先，我们整理定位点与生存时间数据，剔除无有效空间数据的个体。随后，基于前期分析得到的中位迁徙日期划定季节性风险期：将冬季定义为11月11日至次年3月21日，春季为3月22日至4月10日，夏季为4月11日至10月30日，秋季为10月31日至11月10日。我们使用R语言中的‘amt’包，以30分钟为容差，将定位数据稀疏化为统一的4小时间隔。随后通过R语言‘adehabitatHR’包，基于95%最小凸多边形估算季节性活动区。我们通过对个体活动区内的协变量值取平均，构建年度与季节特异性的风险协变量。具体提取的变量包括线性地表特征（道路与围栏密度）、雪深替代指标SWE，以及采食生产力指标NDVI。我们将所有栅格数据重采样至统一的1 km²分辨率。鉴于围栏密度模型刻画了区域尺度的围栏密度变异（即1.5 km²），该分辨率适配本研究的风险分析需求。

我们采用生存时间分析方法拟合贝叶斯比例风险（Proportional Hazards, PH）模型，以探究时空协变量对叉角羚死亡风险的影响。本研究旨在构建模型以理解NSS生态系统中叉角羚风险协变量的相对效应。围栏或道路密度的效应可能依赖于SWE，即二者存在交互作用共同影响死亡风险。因此，我们的完整候选模型包含4个主效应项与2个交互项。我们使用可逆跳跃马尔可夫链蒙特卡洛（Reversible-jump Markov Chain Monte Carlo, RJMCMC）方法，对一系列嵌套的贝叶斯PH模型进行相对支持度评估。通过R包‘nimble’提供的两种自定义采样器，我们可以一步完成贝叶斯模型选择与模型平均。

为简洁起见，本研究仅提供最终的生存时间数据集与分析代码，而非用于估算季节性活动区、提取并整理个体空间风险协变量的全部代码、GIS工具等文件。我们仅提供可复现论文中风险回归结果的数据与完整代码。