监测预警基础45-这次学透 ARIMA之第6讲：ARMA建模R语言实操

系列简介

这是我们一系列原创技术贴，从易到难，每天学习一点。所有内容均为疾控数据分析、科研论文相关，或者说很多和现在的热门监测预警相关，所以我们这个系列就叫“监测预警基础”。

今天是第45节，这算是我们一个小专题吧，大概11讲，我们细化，短篇化，彻底把ARIMA模型搞懂。

今天是学透ARIMA专题第6讲的内容，R语言实操。

前面五篇我们把 ARIMA 的核心概念、建模流程全部讲透了，很多同行留言说：理论都懂了，就差上手跑一遍代码。

这一篇我们直接进入实操环节，全程用疾控最常用的周流感病例数做案例，所有代码逐行注释，复制粘贴就能运行。

哪怕你是 R 语言零基础，跟着步骤走，也能跑出自己的第一个 ARIMA 预测模型。

提前做好 2 项准备，就能直接跟着操作：

第一，安装好 R 与 RStudio（第 1 篇有安装地址，免费开源）；安装并加载必备分析包，打开 RStudio 在控制台运行以下代码即可：

第二，准备好数据，比如某传染病周数据

library(tseries)    # 时间序列基础工具包

我们用 Excel 格式的手足口病周数据演示，数据放在 R 工作目录下，或填写完整文件路径：

#dataset <- read.csv("文件名称，比如流感周监测数据.csv")

ARIMA 建模必须用ts时间序列格式，普通数据框不能直接建模，两个核心参数：

• start：数据起始时间，周度数据写c(年份, 起始周次)
• frequency：周期长度，周度数据设 52，月度数据设 12

     col = "#1E90FF")

拿到图先做业务判断：有没有年度流行高峰？整体趋势是升是降？和你的工作认知是否一致，避免数据导入顺序错误，我们运行完以后大概是下面这个样子：

稳是 ARIMA 建模的前提，我们用 ADF 单位根检验，先检验原始序列，不平稳再做差分。

# adf.test(dif2_hfmd)

我们的检测结果实质上是平稳的，不需要差分了：

（1）第一种方式： 绘制 ACF/PACF 图，手动定阶
tsdisplay()函数一行代码同时输出「时序图 + ACF 自相关图 + PACF 偏自相关图」
（2）那么，怎么看图定 p 和 q？
必须先明确两个基础概念：截尾是到某一阶之后，所有柱子突然全部落入两条虚线（95% 置信区间）里，相关性直接消失。拖尾是柱子慢慢衰减、波动下降，不会突然掉到 0 里，一直拖着尾巴
（3）定阶核心规则
AR (p) 模型：PACF 图 p 阶截尾，ACF 图拖尾 → p 就是 PACF 最后一个显著超出虚线的阶数。MA (q) 模型：ACF 图 q 阶截尾，PACF 图拖尾 → q 就是 ACF 最后一个显著超出虚线的阶数。ARMA (p,q) 模型：ACF 和 PACF 都是拖尾 → p 和 q 都大于 0，需要组合模型。
（4）新手看图提醒
真实疾控数据很少有教材里那样完美的截尾，只要能看出大致范围就可以，比如 PACF 前 2 阶显著，第 3 阶开始基本落入区间，那 p 大概率是 1 或 2，不用死抠完美截尾。
6.2 自动定阶：AIC 准则选最优模型
手动定阶给出范围后，我们用auto.arima()函数，基于 AIC 准则自动筛选最优的 p、d、q 组合，是日常工作最高效的方法：

best_model

我们看看我们的运行结果

左下 ACF 自相关图：前几阶自相关系数很高，随后缓慢衰减、呈波浪式波动，属于拖尾特征；每隔约 52 阶（也就是 1 年）会出现一个明显的高峰，有清晰的周期性波动。ACF 拖尾是 AR（自回归）模型的典型特征。
右下 PACF 偏自相关图：前 3 阶的偏自相关系数显著超出 95% 置信区间（两条虚线），第 3 阶之后，所有柱子基本都落入置信区间内，快速衰减到 0 附近，属于3 阶截尾

best_model 输出结果逐行详细解释

该对结果详细解释如下：

拟合完模型第一步不是预测，而是验证模型合不合格，先做基础的残差诊断：

Box.test(residuals(final_model), type = "Ljung-Box", lag = 10)

我们看一下我们的运行结果：

第一张图：Standardized Residuals（标准化残差时序图）

这是残差的原始波动走势图，看 3 个核心点：

第一，均值是否稳定：你的残差整体围绕 0 水平线上下随机波动，没有持续的上升 / 下降趋势，也没有固定的周期性高峰低谷，符合白噪声「均值恒定为 0」的特征。

