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 decompositionfit <- 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)fit2 <- HoltWinters(myts, gamma=FALSE)fit3 <- HoltWinters(myts)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 modelfit4 <- ets(myts)# Automated forecasting using an ARIMA modelfit5 <- auto.arima(myts)# 预测后三个月的数值forecast(fit4,3)forecast(fit5,3)
不同模型的预测结果是一致的！
关于时间序列分析的内容就先简单讲到这里，它不是我们进阶阶段的重点内容，感兴趣的朋友可以自学相关原理。