关于数据挖掘:R语言分布滞后非线性模型DLNM研究发病率死亡率和空气污染示例附代码数据

全文下载链接：http://tecdat.cn/?p=21317

最近咱们被客户要求撰写对于散布滞后非线性模型（DLNM）的钻研报告，包含一些图形和统计输入。

本文提供了运行散布滞后非线性模型的示例，同时形容了预测变量和后果之间的非线性和滞后效应，这种互相关系被定义为裸露 - 滞后 - 反馈关联

数据

数据集蕴含 1987-2000 年期间每日死亡率（CVD、呼吸道），天气（温度，相对湿度）和净化数据（PM10 和臭氧）。数据是由衰弱影响研究所资助的《国家发病率，死亡率和空气污染钻研》（NMMAPS）的一部分[Samet et al.,2000a,b]。

该钻研是对于随工夫变动的职业裸露与癌症之间的关系。该钻研包含 250 个危险集，每个危险集都有一个病例和一个对照，并与年龄相匹配。裸露数据以 15 岁至 65 岁之间的 5 岁年龄区间收集。

数据集药物蕴含模仿数据，来自一个假如的随机对照试验，对随工夫变动剂量的药物的影响。该钻研包含 200 名随机受试者，每人每天承受药物剂量，继续 28 天，每周都有变动。每隔 7 天报告一次。

DLNM 办法

在这里，我提供了一个简短的摘要来介绍概念和定义。

裸露 - 滞后 - 反馈关联

DLNM 的建模类用于形容关联，在该关联中，裸露和后果之间的依赖关系会在工夫上滞后。能够应用两个不同且互补的观点来形容此过程。咱们能够说，在工夫 t 处的裸露事件确定了在工夫 t + l 处的将来危险。应用后向视角，工夫 t 的危险由过来在工夫 t - l 经验的一系列危险确定。这里的 l 是滞后，示意裸露和测得的后果之间的滞后。

DLNM 统计模型

DLNM 类提供了一个概念和剖析框架，用于形容和预计裸露 - 滞后 - 反馈关联。DLNM 的统计倒退基于以下抉择：DLNM 类为形容和预计裸露 - 滞后 - 反馈关联提供了一个概念和剖析框架。DLNM 的统计倒退基于该抉择。

裸露 - 滞后 - 反馈关联的一个简略状况是，预测变量空间中的关系（即裸露 - 滞后关系）是线性的。能够通过 DLM 对这种类型的关系进行建模。在这种状况下，关联仅取决于滞后反馈函数，该函数模仿线性危险如何随滞后变动。滞后反馈函数的不同抉择（样条曲线，多项式，档次，阈值等）导致指定了不同的 DLM，并暗示了滞后反馈关系的代替假如。

DLNM 解释

DLNM 的后果能够通过应用 3 - D 绘图提供沿两个维度变动的关联，通过为每个滞后和预测变量的拟合值构建预测网格来解释。

第一是与特定裸露值相关联的滞后反馈曲线，定义为预测变量特定性关联。这被解释为与工夫 t 危险相干的工夫 t + l 的危险奉献序列。

第二是与特定滞后值相关联的裸露 - 反馈曲线，该特定滞后值定义为滞后特定关联。这被解释为与在工夫 t 处产生的裸露值相关联的在工夫 t + l 处的裸露 - 反馈关系。

第三个也是最重要的是与在思考的滞后期内经验的整个裸露历史相干的裸露反馈曲线，定义为总体累积关联。应用正向视角，这被解释为示意工夫 t 产生的给定裸露期间 [t，t+L] 期间经验的净危险的裸露反馈关系。

工夫序列之外的利用

散布滞后模型首先是在很久以前的计量经济工夫序列剖析中提出的 [Almon，1965]，而后在环境流行病学 Schwartz [2000] 的工夫序列数据中从新提出。DLNM 的扩大是由 Armstrong [2006]构想的。Gasparrini 等人对工夫序列数据的建模框架进行了从新评估。[2010]。乏味的是，曾经在不同的钻研畛域中提出了这种裸露 - 滞后 - 反馈关联的模型。个别的想法是通过特定函数加权过来的裸露，这些函数的参数由数据估算。在癌症流行病学 [Hauptmann 等，2000；Langholz 等，1999；Richardson，2009；Thomas，1983；Vacek，1997] 和药物流行病学 [Abrahamowicz 等] 中，阐明了相似于 DLM 的线性 - 裸露 - 反馈关系模型。

