肾综合征出血热（hemorrhagic fever with renal syndrome，HFRS）是一种因接触感染汉坦病毒（Hantavirus，HV）的啮齿动物排泄物而引起的疾病^[1]。主要宿主动物为啮齿类动物，如褐家鼠（Rattus norvegicus）、黑线姬鼠（Apodemus agrarius）。感染汉坦病毒后部分患者出现发热、休克、出血、急性肾损伤等症状，严重者可导致并发症、多器官损害，甚至危及生命，病死率可达12%^[2]。及时预警有助于及早发现和控制疫情，降低疾病负担。疾病预测预警模型包括统计模型、机器学习模型、深度学习模型、多源数据融合模型和集成模型等类型，其中自回归差分移动平均（autoregressive integrated moving average，ARIMA）模型、季节性差分自回归移动平均（seasonal autoregressive integrated moving average，SARIMA）模型、指数平滑法（exponential smoothing，ETS）模型等为应用较为成熟的时间序列统计模型，能够解释发病与时序的关系。Zhang等^[3]运用差分自回归移动平均模型和支持向量积模型（ARIMA-SVM）的模型纳入平均气温、相对湿度、降雨量和日照时数等4个气象因素预测HRFS发病率，预测效果优于单一ARIMA模型。刘天等^[4]利用7种时间序列模型对全国HFRS发病率进行预测，结果显示指数平滑空间状态模型（the exponentially smoothed spatial state model，TBATS）为最优预测预警模型，适于优化HFRS预警模型。然而受疾病传播途径、地域特征、疾病相关因素（如媒介因素、社会因素和气象因素等^[5]）影响，预测结果间差异较大。基于区域特点建立区域内预测模型在HFRS防控中有重要意义，本研究拟探索基于最小绝对值收缩与选择算子算法回归联合引入自变量的季节性差分自回归移动平均模型（least absolute shrinkage and selection operator-seasonal autoregressive integrated moving average with exogenous regressors，LASSO-SARIMAX）^[6-7]、SARIMA模型和ETS模型对广州市2006－2022年HFRS发病中的预测效果，为制定区域内科学合理的疾病防控策略和措施提供科学建议。

1 材料与方法 1.1 资料来源 1.1.1 病例数据

病例数据来源于中国疾病预防控制信息系统的传染病报告信息系统。包括病例基本信息（年龄，性别，职业，家庭住址）、发病日期、诊断时间和死亡日期等，仅统计2006－2022年地区归属为广州市病例。

1.1.2 宿主动物数据

来源于2006年1月－2022年12月的广州市疾病预防控制中心鼠类监测数据。采用鼠密度指标衡量宿主动物的月度变化。

1.1.3 社会经济学数据

数据来源于广州市统计局（http://tjj.gz.gov.cn/）。主要为广州市常住人口密度、绿地覆盖面积以及常用耕地面积。

1.1.4 气象数据

来源于中国气象数据网（http://data.cma.cn/），包括日平均气温、最高气温、最低气温、相对湿度、风速、降雨量、气压和日照时数，通过数据变换将日数据转换为月度数据。

1.2 方法

基于建模复杂程度和预测结果对防控的指导意义考虑，本研究选用ETS、SARIMA和LASSO-SARIMAX模型进行广州市2006－2022年间HFRS发病的预测研究。其中，ETS模型利用历史数据解释短期的趋势和波动效果较好，选作建模复杂程度较低的模型代表，SARIMA模型建模过程中考虑时序的季节性、自相关性等因素，在模型参数调整上自由度更高，选作建模复杂程度中等的模型代表，LASSO-SARIMAX模型可纳入自变量进行预测，预测结果一定程度上还可揭示自变量与疾病之间的关联，有助于寻找防控中的着力点，选作建模复杂程度较高的模型代表。

1.2.1 SARIMA模型

