R语言进阶之时间序列分析

时间序列分析虽然主要应用于经济领域，但它作为一种分析时间依赖性变量之间关系的重要方法，值得我们去学习。就像孟德尔随机化里的工具变量方法那般，虽然它起自计量经济学，但在流行病学和遗传学上得到了广泛应用，所以我们做研究时需要有学科交叉思维，学科交叉往往能带来突破。

在这一期内容中，我主要会和大家讲解时间序列数据的创建、季节性分解、指数模型与ARIMA模型。

1. 创建时间序列

R语言的内置函数ts()可将数值型向量转换成R里的时间序列对象，其使用形式如下

ts(vector, start=, end=, frequency=)

这里start是指第一个观测值的时间，也即起始时间；而end相应地就是最后一个观测值的时间，即结束时间。另外，frequency则是指单位时间的长度，比如1代表以年为单位计录数据, 4则代表季度，12则代表月。

# 随机创建一个72个月内某医院结核患者数目的向量
# 这里的数据是人为创建的，不代表真实情况
# 起始时间点为2009年1月，终止时间点是2014年12月
myvector <- c(200:271) # 生成72个数据
myts <- ts(myvector, start=c(2009, 1), end=c(2014, 12), frequency=12) # frequency=12代表以月为单位，start和end里的第一个数代表年份，第二个代表月数
# 利用window()函数提取从2014年6月到2014年12月这部分的时间序列数据
myts2 <- window(myts, start=c(2014, 6), end=c(2014, 12)) #start和end分别代表提取数据的起止点
# 绘制时间序列图
plot(myts)

时间序列图的横坐标代表的是时间，纵坐标代表的是观测值。

2. 季节性分解

一个季节性时间序列中会包含三部分，趋势部分、季节性部分和无规则部分，我们可以在R中使用stl()函数来对时间序列进行季节性分解。如果数据间是相乘的关系，我们可以利用log()函数来将其转换成加性的。

# Seasonal decomposition
fit <- stl(myts, s.window="period") # 季节性分解
plot(fit)

从上图我们可以看出：时间序列数据被分解成三部分，季节性部分（seasonal）、趋势部分（trend）和剩余无规则部分（remainder）。从图中可以看出数据是有一定季节性的（以年为单位重复波动），但是由于季节性数据比趋势小很多，我们其实可以忽略这个季节性。

# 使用forcast包绘图
library(forecast)
seasonplot(myts)

上图是将每一年的数据单独绘制在一张图上，比如最底端的直线代表2009年数据，最顶端代表2014年数据。

3.指数平滑模型

R语言的内置函数 HoltWinters()和“forecast”包的ets()都可以用来拟合指数模型，这里我们主要使用的是HoltWinters()函数。

# 一次指数平滑模型，代表数据没有趋势和季节性
fit1 <- HoltWinters(myts, beta=FALSE, gamma=FALSE) # beta代表趋势，gamma代表季节性
# 二次指数平滑模型，指数据有趋势但没有季节性
fit2 <- HoltWinters(myts, gamma=FALSE)
# 三次指数平滑模型，数据既有趋势又有季节性
fit3 <- HoltWinters(myts)
# 利用模型3预测后三个观测值
library(forecast)
forecast(fit3, 3)
plot(forecast(fit3, 3))

可以看出后三个月的预测值分别为272、273和274，与事实是相符的（每月增加一个感染者）。

4. ARIMA模型

ARIMA模型中文全称是自回归积分滑动平均模型（autoregressive integrated moving average），在R中我们可以使用“forecast”包的auto.arima()函数来实现该模型。

library(forecast)
# Automated forecasting using an exponential model
fit4 <- ets(myts)
# Automated forecasting using an ARIMA model
fit5 <- auto.arima(myts)
# 预测后三个月的数值
forecast(fit4,3)
forecast(fit5,3)

不同模型的预测结果是一致的！

关于时间序列分析的内容就先简单讲到这里，它不是我们进阶阶段的重点内容，感兴趣的朋友可以自学相关原理。