第二，方差是否稳定：整体波动幅度基本一致，只有个别时间点出现小幅尖峰（对应手足口流行高峰的波动），没有出现某一阶段暴涨暴跌的情况，方差整体稳定。

第三，异常值情况：没有极端异常的离群点，残差都在合理范围内。

最后这张图结论：残差的基础波动特征符合平稳白噪声的要求。

第二张图：ACF of Residuals（残差的自相关图）

这是判断残差有没有自相关的核心图，判断规则非常明确：除了滞后 0 阶（lag=0）之外，所有滞后阶数的柱子，全部落在两条蓝色虚线（95% 置信区间）以内，就说明残差没有显著自相关性。

我们的图里，lag=1 到 lag=10 + 的所有柱子，全部都在蓝色虚线范围内，没有任何一阶显著超出置信区间，所以这张图结论是残差不存在显著的自相关性，序列里的自回归规律已经被模型提取干净了。

第三张图：p values for Ljung-Box statistic（Ljung-Box 检验 P 值图）

这是白噪声检验的统计验证图，对应不同滞后阶数的检验 P 值，横轴是滞后阶数（lag），纵轴是对应阶数的检验 P 值。判断标准只有一个：图中所有的点，全部都在 0.05 的蓝色虚线以上（也就是 P 值 > 0.05），就通过白噪声检验。

我们的图中所有点的 P 值都在 0.8 以上，远高于 0.05 的阈值，说明从滞后 1 阶到滞后 10 阶，残差都不存在显著的自相关性。结论就是统计层面验证残差为白噪声，模型信息提取充分。

模型通过基础检验后，就可以做预测了，周度传染病数据建议预测 1-4 周，最多不超过一个季度：当我们举例为了作图明显，预测了40周。

     col = "#ff1e1e")

看不看我们最终的结果，准确不准确姑且不论，我们只讲怎么看：

右边的图其实很清楚，左边的数值每一列分别是什么呢？

第1列：时间轴，对应预测的周次，匹配监测周期

第2列：点预测值，模型预测的该周病例数的「最可能数值」，是预测的核心结果，作为趋势预判的核心参考值.

第3-4列：80% 置信区间下限 / 上限，有 80% 的概率，该周的真实病例数会落在这个区间内，区间范围更窄。

第5-6列：95% 置信区间下限 / 上限列，有 95% 的概率，该周的真实病例数会落在这个区间内，区间更宽、可信度更高，传染病预警最常用：把 Hi 95（95% 置信区间上限）作为预警阈值，真实病例超过上限就提示发病超出正常波动范围，可能触发预警。

另外需要特别说明，结果中出现负数是正常现象，ARIMA 是线性统计模型，默认没有「病例数不能为负数」的约束，当预测期变长、不确定性变大，区间下限就会跌破 0；此外，我们用的是无季节性的基础 ARIMA，没有捕捉手足口病的年度周期规律，模型对低谷期的预测偏差更大，更容易出现负数。

实际应用处理方法就是直接把负数的下限替换为 0 即可，病例数不可能为负；后续换成 SARIMA 季节性模型后，预测区间会更贴合真实发病规律，负数出现的概率会大幅降低。

完整代码如下：

     col = "#ff1e1e")

跑通第一个模型只是开始，模型合不合格、预测精度够不够、能不能用于预警，必须靠系统的诊断和评价来判断。

下一篇我们会专门讲ARIMA 模型诊断与效果评价，拆解残差白噪声检验、训练集 / 测试集误差评价的完整方法，以及诊断不通过的问题排查方案，确保你建出来的模型可信、可用。

参考：

《时间序列分析-基于R》. [M] .王燕.中国人民大学出版社出版

传染病预测预警技术及实践案例分析. [M]. 杨鹏, 王小莉. 人民卫生出版社

编辑：普通疾控人 | 审核：诗酒趁年华

文章来源 | 原创

说明 | 转载只为分享，如有侵权联系删除

©版权声明 | 部分信息和图片来自公开网络

转载请注明

再次转载请注明出处

科普健康 | 宣传疾控

本号为多位疾控机构从业者运营

重点关注国内外健康事件

致力于疾控科普

在做好科普服务大众的同时

做好疾控机构的宣传

让更多的人了解疾控，拥抱健康

欢迎加「小编」微信(cdcjkr126com)

本文具体说明

本文为原创内容，文章为个人理解所学，不涉及疫情信息及内部保密数据，发表的目的为自我总结及给有需求的人士学习使用。如有不妥之处，欢迎联系小编修改、删除。

更多精彩视频，尽在“CDC疾控人”视频号