SARIMA模型一般形式为：SARIMA（p，d，q）（P，D，Q）s。其中，d、D分别表示一般差分、季节性差分阶数，p、q表示自回归阶数、移动平均阶数，P、Q表示季节性自回归阶数和季节性移动平均阶数，s表示季节周期长度，存在明显季节性的疾病监测月数据，一般取12^[8]。建模的关键在于确定上述参数，上述参数取值一般不超过5。本研究将p、q、P、Q最大值设为5，d、D最大值设为2，采取逐步测试方法，利用最小赤池信息量准则（Akaike information criterion，AIC）原则筛选最优模型^[9]。利用Box-LjungQ检验最优模型残差是否为白噪声序列。

1.2.2 ETS模型

可分为1、2、3次指数平滑模型，4个参数分别为水平分量（α）、趋势分量（β、ϕ）和季节分量（γ）。在实际应用中，根据数据的趋势性、季节性特点，结合模型拟合优度[如对数似然比（log-likelihood ratio，LLR）]等指标选择最优模型。在R 4.2.2软件中利用“forecast”包中ets函数自动建模，筛选最优模型，筛选指标选择LLR^[10]。

1.2.3 LASSO回归^[4]

通过在回归模型系数的求解公式中添加各系数的绝对值之和作为惩罚函数，以发挥收缩系数和筛选特征变量，从而得到稀疏化模型的作用。LASSO回归复杂度调整的程度由参数λ来控制，λ越大对变量较多的线性模型的惩罚力度就越大，从而最终获得一个变量较少的模型。

1.2.4 SARIMAX模型^[11]

SARMIAX模型为带输入变量的SARIMA模型，其构建思想是：假设响应序列y_t{}和输入变量序列（即自变量序列）x_1t{}，x_2t{}，…，x_kt{}均平稳，首先构建响应序列和输入变量序列的回归模型，使用SARIMA模型继续提供残差序列ε_t{}中的相关信息，得到最终模型。

将与HFRS发生密切相关的因素引入自变量，可提高预测模型的预测精准度。然而，在引入多个变量，特别是气象因素时，由于潜在的相关效应，可能导致拟合过度或拟合精度不足，基于LASSO回归对自变量予以筛选，可以有效地进行特征选择和系数压缩，减少模型的复杂度，提高泛化能力，具有进一步提升预测精准度的能力。

1.2.5 模型比较

根据AIC判定模型的拟合优度（AIC值相对较小的模型较好）。将模型拟合值与实际值进行对比，计算平均误差（mean error，ME）、均方根误差（root mean square error，RMSE）、平均绝对误差（mean absolute error，MAE）、平均百分比误差（mean percentage error，MPE）和平均绝对百分比误差（mean absolute percentage error，MAPE）用于对比模型拟合和预测情况；上述指标均为评价模型拟合和预测的常用指标，其中ME可以反映模型的整体偏差；RMSE和MAE消除了量纲的影响，可以反映模型的总体精度；MPE为模型预测值与真实值之间百分比差异的平均值，可以反映模型对于数据变化的敏感程度；MAPE为模型预测值与真实值之间绝对百分比差异的平均值，同样用于反映模型对于数据变化的敏感程度^[12]。同时考虑预测时长对预测准确性的影响，本研究中利用上述3种模型进行1~24期的月发病预测，通过MAPE比较模型的预测效果。以上模型均在R 4.2.2软件中完成。

2 结果 2.1 广州市2006－2022年HFRS发病率基本情况

2006－2022年间广州市共报告HFRS病例1 895例，年均发病率0.06/10万，整体发病率呈现3~5年周期型分布，其中2006－2009年间发病率持续降低，由2006年的1.12/10万下降至2009年的0.61/10万，2009年到2013年持续攀升至1.09/10万，随后，2017年呈现第3个高峰，发病率达到0.90/10万，2017年后，发病率持续下降，2022年发病率下降至0.26/10万。月份分布中，每年发病集中于1－4月，共报告873例病例，占总病例数的46.56%。广州市2006－2022年逐月病例数见图 1，将2006－2022年广州市HFRS发病趋势分解为趋势项、季节项目和随机项，结果显示2006－2022年间发病呈现双峰型分布，存在2013和2017年2个高峰，季节性周期明显，但每个周期内波动幅度变化较小，且不随时间趋势变化，随机波动时序内较为稳定。见图 2。

