第13章 流感预测

预测特定区域每周的流感发病率一直是重要的问题。传染病专家及全球安全专家均认为，传染病对人类福祉构成了重大危害。尤以流感为甚，它每年攻击全世界的易感人群，造成数百人死亡，其中大多数是幼儿和老年人。从公共卫生和国家安全的角度，根据特定季节的流感发展规律开展精准的预测模型至关重要。无论是对于预测流感，还是对于帮助研究人员探索传染病如何传播的通用理论，流感预测模型都是有用的。

案例研究：发生在大都市的流感疫情

我们将分析从2004年到2013年法国各行政区域的每周流感报告得到的数据集，来预测法兰西岛大区的流感发病率。

1.探索和清洗数据

我们从熟悉原始数据开始，先以表格形式来查看数据。

flu = fread('french_flu.csv')
flu[,flu.rate:=as.numeric(TauxGrippe)]
head(flu)

我们同时做一些基本的质量检查工作，比如查找感兴趣的变量中是否有NA值。虽然我们不知道NA值从何而来，但仍需要留意它。

# 计算flu.rate列中缺失值所占的百分比
nrow(flu[is.na(flu.rate)]) / nrow(flu)
# 识别哪些地区在flu.rate数据中存在缺失值
unique(flu[is.na(flu.rate)]$region_name)

上述代码的作用是： 1.计算flu.rate列中缺失值所占的百分比 2.在数据框flu中找出flu.rate列有缺失值的那些行，并提取这些行中的region_name列，然后返回这些region_name值中的唯一值。这通常用于识别哪些地区在flu.rate数据中存在缺失值。

总体来看，NA数据点的占比不是很高。此外，我们感兴趣的区域（法兰西岛大区）不存在NA值。

我们再执行一些数据清洗操作，从时间戳字段中分别获取周和年部分（时间戳字段目前是字符格式，不是数字格式或时间戳格式）

flu[,year := as.numeric(substr(week,1,4))]
flu[,wk := as.numeric(substr(week,5,6))]

我们添加了一个Date类型的字段，这样就可以在绘图时更好地设置时间轴，而不是讲数据视为没有时间戳。

flu[,date:=as.Date(paste0(as.character(flu$week),"1"),"%Y%U%u")]

警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2009 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid警告: (0-based) yday 368 in year 2004 is invalid

然后，我们查看与法兰西岛大区相关的数据，并按日期排序。

# 我们仅关注法兰西岛大区
paris.flu = flu[region_name == "ILE-DE-FRANCE"]
paris.flu = paris.flu[order(date,decreasing = FALSE)]

paris.flu[,.(week,date,flu.rate)]

如果留意行数，你也许会感到吃惊。如果说一年有52周，我们有10年的数据，为什么会有522行？我们原本期望得到520行（52×10=520）。此外，为什么会有两个NA日期？如果再看看原始数据，你就会得到答案。2004年和2009年似乎都有第53周。每隔几年就有一个53周而非52周的年份，这不是错误，而是公历系统的一部分。

我们首先要检查数据采样是完整且规律的，并确保每年有相同数量的数据点。

paris.flu[,.N,year]

我们可以看到，数据符合预期，也就是说，每年（除了刚才讨论的那两年）有52周，且每周有10个数据点（第53周除外）

现在我们已经考虑了数据的时间戳，接下来要检查时间序列的实际值（到目前为止，我们只考虑过时间索引）。是否存在某种趋势或季节性？让我们找找看。

paris.flu[,plot(date,flu.rate,
                type = "l",xlab = "Date",ylab="Flu rate"
                )]

NULL

透过简单的折线图，可以清楚地看出存在明显的季节性。虽然上图展示了明显的季节性成分，却没有显示出在季节性之外存在时间偏移。

季节性行为使第53周的情况变得复杂。如果我们想拟合一个季节性模型，就需要根据一年中的周数来定义季节性，并且季节长度不可变（第3章讨论过，这是区分周期和季节的依据）。尽管我们可以设想一些解决第53周问题的创造性方案，但本例将做简单处理，即清除这些数据。

paris.flu <- paris.flu[week != 53]

2.拟合季节性差分自回归移动平均模型