根本函数

指定规范裸露反馈和滞后反馈关系的根本函数，例如多项式，分层或阈值函数。例如，样条线由举荐的包样条线中蕴含的函数 ns（）和 bs（）指定。多项式是通过函数 poly（）取得的。这是一个简略向量的转换示例：

poly(1:5,degree=3)
1 2 3
[1,] 0.2 0.04 0.008
[2,] 0.4 0.16 0.064
[3,] 0.6 0.36 0.216
[4,] 0.8 0.64 0.512
[5,] 1.0 1.00 1.000
attr(,"degree")
[1] 3
attr(,"scale")
[1] 5
attr(,"intercept")
[1] FALSE
attr(,"class")
[1] "poly" "matrix"

第一个未命名的参数 x 指定要转换的向量，而参数度设置多项式的度。定义分层函数是通过 strata（）指定的。

strata(1:5,breaks=c(2,4))[,]
1 2
[1,] 0 0
[2,] 1 0
[3,] 1 0
[4,] 0 1
[5,] 0 1

后果是带有附加类别“层”的根底矩阵。转换是定义比照的虚构参数化。参数 break 定义了层的右凋谢区间的下边界。

阈值函数通过 thr（）指定。一个例子：

thr(1:5,thr.value=3,side="d")[,]
1 2
[1,] 2 0
[2,] 1 0
[3,] 0 0
[4,] 0 1
[5,] 0 2

后果是具备附加类别“thr”的根底矩阵。参数 thr.value 定义一个带有一个或两个阈值的向量，而 side 用于指定高（“h”，默认值），低（“l”）或双精度（“d”）阈值参数化。

点击题目查阅往期内容

【视频】R 语言中的散布滞后非线性模型（DLNM）与发病率，死亡率和空气污染示例

左右滑动查看更多

01

02

03

04

根本转换

此函数代表以 dlnm 为单位进行根本转换的次要函数，实用于指定裸露 - 反馈和滞后 - 反馈关系。它的作用是利用选定的转换并以实用于其余函数（例如 crossbasis（）和 crosspred（））的格局生成根本矩阵。以下示例复制了该局部中显示的多项式变换：

onebasis(1:5,fun="poly",degree=3)
b1 b2 b3
[1,] 0.2 0.04 0.008
[2,] 0.4 0.16 0.064
[3,] 0.6 0.36 0.216
[4,] 0.8 0.64 0.512
[5,] 1.0 1.00 1.000
attr(,"fun")
[1] "poly"
attr(,"degree")
[1] 3
attr(,"scale")
[1] 5
attr(,"intercept")
[1] FALSE
attr(,"class")
[1] "onebasis" "matrix"
attr(,"range")
[1] 1 5

后果是带有附加类“onebasis”的根底矩阵。同样，第一个未命名参数 x 指定要转换的向量，而第二个参数 fun 将字符转换定义为利用转换而调用的函数的名称。具体来说，根本矩阵包含 fun 和 range 属性，以及定义转换的被调用函数的参数。如前所述，onebasis（）还能够依据特定要求调用用户定义的函数。一个简略的例子：

> mylog <- function(x) log(x)
> onebasis(1:5,"mylog")
b1
[1,] 0.0000000
[2,] 0.6931472
[3,] 1.0986123
[4,] 1.3862944
[5,] 1.6094379
attr(,"fun")
[1] "mylog"
attr(,"range")
[1] 1 5
attr(,"class")
[1] "onebasis" "matrix"

穿插基

这是 dlnm 软件包中的次要函数。它在外部调用 onebasis（）来生成裸露 - 反馈和滞后 - 反馈关系的基矩阵，并通过非凡的张量积将它们组合起来，以创立穿插基，该穿插基在模型中同时指定了裸露 - 滞后 - 反馈关联性。它的第一个参数 x 的类定义如何解释数据。能够应用第二个变量 lag 批改滞后期。

作为一个简略的示例，我模仿了 2 - 5 个滞后期内 3 个对象的裸露历史矩阵：它们中的每一个都将传递给 onebasis（）来别离构建裸露 - 反馈和滞后 - 反馈关系的矩阵。仅用于工夫序列数据的附加参数组定义了被视为独自无关序列的察看组，例如在季节性剖析中可能有用。作为一个简略的示例，我模仿了 2 - 5 个滞后期内 3 个对象的裸露历史矩阵：它们中的每一个都将传递给 onebasis（）来别离构建裸露 - 反馈和滞后 - 反馈关系的矩阵。作为一个简略的示例，我模仿了 2 - 5 个滞后期内 3 个对象的裸露历史矩阵：

 > hist