2.2 ETS模型

利用2006－2020年广州市HFRS发病数据作为训练集建立ETS模型，最终2次指数平滑模型为最优模型，残差项为相乘模型，无趋势项，季节项为相加模型。其中α=0.121，γ=1×10^-4，AIC=1 331.384。训练集的自相关函数（autocorrelation function，ACF）、ME、RMSE、MAE、MPE和MAPE分别为-0.015、-0.077、4.261、3.325、-22.481和45.066（表 1）。残差QQ图中残差点均围绕45°线分布，呈现正态分布（图 3）。利用建立的ETS模型分别进行步长为12和24月的预测，结果显示，预测12和24月的MAPE值分别为58.532和104.016。

2.3 SARIMA模型建立

利用R语言auto.arima函数自动拟合模型，训练集选用2006－2020年广州市HFRS发病数据，得到最优模型为SARIMA（2，0，0）（1，0，1）_[12]。AIC=1 055.68，贝叶斯信息准则（Bayesian information criterion，BIC）=1 074.84，训练集的ACF1、ME、RMSE、MAE、MPE和MAPE分别为0.002、4.353、3.320、-30.443和51.403（表 1）。整体残差呈现正态分布（图 3），Box-Ljung c2=0.069 3，P=0.792 3。利用建立的SARIMA模型分别进行步长为12和24月的预测，结果显示，预测12和24月的MAPE分别为124.572和194.269（表 1）。

2.4 SARIMAX模型建立 2.4.1 影响因素选择

基于样本量和预测因子数量，本研究选择用LASSO回归，在14个预测因子（年份、月份、鼠密度、平均气温、最高气温、最低气温、相对湿度、风速、降雨量、气压、日照时数、常住人口数、绿地覆盖面积和常用耕地面积）中筛选出有意义的预测因子组合，剔除共线性较强变量，并进行CV-LASSO验证。结果显示随着λ值增大，14个预测因子对应系数逐步被压缩至0，其中lgλ＜1后，仅有4个变量的系数未被压缩至0（图 4A）。交叉验证结果显示最优结果（距离MSE 1个标准误时）对应λ=0.022，lse=0.701（图 4B）。选取1se时预测因子数，最终留下4个预测因子用于建模，4个因子分别是最高气温、风速、日照时数和常住人口数。

2.4.2 LASSO-SARIMAX模型建立

利用2006－2020年广州市HFRS发病数据作为响应因子，同期最高气温、风速、日照时数和常住人口数作为自变量，建立LASSO-SARIMAX模型，通过对比AIC，选择AIC最小的模型（p，d，q）（P，D，Q）作为最终模型参数，最终建立模型为（1，0，1）（1，1，1）_[12]，训练集的ACF1、ME、RMSE、MAE、MPE和MAPE分别为-0.021、-0.074、3.968、2.920、-18.332和39.466（表 1）。整体残差呈现正态分布（图 3），利用建立的LASSO-SARIMAX模型分别进行步长为12和24月的预测，结果显示，预测12和24月的MAPE分别为49.316和120.246（表 1）。模型参数中，风速的参数值有统计学意义（t=2.512，P=0.012），最高气温、日照时数和常住人口数均无统计学意义。见表 2。

2.5 3个预测模型间比较

对3个模型拟合数据集、12月预测和6月预测结果进行比较，结果显示除预测24月时低于ETS模型外，LASSO-SARIMAX模型训练数据集、预测12月的MAPE值均最低。对3个模型分别进行预测步长为1~24月的预测（LASSO-SARIMAX模型预测步长为2~24月），预测数据集的MAPE结果显示，在预测步长为6个月以内，3个模型之间差异不明显，预测时长超过9个月后，SARIMA模型的MAPE值陡然升高，在步长为20个月时，MAPE值超过100。ETS和LASSO-SARIMAX模型MAPE值相对较低，其中8~18个月的中期预测中，LASSO-SARIMA模型预测效果最优。见图 5、6。

