《Ecological Modelling》：A Censored Seasonal Trend Model for Marine Environmental Time Series with R-INLA
编辑推荐：
摘要在全球变化背景下，分析海洋环境时间序列对于识别趋势、季节性变化以及超出临界阈值的情形至关重要，从而为决策提供支持。然而，这类时间序列往往存在量化限制（即截断数据）、不规律的采样、缺失值以及多模态季节性特征。若用替代值或忽略截断值，将会在后续分析中引入显著偏差。本研究提出了截断
摘要
在全球变化背景下，分析海洋环境时间序列对于识别趋势、季节性变化以及超出临界阈值的情形至关重要，从而为决策提供支持。然而，这类时间序列往往存在量化限制（即截断数据）、不规律的采样、缺失值以及多模态季节性特征。若用替代值或忽略截断值，将会在后续分析中引入显著偏差。本研究提出了截断趋势季节性模型（CSTM），这是一种专门用于分析截断环境时间序列的结构模型。该模型结合二阶随机游走模型来表示趋势，同时运用三角函数来分解季节性成分，进而将时间序列分为趋势、季节性和噪声三部分。模型假设观测值服从对数正态分布，并利用R-INLA软件包进行高效的贝叶斯推断，无需特殊插补即可直接处理左截断数据。模型拟合效果通过随机分位数残差法进行严格评估，该方法可对结果进行标准化处理，以便使用传统的拟合优度指标。该模型在模拟数据以及法国IFREMER的REMI监测网络所收集的贝类中大肠杆菌浓度这一真实数据集上均得到了验证。实验结果表明，该模型能够准确估算浓度——即便是在量化限制之下——同时还能识别出趋势和季节性模式。CSTM的灵活性、计算效率及稳定性使其非常适合大规模、自动化的环境时间序列分析，可为生态研究、监管监测以及公共卫生决策提供有力支持。此外，该模型还为分析其他具有截断特征的环境变量带来了新的可能性。
引言
在全球变暖的背景下，监测海洋环境变量再次受到重视（Koslow和Couture，2013）。在法国，法国海洋开发研究院（IFREMER）数十年来一直运营着国家和区域的实地监测网络，产生了数以万计的持续时间序列数据。这些数据通常用于识别趋势、季节性变化、物候现象或临界阈值被突破的情况，从而为决策提供依据。然而，不规律的采样、缺失数据以及多模态季节性特征给数据处理带来了困难。结构分析方法（Harvey，1990；West和Harrison，1997）可以有效解决这些问题（Ratmaya等人，2019）。另一个难题是，测量过程常常会给出诸如“小于值”这样的结果，这类数值代表了量化限制值或检测极限。1在海洋环境中，检测贝类中的大肠杆菌、沉积物与贝类中的污染物，或是海水中的叶绿素a和营养物质，需要采用不同的方法、设备和操作人员，因此量化限制值也各不相同。这类结果甚至可能是隐性的。例如，如果专家在海水样本中统计微浮游植物种类，而某种专家非常熟悉的常见物种并未出现，这并不意味着浓度为零，而是说明其浓度低于量化限制值。如果在10毫升的沉淀室中进行计数，那么量化限制值就是100个细胞/升。此外，随着时间推移，方法、实验室、设备以及操作人员都会发生变化，从而导致量化限制值改变，进而形成具有多个量化限制值的时间序列。截断数据的统计难题早已为人所知。忽略或替换截断测量值会引入偏差，影响对位置、趋势、季节性、变异程度以及统计检验结果的准确判断（Helsel，2006）。在时间序列分析中，关键在于如何恰当地处理无法检测到的数据，从而客观地分析季节性调整后的趋势等结构成分。
自20世纪90年代以来，Dennis Helsel在水质领域做出了重要贡献，他不仅强调了忽略这些测量数据或忽视其特殊性质的潜在风险，还提出了有效管理这些数据的方法。为了估算截断时间序列的趋势，Helsel等人（2020）推荐使用非参数方法（如Mann–Kendall、Akritas–Theil–Sen）以及基于最大似然估计的Tobit回归方法来检测单调趋势。尽管关于截断回归的文献较多，但正如Wang和Chan（2017）所指出的，针对截断时间序列的研究相对较少。他们提出了在截断变量与具有自回归误差的协变量之间建立线性回归关系，但这需要存在相关的协变量。Mohammad（2014）开发了一种准EM算法，用于将ARMA、平稳ARFIMA以及其他线性模型应用于截断时间序列。Durbin和Koopman（2012）指出，与ARIMA相比，状态空间模型在分解趋势、季节性和噪声等结构成分方面更具优势。Allik等人（2015）将卡尔曼滤波器应用于截断数据，但并未发展出具备平滑功能的反向滤波器。Courbariaux等人（2022）提出了一种卡尔曼平滑器，可用于处理检测阈值（即数据的左截断）和异常值，但由于采用了状态空间离散化方法，导致计算成本过高，无法体现海洋生态学中至关重要的二阶依赖关系和动态季节性特征（如Ratmaya等人，2019）。Lewis-Borrell等人（2023）使用RStan为该模型提供了更为灵活的贝叶斯版本，支持二阶时间依赖关系，但对于长时间序列而言，其计算时间仍然过长。最后，Petris将动态线性模型视为状态空间模型的一种特例（Petris等人，2009），认为该模型在生态研究中具有应用前景（Levy等人，2014；Auger-Méthé等人，2021）。虽然Petris（2010）开发了用于贝叶斯动态线性模型分析的“dlm”R软件包，但目前尚不存在用于分析截断时间序列的相应操作方法。
集成嵌套拉普拉斯近似法（INLA，Rue等人，2009）是一种用于近似贝叶斯推断的方法。由于其速度较快且可通过R-INLA软件包轻松使用，近年来它已逐渐成为马尔可夫链蒙特卡洛等方法的一种替代方案（Lindgren和Rue，2015）。Ruiz-Cárdenas等人（2012）提出了一个基于INLA的计算框架，用于在线性模型和广义动态线性模型中实现近似完全贝叶斯推断。Ravishanker等人（2022）则提供了使用R-INLA软件包在动态贝叶斯框架下分析各类时间序列的指南。相较于“dlm”软件包，R-INLA在模型设定方面具有更大灵活性，不过代码实现的复杂度也相应更高。具体而言，R-INLA允许用户指定截断观测值，包括那些低于量化限制值的数据。本文旨在在Courbariaux等人（2022）和Lewis-Borrell等人（2023）的工作基础上，利用R-INLA软件包开发一种新框架，从而扩展传统上在动态线性模型中实现的结构性时间序列分析方法，使其能够妥善处理截断数据。该方法需能够在合理时间内自动处理成千上万的时间序列。此外，本文还特别关注如何评估模型假设的符合程度以及衡量模型的拟合优度。
章节节选
方法、真实数据与模拟
对监测数据的时间序列分析是使用动态线性模型（DLMs，West和Harrison，1997）进行的。作为状态空间模型的一种特例（Durbin和Koopman，2012），动态线性模型以线性方程和正态分布为特征，其核心目的是对时间序列进行结构分析，通常是通过将时间序列分解为趋势成分（即经过季节性调整后的趋势）、季节性成分以及观测误差项。这些方法的有效性
模拟时间序列
如图2所示，在100次模拟过程中，估计得到的超参数——包括观测值、演变参数（图中的斜率）、一阶和谐波方差及二阶和谐波方差，此处以标准差形式呈现——的平均值都与用于模拟数据的预期值十分接近。在这些具有现实特征（就海洋环境数据而言）的模拟数据集中，有一组数据后来被选为“模拟数据”，用于进一步
在模拟数据与真实数据上的应用
将CSTM应用于模拟数据后获得了大量有用信息。观测值方差和演变误差的估计值非常重要，因为它们与信噪比密切相关，而信噪比决定了模型各组成部分的适应性。因此，确定先验参数尤为重要。除RW2过程的参数外，其余参数均使用了R-INLA中的默认值：斜率方差的估计值则采用了从动态线性模型中得出的结果。这个数值，
结论与展望
通过这项工作，我们开发并验证了截断季节性趋势模型（CSTM），这是一种专为分析具有截断特征和季节性特征的时间序列而设计的状态空间模型。尽管该模型最初是为与海洋环境相关的时间序列开发的，但它可应用于任何需要将截断时间序列分解为趋势和季节性成分的领域。通过模拟的截断季节性时间序列以及真实的
CRediT作者贡献说明
Dominique Soudant：写作——审稿与编辑、写作——初稿撰写、可视化、方法论、数据整理、概念构建。Marie Courbariaux：写作——审稿与编辑、写作——初稿撰写、可视化、方法论。Yann Reynaud：写作——审稿与编辑、写作——初稿撰写、项目管理、资金筹集、数据整理、概念构建。
利益冲突声明
作者声明，他们不存在任何可能影响本文研究结果的已知财务利益冲突或个人关系。
致谢
感谢参与本项目开发的GIS Obépine团队，以及设计该工具、开发R-INLA软件包并在R-INLA论坛上提供用户支持的团队。同时也要感谢Nalini Ravishanker、Balaji Raman和Refik Soyer，他们解答了我们关于在R-INLA中实现动态线性模型的疑问。我们还要特别感谢Carlo Gaetan，是他让我们了解了分位数残差的概念。此外，也要感谢负责协调REMI项目的人士，以及那些采集样本、进行测量并录入数据的工作人员。
Dominique Soudant|Marie Courbariaux|Yann Reynaud
IFREMER, COAST, Nantes, F-44000, France
