时间序列分析虽然主要应用于经济领域,但它作为一种分析时间依赖性变量之间关系的重要方法,值得我们去学习。就像孟德尔随机化里的工具变量方法那般,虽然它起自计量经济学,但在流行病学和遗传学上得到了广泛应用,所以我们做研究时需要有学科交叉思维,学科交叉往往能带来突破。
在这一期内容中,我主要会和大家讲解时间序列数据的创建、季节性分解、指数模型与ARIMA模型。
1. 创建时间序列
R语言的内置函数ts()可将数值型向量转换成R里的时间序列对象,其使用形式如下
代码语言:javascript复制ts(vector, start=, end=, frequency=)
这里start是指第一个观测值的时间,也即起始时间;而end相应地就是最后一个观测值的时间,即结束时间。另外,frequency则是指单位时间的长度,比如1代表以年为单位计录数据, 4则代表季度,12则代表月。
代码语言:javascript复制# 随机创建一个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()函数来将其转换成加性的。
代码语言:javascript复制# Seasonal decomposition
fit <- stl(myts, s.window="period") # 季节性分解
plot(fit)
从上图我们可以看出:时间序列数据被分解成三部分,季节性部分(seasonal)、趋势部分(trend)和剩余无规则部分(remainder)。从图中可以看出数据是有一定季节性的(以年为单位重复波动),但是由于季节性数据比趋势小很多,我们其实可以忽略这个季节性。
代码语言:javascript复制# 使用forcast包绘图
library(forecast)
seasonplot(myts)
上图是将每一年的数据单独绘制在一张图上,比如最底端的直线代表2009年数据,最顶端代表2014年数据。
3.指数平滑模型
R语言的内置函数 HoltWinters()和“forecast”包的ets()都可以用来拟合指数模型,这里我们主要使用的是HoltWinters()函数。
代码语言:javascript复制# 一次指数平滑模型,代表数据没有趋势和季节性
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()函数来实现该模型。
代码语言:javascript复制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)
不同模型的预测结果是一致的!
关于时间序列分析的内容就先简单讲到这里,它不是我们进阶阶段的重点内容,感兴趣的朋友可以自学相关原理。