3 讨论 3.1 广州市HFRS整体趋势

研究结果显示，广州市HFRS发病整体呈现下降趋势，这与居民生活水平改善、灭鼠、疫苗接种和健康宣教等干预措施的推广应用有关^[13]。HFRS表现出的3~5年流行周期可能与气象因素及媒介因素的周期性波动有关^[14]。2019年以后HFRS的发病持续降低，可能与新型冠状病毒感染疫情（新冠疫情）期间人群减少户外活动，降低了与鼠类等宿主动物的接触有关，随着野外作业、旅游等的恢复，HFRS的发病可能存在升高的潜在风险，提示应加强HFRS的监测。

3.2 不同模型间总体比较

本研究中建立的2次指数平滑ETS模型中无趋势项，与研究周期内，HFRS呈现双峰型，无递增或递减趋势有关。基于最小AIC选择的SARIMA模型与LASSO-SARIMAX模型中的差分项以及季节性差分项有所不同，提示建模时在数据整体趋势中，纳入影响因素时仍需在单纯时间序列分析中进一步优化模型参数^[15]。LASSO-SARIMAX模型自变量选择中未纳入鼠密度可能与后者受气象因素、经济因素、耕地面积等社会因素影响较大有关^[16]。

3.3 不同模型间拟合效果比较

本研究中模型拟合效果由高到低依次为LASSO-SARIMAX、ETS和SARIMA模型，其中SARIMA模型预测效果较差（随着预测步长增加，MAPE值快速突破100）与其他已有研究报告不同^[17]，分析其原因可能与区域内发病特征、研究时段跨度有关。LASSO-SARIMAX模型拟合效果最佳可能与通过LASSO回归对影响因子进行聚合压缩后纳入影响因素有关，该建模策略在通过影响因素提升对时序内发病变化解释度的同时又有效避免了过度拟合，使得其拟合效果优于单纯时间序列分析，目前尚无LASSO-SARIMAX模型应用于疾病预测的报道，李凤灵^[18]以居民消费价格指数为样本的预测模型选择中，通过LASSO回归方法所提取的变量进行ARIMAX的拟合，效果同样优于单纯ARIMA模型。单纯时间序列分析中，ETS的拟合效果优于SARIMA模型，可能与较少数据量情形下，SARIMA模型捕捉时间序列数据特征能力低于ETS模型有关。

3.4 不同模型间预测效果比较

从预测效果来看，12个月预测效果由高到低依次为LASSO-SARIMAX、ETS和SARIMA模型，预测24个月时ETS模型预测效果反超LASSO-SARIMAX模型，但两者MAPE值均超过100，LASSO-SARIMAX短期内预测效果更好，可能与纳入影响因素后，对异常事件（如大规模灭鼠活动、新冠疫情期间人群外出减少）造成的短期发病波动预测效果更好。在实际建模中，尽管LASSO-SARIMAX模型建模复杂程度较高，但建模结果可以一定程度上解释自变量对HFRS发病的影响，与ETS模型相比，可用于评估异常事件、防控措施等的影响，在HFRS防控中现实意义更高。

3.5 结论及局限性

综上所述，基于LASSO-SARIMAX模型在广州市HFRS的中短期预测中有着较好效果，具有一定的推广应用空间，在长期预测方面，仍需进一步纳入影响因素和联合模型参数。本研究也存在一定的局限性，首先本研究选取的常见时间序列模型中，仅LASSO-SARIMA模型应用到了影响因素，其应用效果仍需与广义线性模型和分布滞后非线性模型等纳入影响因子分析的模型进行比较。其次，数据选择中被研究纳入的影响因素考量仍较少，例如对HFRS影响较强的户外活动情况去进行量化衡量且受限于数据规整尺度，部分因素如人口数据和耕地数据等无法展示年度内变化，使得模型预测所需数据的可获取性和模型应用可及性受到一定程度影响。

利益冲突  无