lag2 lag3 lag4 lag5
sub1 1 3 5 6
sub2 7 8 9 4
sub3 10 2 11 12

而后，我利用穿插基参数化，将二次多项式作为裸露反馈函数，并将分层函数 2 - 3 和 4 - 5 定义为滞后反馈函数的分层函数：

lag=c(2,5),argvar=list(fun="poly",degree=2),
arglag=list(fun="strata",breaks=4))[,]
v1.l1 v1.l2 v2.l1 v2.l2
sub1 1.250000 0.9166667 0.4930556 0.4236111
sub2 2.333333 1.0833333 1.4583333 0.6736111
sub3 2.916667 1.9166667 2.5625000 1.8402778

该函数返回“crossbasis”类的矩阵对象。它首先应用 argvar 和 arglag 列表中的参数调用 onebasis（），以建设裸露反馈空间和滞后反馈空间的矩阵根底。在另一个示例中，我将 crossbasis（）利用于数据集中的变量 temp，该数据集示意 1987-2000 年期间日平均温度序列：

 > summary(cb)
CROSSBASIS FUNCTIONS
observations: 5114
range: -26.66667 to 33.33333
lag period: 0 30
total df: 10
BASIS FOR VAR:
fun: thr
thr.value: 10 20
side: d
intercept: FALSE
BASIS FOR LAG:
fun: ns
knots: 1 4 12
intercept: TRUE
Boundary.knots: 0 30

此处，将裸露反馈建模为阈值为 10 和 20 的双阈值函数。滞后工夫设置为 0 到 30。滞后反馈函数留给默认的天然三次样条（fun =“ns”），其滞后值为 1、4 和 12。

预测

crossbasis（）生成的穿插基矩阵须要蕴含在回归模型公式中能力拟合模型。关联通过函数 crosspred（）进行汇总，该函数针对默认值或用户间接抉择的预测值和滞后值的组合的网格进行预测。例如，我应用创立的穿插基矩阵 cb，应用数据集工夫序列数据来钻研温度与心血管疾病死亡率之间的关联。首先，我将一个简略的线性模型与模型公式中蕴含的穿插基矩阵拟合。而后，我通过应用 cross-basis 和回归模型对象作为前两个参数调用 crosspred（）来取得预测：

crosspred(cb,model,at=-20:30)

后果是“crosspred”类的列表对象，其中的存储预测和无关模型的其余信息，例如系数和与穿插基参数无关的关联（协）方差矩阵的一部分。能够为特定的预测器 - 滞后组合抉择预测的网格。例如，我提取温度为 -10°C 且滞后 5 的预测和置信区间，而后提取 25°C 的整体累积预测：

 > pred$allfit["25"]
25
1.108262

第一个结果表明，在给定的一天中，-20°C 的温度会在五天后导致 0.95 例心血管死亡的减少，或者在给定的一天中，温度为 - 6 摄氏度时，心血管死亡的数目减少 0.95。其余类型的预测能够通过 crosspred（）取得。特地是，如果模型链接等于 log 或 logit，则将主动返回取幂的预测。如果参数 cum 设置为 TRUE，则是累积预测的矩阵 cum。

crosspred（）的另一种用法是预测特定的裸露历史记录集的影响。这能够通过输出裸露历史矩阵作为参数来实现。例如，咱们能够从拟合模型中预测出，在过来 10 天裸露于 30°C 和在滞后期的其余工夫裸露于 22°C 之后，心血管死亡的总体累积减少：如果参数 cum 设置为 TRUE，则包含增量累积预测的矩阵 cum，并将其存储在组件 cum 中。crosspred（）的另一种用法是预测特定的裸露历史记录集的影响。这能够通过输出裸露历史矩阵作为参数来实现。例如，咱们能够从拟合模型中预测出，在过来 10 天裸露于 30°C 和在滞后期的其余工夫裸露于 22°C 之后，心血管死亡的总体累积减少：

 > crosspred(cb,model,at=histpred)$allfit
1
5.934992