由于数据具有很强的季节性，我们首先考虑对数据拟合一个季节性差分自回归移动平均（SARIMA）模型。在本例中，由于数据每周采样一次，因此数据的周期是52.我们想选择一个简洁的模型，这是因为我们的时间序列不是特别长，只有520个数据点。

这个例子很好地说明，如果过于依赖自动化方法拟合模型，那我们将误入歧途。举例来说，我们可能首先考虑是否对数据进行差分，因此可以绘制流感发病率的自相关函数图和差分后的时间序列的自相关函数图，如下所示：

acf(paris.flu$flu.rate)

acf(diff(paris.flu$flu.rate,52))

流感发病率和差分后的时间序列的自相关函数图。我们只关注有限范围的滞后值。

如果我们粗心一些，可能会庆幸自己通过一次差分就解决了这个时间序列的平稳性问题。但这一点儿也说不通。这是每周的数据，并且我们观察到数据具有很强的季节性。为什么我们在自相关函数图中无法看到这一点？这是因为我们使用了acf()函数的默认参数。这些参数没有为我们展示足够远的滞后距离来让我们观察到季节性影响，而已知季节性开始于滞后值52（一年）。让我们用一个足够大的窗口来重新执行acf()，如下所示。

acf(paris.flu$flu.rate,lag.max = 104)

acf(diff(paris.flu$flu.rate,52),lag.max = 104)

流感发病率和差分后的时间序列的自相关函数图。现在我们扩大了滞后值范围。

这让我们对时间序列的自相关性有了更真实的了解。正如我们所看到的，在不同的滞后值处存在很强的自相关性。对于生活在四季变换的世界中的我们来说，这张图才是合理的。而流感发病率与临近的几周有很强的相关性，也就是与测量时间接近的几周。

考虑到季节性，流感发病率与滞后值为52或104的时间周期也有很强的相关性，这表示每年的季节性。不过，流感发病率也与滞后量为中间值的时间周期（比如半年，即26周）有相当强的关系，因为这种滞后还与季节差异和可预测的天气变化有关。举例来说，我们知道在半年内流感发病率可能会有很大的变化。如果早些时候流感发病率很高，现在应该很低，且反之亦然，这同样是由于季节性。

然后我们查看差分序列。我们可以看到，大量的时间序列自相关性已经被削弱。不过，在取值范围内仍然存在明显的自相关性，不仅是52或104（一年或两年），中间值也是如此。

尽管我们可能想继续差分，但需要记住，现实世界的数据永远不可能完全拟合季节性差分自回归移动平均模型。但我们要寻求最合理的方法来对数据建模，可以考虑根据季节性进行差分，或者采取不同的策略，指定某个时间范围进行差分。如下所示：

plot(diff(diff(paris.flu$flu.rate,52),52))

plot(diff(diff(paris.flu$flu.rate,52),1))

以两种方式对时间序列进行差分，以了解数据的季节性行为。

虽然两个结果都不理想，但相比之下，后者（季节性差分又再进行一阶差分）更令人满意。

关于拟合的决定或参数的选择是一个判断问题，与应用测试的问题类似。我们对清楚观察到的季节性赋予权重，但也不要使模型过于复杂或难懂。因此，在SARIMA(p,d,q),(P,D,Q)模型中，我们将分别采用d=1和D=1。然后，我们为标准ARIMA参数p和q选择AR参数和MA参数。使用以下代码，我们来进行标准可视化。如下所示：

par(mfrow=c(2,1))
acf (diff(diff(paris.flu$flu.rate,52),1),lag.max = 104)
pacf (diff(diff(paris.flu$flu.rate,52),1),lag.max = 104)

我们选择的差分序列的部分自相关函数图。

我们的数据集规模有限，如果使用太简单的模型就容易出错。PACF模型表明AR(2)模型是合适的，因此我们简单地将数据建模为SARIMA(2,1,0),(0,1,0) *********上面这步是怎么得到的？**********

我们想知道，如果不断地让模型拟合新数据（大多数为现实世界构建的模型正是如此），那么这个模型的性能会如何。也就是说，如果我们在几年前就开始对流感数据建模，每周只用当时可用的数据建模，那么模型的效果会怎样？我们通过前向滚动拟合和评估的过程来回答这个问题。