关于算法:R语言分布滞后线性和非线性模型DLMs和DLNMs分析时间序列数据

62次阅读

共计 7151 个字符,预计需要花费 18 分钟才能阅读完成。

原文链接:http://tecdat.cn/?p=20953

序言

本文演示了在工夫序列剖析中利用散布滞后线性和非线性模型(DLMs 和 DLNMs)。Gasparrini 等人 [2010] 和 Gasparrini[2011]论述了 DLMs 和 DLNMs 的倒退以及工夫序列数据的实现。本文形容的示例涵盖了工夫序列数据 DLNM 办法的大多数规范利用,并探讨了 DLNM 包用于指定、总结和绘制此类模型。只管这些例子在空气污染和温度对衰弱的影响方面有具体的利用,但它们很容易被推广到不同的主题,并为剖析这些数据集或其余工夫序列数据源奠定了根底。

数据

示例应用工夫序列数据集(包含 1987-2000 年期间每日观测数据)摸索了空气污染和温度与死亡率之间的关系。在 R 会话中加载后,让咱们看一下前三个察看后果:

 date time year month doy dow death cvd resp temp dptp
1 1987-01-01 1 1987 1 1 Thursday 130 65 13 -0.2777778 31.500
2 1987-01-02 2 1987 1 2 Friday 150 73 14 0.5555556 29.875
3 1987-01-03 3 1987 1 3 Saturday 101 43 11 0.5555556 27.375
rhum pm10 o3
1 95.50 26.95607 4.376079
2 88.25 NA 4.929803
3 89.50 32.83869 3.751079

数据集由 1987-2000 年期间每天进行观测的序列组成。

示例 1:一个简略的 DLM

在第一个例子中,我指定了一个简略的 DLM,评估 PM10 对死亡率的影响,同时调整温度的影响。我首先为这两个预测值建设两个穿插基矩阵,而后将它们蕴含在回归函数的模型公式中。假如 PM10 的影响在预测因子的维度上是线性的,因而,从这个角度来看,咱们能够将其定义为一个简略的 DLM,即便回归模型也预计了温度的散布滞后函数,这是一个非线性项。首先,我运行 crossbasis()来构建两个穿插基矩阵,将它们保留在两个对象中。两个对象的名称必须不同,以便别离预测它们之间的关联。代码如下:

cb(pm10, lag=15, argvar=list(fun="lin",
arglag=list(fun="poly",degree=4 

在具备工夫序列数据的程序中,第一个参数 x 用于指定向量序列。在这种状况下,咱们假如 PM10 的影响是线性的(fun=“lin”),同时通过一个具备 5 个自由度的天然三次样条曲线(fun=“ns”,默认抉择)来模仿与温度的关系。外部结点(如果未提供)由 ns()搁置在默认的等距分位数处,而边界节点位于温度范畴处。对于滞后空间的基数,我用 4 次多项式函数(设置次数 =4)指定 PM10 长达 15 天的滞后效应(最小滞后默认为 0)。温度的滞后效应由两个滞后层(0 和 1 -3)定义,假如每个层内的效应为常数。参数 breaks= 1 定义了第二个区间的下边界。此类的办法函数 summary()提供了穿插基(以及二维中的相干基)的概述:

 CROSSBASIS FUNCTIONS
observations: 5114
range: -3.049835 to 356.1768
lag period: 0 15
total df: 5
BASIS FOR VAR:
fun: lin
intercept: FALSE
BASIS FOR LAG:
fun: poly
degree: 4
scale: 15
intercept: TRUE

当初,在回归模型的模型公式中能够蕴含这两个穿插基对象。在这种状况下,我拟合工夫序列模型,假如泊松散布,工夫的润滑函数,7 df/ 年(为了校对季节性和长时间趋势)和星期几作为因子:

glm(death ~ cb1.pm + cb1.temp + ns(time, 7*14) + dow,
family=quasipoisson()

通过上述模型预测的 PM10 与死亡率的特定程度的预计关联可通过函数 crosspred()进行总结,并保留在具备雷同类的对象中:

pred(cb1.pm, model1, at=0:20, bylag=0.2, cumul=TRUE) 

该函数包含用来预计参数的 base1.pm 和 model1 对象作为前两个参数,而 at = 0:20 示意必须为从 0 到 20 µgr / m3 的每个整数值计算预测。通过设置 by lag = 0.2,沿着滞后空间以 0.2 的增量计算预测。绘制后果时,此更精密的网格产生更平滑的滞后曲线。参数 cumul(默认为 FALSE)批示还必须包含沿滞后的增量累积关联。没有通过参数 cen 定义核心,因而默认状况下将参考值设置为 0(这种状况产生在函数 lin()上)。当初,这些预测已存储在 pred1.pm 中,能够通过特定的办法对其进行绘制。例如:

> plot(pred1, "slices",
main="与 PM10 减少 10 个单位的关联性")
> plot(pred1,ylab="累计 RR",
main="PM10 减少 10 个单位的累积关联")

该函数蕴含带有存储后果的 pred1.pm 对象,参数“slices”定义了咱们要绘制对应于相干维度中 predictor 和 lag 的特定值的关系图。var=10 时,我显示 PM10 特定值的滞后响应关系,即 10µgr/m3。该关联应用 0µgr/m3 的参考值来定义,从而为 10 个单位的减少提供预测特定关联。我还为第一个图抉择了不同的色彩。PM10 的特定值,即 10 µgr / m3。应用 0 µgr / m3 的参考值定义此关联,从而为减少 10 个单位提供了特定于预测变量的关联。我还为第一个图抉择了不同的色彩。参数 cumul 批示是否必须绘制以前保留在 pred1.pm 中的增量累积关联。后果如图 1a-1b 所示。置信区间被设置为参数 ci 的默认值“area”。在左面板中,其余参数通过 ci.arg 传递给绘图函数 polygon(),绘制暗影线作为置信区间。对这些曲线图的解释有两个方面:滞后曲线示意特定日期 PM10 减少 10µgr/m3 后将来每一天的危险减少(正向解释),或者过来每一天雷同 PM10 对特定日期危险减少的奉献(反向解释)。图 1a-1b 中的曲线图表明,PM10 危险的初始减少在较长的滞后工夫被逆转。PM10 在 15 天滞后工夫内减少 10 个单位的总体累积效应(行将所有奉献相加至最大滞后工夫)及其 95% 置信区间可通过 pred1.pm 中蕴含的对象 allRRfit、allRRhigh 和 allRRlow 提取,键入:

> pred1
10
0.9997563
> cbind(pred1.p
[1] 0.9916871 1.0078911

例 2:节令剖析

第二个例子的目标是阐明数据仅限于特定节令的剖析。这种剖析的独特之处在于,假如数据是由不同年份的多个等距有序的多个节令序列组成,而不是一个繁多的间断序列。在本例中,我应用第 3 节中曾经看到的雷同步骤,别离评估了臭氧浓度和温度对滞后 5 天和 10 天死亡率的影响。首先,我创立一个季节性工夫序列数据集,只取冬季的区间,并将其保留在数据框 ChicagonMaps 中:

Sseas <- subset(NMMAPS, month %in% 6:9) 

同样,我首先创立穿插基矩阵:

cb(o3, lag=5,
argvar=list(fun="thr",thr=40.3), arglag=list(fun="integer"),
group=year) 

参数组批示定义多个序列的变量。每个序列必须是间断的、残缺的和有序的。在这里,我假如 O3 的影响在达到 40.3 µgr / m3 之前为零,而后呈线性,并利用了高阈值参数化(fun =“thr”)。对于温度,我应用一个双阈值,并假如在低于 15°C 且高于 25°C 时成果是线性的,并且在两者之间为零。阈值应用自变量 thr.value(缩写为 thr)进行抉择,而未指定的自变量侧则将第一个穿插基准的默认值设置为“h”,将第二个穿插基准的默认值设置为“d”(给定)提供了两个阈值)。对于滞后维度,我为 O3 指定了一个不受约束的函数,为最多 5 天的每个滞后(fun =“integer”)利用一个整数(默认状况下,最小滞后等于 0)。对于温度,我定义了滞后 0 -1、2-5、6-10 的 3 个区间。回归模型包含一年中某天和某天的天然样条,以便别离形容每年的季节性影响和长期趋势。特地是,与以前的剖析相比,后者的自由度要小得多,因为它只须要捕捉安稳的年度趋势即可。除此之外,以雷同的形式进行预计和预测。代码为:

glm(death ~ cb2.o3 + cb2.temp + ns(doy, 4) + ns(time,3) + dow,
family=quasipoisson()) 

必须在其中指定要进行预测的值:在这里,我定义了 0 到 65 µgr / m3(大概是臭氧散布范畴)的整数,加上阈值和 50.3µgr/m3 的值,对应于阈值以上 10 个单位的减少。向量将主动排序。将主动抉择由 thr()建模的参考 - 反馈曲线,并且能够不定义核心参数。我绘制了 O3 减少 10 个单位的预测因子特定滞后反馈关系,但置信区间为 80%,并且还绘制了总体累积裸露反馈关系。

> plot(pred2.o3, "slices", main="滞后响应 超过阈值 10 个单位"(80 置信区间))
> plot(pred2.o3,"overall",xlab="臭氧", ci="l", main="5 个滞后的总体累积关联")

在第一个语句中,参数 ci =“bars”批示与图 1a-1b 中的默认“区域”不同,置信区间用条形图示意。此外,参数 ci.level = 0.80 指出绘制 80%的置信区间。最初,我依据参数类型和 pch 抉择了带有特定符号的点。在第二个语句中,参数 type =“overall”示意必须绘制整体累积关联,置信区间为线,ylim 定义 y 轴的范畴,lwd 示意直线的厚度。与后面的示例相似,通过 ci.arg 指定的参数列表来欠缺置信区间的显示,在这种状况下,将其传递给低级函数 lines()。与上一个示例相似,咱们能够从 pred2.o3 中提取臭氧浓度超过阈值(50.3−40.3µgr/m3)10 个单位时的预计总体累积效应,以及 95% 置信区间:

> pred2.o3$allRRfit["50.3"]
50.3
1.047313

> cbind(allRRlow, allRRhigh)["50.3",]
[1] 1.004775 1.091652

能够将雷同的图和计算利用于温度的冷热效应。例如,咱们能够形容超过低阈值或高阈值 1°C 的危险减少。用户能够反复上述步骤执行此剖析。

示例 3:二维 DLNM

在后面的例子中,空气污染(别离为 PM10 和 O3)的影响被假设为齐全线性或高于阈值的线性。这一假如有助于解释和示意这种关系:从不思考预测因子的维度,并且很容易绘制出 10 个单位减少的特定或总体累积关联。相同,当思考到温度的非线性相关性时,咱们须要采纳二维透视图来示意沿预测变量空间和滞后量非线性变动的关联。在此示例中,我指定了一个更简单的 DLNM,其中应用两个维度的平滑非线性函数来预计相关性。只管关系的复杂性更高,但咱们将看到指定和拟合模型以及预测后果所需的步骤与之前看到的简略模型完全相同,只须要抉择不同的绘图即可。用户能够采纳雷同的步骤来钻研先前示例中的温度影响,并扩大 PM10 和 O3 的图。在这种状况下,我运行 DLNM 来钻研温度和 PM10 对死亡率的影响,别离达到滞后 30 和 1。首先,我定义了穿插基矩阵。特地是,温度的穿插基是通过天然和非天然样条曲线指定的,应用来自软件包样条曲线的函数 ns()和 bs()。代码如下:

 > varknots <- equalknots(temp,fun="bs",df=5,degree=2)
> lagknots <- logknots(30, 3) 

预测空间的抉择基函数是 PM10 效应的线性函数和温度 5 自由度的二次 B 样条(fun=“bs”),通过函数 equalknots()抉择,默认状况下,节点搁置在预测器空间中的等间距值。对于滞后空间,我假如 PM10 的简略滞后 0 - 1 参数化(即直到滞后 1 的单个层,最小滞后默认等于 0,放弃默认值 df=1),而我定义了另一个三次样条曲线,这一次温度滞后维度具备天然束缚(fun=“ns”默认)。应用函数 logknots(),将滞后样条曲线的节点搁置在滞后对数比例中的等间距值处。温度和死亡率之间关系的预计、预测和绘图通过以下形式进行:

 > plot(pred3.temp, xlab="温度" lphi=30,
main="温度效应的 3D 图")
> plot(pred3.temp, "contour", xlab="温度",
plot.title=title("等高线图",xlab="温度",ylab="滞后")) 

请留神,此处的预测值以 21°C 为核心,该点代表解释预计效应的参考。这里须要执行此步骤,因为该关系是应用没有显著参考值的非线性函数建模的。仅在 crosspred()中应用参数 by = 1 来抉择值,这些值定义了预测变量范畴内的所有整数值。第一个绘图表达式生成一个 3D 绘图,如图 3a 所示,其中通过参数 theta-phi-lphi 取得了非默认的视角选项。第二个绘图表达式指定图 3b 中的轮廓图,其中题目和轴标签由参数 plot.title 和 key.title 抉择。图 3a-3b 中的曲线图提供了二维裸露 - 滞后 - 反馈关联的综合总结,但其在预测值或滞后的特定值下提供关联信息的能力无限。此外,因为三维图和等高线图中未报告预计关联的不确定性,因而它们也仅限于推理目标。通过绘制特定预测值和滞后值的效应面“切片”来提供更具体的剖析。代码是:

> plot(pred3.temp, "slices", var=-20, 
main="不同温度下的滞后反馈曲线,参考 21C")
> for(i in 1:3) lines(pred3.temp, "slices", var=c(0,27,33)[i]
> legend("topright",paste("温度 =" 

后果如图 4a-4b 所示。图 4a 阐明了特定于 -20℃、0℃、27℃和 33℃的温和和极其冷热温度(参考 21℃)的滞后反馈曲线。图 4b

形容了滞后 0 和滞后 5 的裸露 - 反馈关系(左列)以及温度 -20℃和 33℃下的滞后 - 反馈关系(右列)。参数 var 和 lag 定义了要在图 3a-3b 中的成果外表上切割的“切片”的温度和滞后值。第一个表达式中的参数 ci =“n”示意不能绘制置信区间。在多面板图 4b 中,列表参数 ci.arg 用于绘制置信区间,将其作为暗影线减少灰色对比度,在此处更加显著。初步解释表明,高温比低温具备更长的死亡危险,但不是立刻的,在滞后 0 时显示出“爱护”效应。这种剖析能力很难用更简略的模型实现,可能会失落关联的重要细节。

=

示例 4:降维 DLNM

在最初一个例子中,我展现了如何应用函数 crossreduce()将二维 DLNM 的拟合度升高到由一维基的参数示意的摘要。首先,我指定一个新的穿插基矩阵,运行模型并以通常的形式进行预测

指定的温度穿插基由双阈值函数和天然三次样条组成,别离以 10°C 和 25°C 的截止点作为预测器的维数,以对数标度中相等间距的节点值作为滞后量,如前一示例所示。能够对 3 个特定的摘要进行归约,即总的累积,滞后特定和预测变量特定的关联。前两个代表裸露 - 反馈关系,而第三个代表滞后 - 反馈关系。代码如下:

credu(cb4, model4) 

在滞后 5℃和 33℃时,别离在两个空间中计算关联的缩小。“crossreduce”类的 3 个对象蕴含相干空间中一维基的批改缩减参数,可与原始模型进行比拟:

> length(coef(pred4))
[1] 10
> length(coef(redall)) ; length(coef(redlag))
[1] 2
[1] 2
> length(coef(redvar))

正如预期的那样,对于预测变量的空间,参数数量已缩小到 2(与双阈值参数化统一),对于滞后空间,参数数量已缩小到 5(与天然三次样条曲线基的维度统一)。然而,原始拟合和简化拟合的预测是雷同的,如图 5a 所示:

> plot(pred4, "overall", xlab="温度", ylab="RR",
ylim=c(0.8,1.6), main="整体累积关联")
> lines(redall, ci="lines",col=4,lty=2)
> legend("top",c("原始","降维"),col=c(2,4),lty=1:2,ins=0.1) 

这个过程也能够通过从新结构原始的一维基和预测给定修改参数的关联来说明。作为一个例子,我应用 onebasis()为滞后空间再现了天然三次样条曲线,并预测后果:

样条基是基于对应于滞后 0:30 的整数值计算的,节点与原始穿插基的值雷同,并且不居中,以截距作为滞后基的默认值。应用修改后的参数对 33℃的预测值进行计算。原始、简化和重建预测值的雷同拟合如图 5b 所示,由以下公式得出:

> plot(pred4, "slices", var=33, 
main="33°C 时特定于预测变量的关联")

> legend("top",c("原始","降维","重构"),


最受欢迎的见解

1. 在 python 中应用 lstm 和 pytorch 进行工夫序列预测

2.python 中利用长短期记忆模型 lstm 进行工夫序列预测剖析

3. 应用 r 语言进行工夫序列(arima,指数平滑)剖析

4.r 语言多元 copula-garch- 模型工夫序列预测

5.r 语言 copulas 和金融工夫序列案例

6. 应用 r 语言随机稳定模型 sv 解决工夫序列中的随机稳定

7.r 语言工夫序列 tar 阈值自回归模型

8.r 语言 k -shape 工夫序列聚类办法对股票价格工夫序列聚类

9.python3 用 arima 模型进行工夫序列预测

正文完
 0