dlnm 软件包的次要长处之一是，用户能够应用规范回归函数执行 DLNM，只需在模型公式中包含穿插基矩阵即可。函数 crosspred（）主动解决来自回归函数 lm（）和 glm（），gam（）（软件包 mgcv），coxph（）的模型。

降维

DLNM 的拟合度能够升高到预测变量或滞后的一个维度，仅以诸如总累积裸露反馈表白。该计算通过函数 crossreduce（）进行，该函数的工作原理与 crosspred（）十分类似。前两个自变量 base 和 model 指定穿插基矩阵和须要对其执行计算的模型对象。缩小的类型由类型定义，带有选项“overall”-“lag”-“var”，用于汇总总体累积裸露反馈，滞后特异性裸露反馈或预测变量特异性滞后反馈。

绘图

一维或二维关联的解释通过图形示意来辅助。通过办法函数 plot（），lines（）和 points（）为类“crosspred”和“crossreduce”提供高级和低级绘图性能。例如，我应用对象 pred 中的预测。plot（）办法能够通过参数 ptype 为“crosspred”对象生成不同类型的图。具体来说，它会生成整个二维裸露 - 滞后 - 反馈关联的图形。二维关联能够绘制为 3 - D 或等高线图，例如：

> plot(pred,ptype="3d",main="3D plot"

能够通过抉择不同的 ptype 取得定义的关联的摘要。

> plot(pred,"overall"

在这种状况下，办法函数 plot（）在外部调用函数 plot.default（），如下面的示例所示，能够将其余特定参数增加到函数调用中。通过设置 ptype =“slices”，能够将滞后特异性和预测因子特异性关联别离绘制为裸露 - 反馈和滞后 - 反馈关系，因为它们是在 3 - D 曲面中沿特定维度切割的切片。例如：

> plot(pred,"slices",lag=5

这两个图别离代表了滞后 5 的裸露反馈和特定于 25°C 温度的滞后反馈。参数 lag 和 var 指定必须别离绘制 lag 和特定于预测变量的关联的值。

点击文末 “浏览原文”

获取全文残缺代码数据资料。

本文选自《R 语言散布滞后非线性模型（DLNM）钻研发病率，死亡率和空气污染示例》。

点击题目查阅往期内容

R 语言散布滞后线性和非线性模型（DLM 和 DLNM）建模
散布滞后线性和非线性模型（DLNM）剖析空气污染（臭氧）、温度对死亡率工夫序列数据的影响
R 语言中的散布滞后非线性模型 DLNM 与发病率和空气污染示例
【视频】R 语言中的散布滞后非线性模型（DLNM）与发病率，死亡率和空气污染示例
R 语言散布滞后线性和非线性模型（DLNM）剖析空气污染（臭氧）、温度对死亡率工夫序列数据的影响
R 语言散布滞后线性和非线性模型（DLMs 和 DLNMs）剖析工夫序列数据
R 语言散布滞后非线性模型（DLNM）空气污染钻研温度对死亡率影响建模利用 R 语言散布滞后非线性模型（DLNM）钻研发病率，死亡率和空气污染示例
R 语言散布滞后线性和非线性模型（DLM 和 DLNM）建模
R 语言狭义相加模型 (GAMs)剖析预测 CO2 工夫序列数据
Python | ARIMA 工夫序列模型预测航空公司的乘客数量
R 语言中生存分析模型的工夫依赖性 ROC 曲线可视化
R 语言 ARIMA，SARIMA 预测道路交通流量工夫序列剖析: 季节性、周期性
ARIMA 模型预测 CO2 浓度工夫序列 -python 实现
R 语言基于递归神经网络 RNN 的温度工夫序列预测
R 语言用多元 ARMA,GARCH ,EWMA, ETS, 随机稳定率 SV 模型对金融工夫序列数据建模
R 语言神经网络模型预测车辆数量工夫序列
卡尔曼滤波器：用 R 语言中的 KFAS 建模工夫序列
在 Python 中应用 LSTM 和 PyTorch 进行工夫序列预测
R 语言从经济工夫序列中用 HP 滤波器，小波滤波和教训模态合成等提取周期性成分剖析
应用 PYTHON 中 KERAS 的 LSTM 递归神经网络进行工夫序列预测
Python 中的 ARIMA 模型、SARIMA 模型和 SARIMAX 模型对工夫序列预测
R 语言 k -Shape 工夫序列聚类办法对股票价格工夫序列聚类
R 语言多元 Copula GARCH 模型工夫序列预测