本讲概述基于趋势预测的时间序列分析方法、时间序列变化的影响因素、时间序列分析的步骤、统计软件和应用案例。
版权归中华医学会所有。
未经授权，不得转载、摘编本刊文章，不得使用本刊的版式设计。
除非特别声明，本刊刊出的所有文章不代表中华医学会和本刊编委会的观点。
本讲概述基于趋势预测的时间序列分析方法、时间序列变化的影响因素、时间序列分析的步骤、统计软件和应用案例。
时间序列是按照时间顺序排列的、随时间变化且相互关联的纵向数据。对时间序列进行观察、研究、找寻其发展变化的规律，预测其未来的走势就是时间序列分析[1]。常用的时间序列分析方法包括自回归滑动平均混合模型（Autoregressive Integrated Moving Average model，ARIMA模型） [2]、泊松回归广义可加模型[3]、分布滞后非线性模型（distributed lag nonlinear model，DLNM）[4]、灰度预测模型（Grey Model，GM模型）[5]以及基于机器学习算法的模型[6]等。其中ARIMA模型是在医学研究领域基于趋势预测的时间序列分析模型。该模型综合考虑了趋势变化、周期变化及随机干扰，展现出了较好的预测效果[7]。时间序列的基本特点有：假设事物发展趋势会延伸到未来；预测所依据的数据具有不规则性；不考虑事物发展之间的因果关系；时间序列数据用于描述现象随时间发展变化的特征。时间序列分析的核心是通过发现既往数据中的变化规律，并利用该规律对未来数据做出估计。时间序列数据每个时刻的值都与之前时刻的值存在着联系，表明过去的数据已经暗示了现在或者将来数据发展变化的规律，这种规律主要包括了趋势性、周期性和不规则性[8]。时间序列分析的主要目的是进行趋势预测，即确定已有的时间序列的变化模式，并假定这种模式会延续到未来，根据已有的时间序列数据预测未来的趋势变化。
我们以医学研究领域最常见的ARIMA模型为例，介绍时间序列分析的步骤。
用观测、调查、统计、抽样等方法取得被观测系统时间序列动态数据，建立用于时间序列分析的数据集。可以将数据集分为两部分：一部分作为训练集用作模型建立，另一部分作为验证集用作模型评估。
根据动态数据作相关图，进行相关分析，求自相关函数和偏自相关函数。相关图能显示出变化的趋势和周期，并能发现跳点和拐点。跳点是指与其他数据不一致的观测值。如果跳点是正确的观测值，在建模时应考虑进去，如果是反常现象，则应把跳点调整到期望值。拐点则是指时间序列从上升趋势突然变为下降趋势的点。如果存在拐点，则在建模时必须用不同的模型去分段拟合该时间序列，例如采用门限回归模型。ARIMA模型的表现形式一般为ARIMA（p，d，q）×（P，D，Q）s，其中7个模型参数p和P分别为非季节性和季节性自回归阶数，q和Q为非季节性和季节性移动平均阶数，d和D为非季节性和季节性差分次数，s为季节周期。我们可以根据所建立的数据集通过统计软件来确定模型参数。
时间序列可以分为平稳序列和非平稳序列。平稳序列是指基本上不存在长期趋势的序列，各观测值基本上在某个固定的水平上波动，或虽有波动，但并不存在某种规律，而其波动可以看成是随机的。或者说只含有随机波动的序列称为平稳序列。非平稳序列是指有长期趋势、季节性和循环波动的复合型序列，其趋势可以是线性的，也可以是非线性的。辨识时间序列平稳性的方法是进行平稳性检验。平稳性指的是期望不变，方差恒定，协方差不随时间改变，协方差只依赖于K这个时间跨度，不依赖于时间点t本身[9]。在建立模型前，需要观察并应用Augmented Dickey-Fuller（ADF）检验来检测时间序列的平稳性，对不稳定的序列做差分处理。待序列平稳后，观察其自相关函数（auto-correlation function，ACF）和偏自相关函数（partial auto-correlation function，PACF）的拖尾或截尾状况来识别模型的参数。自相关函数可以提供具有滞后值的任何序列的自相关值。简单来说，它描述了该序列的当前值与其过去的值之间的相关程度。偏自相关函数可以提供残差与下一个滞后值的相关性。截尾是指时间序列的自相关函数或偏自相关函数在大于某个常数k后快速趋于0为k阶截尾；拖尾是自相关函数或偏自相关函数始终有非零取值，不会在k大于某个常数后就恒等于零（或在0附近随机波动）。我们可以根据参数构建一个或数个模型，对比BIC值，选择BIC最小的模型为最优模型。
用最小二乘法计算最优模型的参数，再对模型的残差序列进行Ljung-Box检验，以检测其是否为白噪声序列，并用QQ图检测残差是否服从正态分布。合适模型的残差应满足均值为零的正态分布，而且任何滞后阶数的残差相关系数都为零。以平均绝对百分误差（Mean absolute percentage error，MAPE）为评价标准来评价所建模型用于趋势预测的准确度。MAPE值在5%~10%认为模型的预测效果高度准确，10%~20%认为模型的预测效果是好的，20%~50%认为模型是合理的，>50%认为模型的预测效果不准确。
时间序列分析方法需要采用统计软件进行分析处理。常用于处理时间序列数据的统计软件有Eviews（Econometrics Views），称为计量经济学软件包。使用EViews软件包可以对时间序列和非时间序列的数据进行分析，建立序列（变量）间的统计关系式，并用该关系式进行预测、模拟等等；SAS（Statistical Analysis System），是美国SAS软件研究所研制的一套大型集成应用软件系统，具有完备的数据存取、数据管理、数据分析和数据展现功能，SAS系统中可提供包括时间序列分析的多种分析功能；SPSS（Statistical Package for the Social Science），称为社会科学统计软件包，包含时间序列分析模块；S-plus软件，是基于S语言，并由MathSoft公司的统计科学部进一步完善，作为统计学家及一般研究人员的通用方法工具箱，S-plus强调演示图形、探索性数据分析、统计方法、开发新统计工具的计算方法，包含时间序列分析模块；Minitab软件，同样是国际上流行的一个统计软件包，其特点是简单易懂，在国外大学统计学系开设的统计软件课程中，Minitab与SAS并列，功能包括：基本统计分析、回归分析、方差分析、多元分析、非参数分析、时间序列分析、试验设计、质量控制、模拟、绘制高质量三维图形等，还具有许多统计软件不具备的功能--矩阵运算；R语言，R语言具有功能强大的程序包，在数据计算、统计分析以及数据挖掘等方面都广受好评，包含时间序列分析模块。
深圳市慢性病防治中心翁榕星研究团队于2021年在国际流行病学传染病学杂志第2期发表了"2005—2020年深圳市淋病疫情时间序列分析"论文[2]。该团队收集中国疾病预防控制信息系统里2005年1月至2020年11月年深圳市淋病月发病统计数据，将序列分为训练集和验证集，其中2005年1月至2020年5月作为训练集用作模型建立，2020年6至11月作为验证集用作模型评估。应用软件R 3.5.0版本建立ARIMA模型。最终确定BIC最小的模型为ARIMA（0，1，1）（2，1，1）12模型（BIC=370.51）。ARIMA（0，1，1）（2，1，1）12模型的参数检验均有统计学意义。因此，该模型可以用于验证集预测2020年6—11月深圳市淋病发病率。该模型预测结果显示2021年深圳淋病发病率将继续下降，区间为0.68/10万~3.70/10万。该模型预测水平较优，MAPE值为18.35%（MAPE值在10%~20%以内可认为模型是好的）。
广州医科大学冉丕鑫教授科研团队于2020至2021年采用时间序列分析方法，开展了多项省级生态研究课题[10,11,12]，旨在探索空气质量变化与因慢性阻塞性肺疾病恶化出现急性呼吸窘迫综合征入院的关系。该团队调查了医院慢性阻塞性肺疾病（AECOPD）急性加重期的入院人数与同时期污染物浓度的变化数据。研究数据包括颗粒物（即PM2.5、PM10和PM粗颗粒物）、二氧化硫的日浓度、2013年至2017年在广东21个城市的二氧化硫（SO2）、二氧化氮（NO2）、一氧化碳（CO）、臭氧（O3）和慢性阻塞性肺疾病入院人数。采用两阶段时间序列分析，估计年度归因分数、数量和直接住院率。研究结果表明：从2013年到2017年，二氧化硫、PM10和PM2.5的日均浓度分别下降了40%、30%、30%和26%。随着平均每日臭氧浓度显著增加，从2015年的103 μg/m3增加到2017年的152 μg/m3，均超过世界卫生组织的空气污染物排放目标，即100 μg/m3。随着污染物浓度的增加，滞后0~3 d，慢性阻塞性肺疾病患者因PM2.5污染入院的相对危险度为1.093（95%置信区间：1.06~1.13）；慢性阻塞性肺疾病患者因臭氧污染入院的相对危险度为1.092（95%CI：1.08~1.11）；慢性阻塞性肺疾病患者因二氧化硫污染入院的相对危险度为1.092（95%CI：1.05~1.14）。由于空气污染导致的慢性阻塞性肺疾病入院率从2013年的9.5%下降到2016年的4.9%，慢性阻塞性肺疾病住院人数下降可能与广东省PM2.5、PM10和二氧化硫浓度降低有关，而臭氧污染已成为慢性阻塞性肺疾病的一个重要的危险因素。
所有作者均声明不存在利益冲突
