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

原文出处：拓端数据部落公众号

前言

本文阐明了 R 语言中实现散布滞后线性和非线性模型（DLM 和 DLNM）的建模。首先，本文形容了除工夫序列数据之外的 DLM / DLNM 的一般化办法，在 Gasparrini [2014] 中有更具体的形容。本文中蕴含的后果并不代表迷信发现，而仅出于阐明目标进行报告。

数据

次要通过两个示例来阐明软件的利用，应用药物数据作为数据对象。数据集别离蕴含一项对于药物的假如试验和嵌套病例对照钻研的模仿数据，两者均包含随工夫变动的裸露量度。

让咱们看一下数据框的前 2 个察看样本：

> head(data, 2)
id out sex day1 day8. day15. day22.
1 1 46 M 0 0 40 37
2 2 50 F 0 47 55 0

数据集蕴含来自一项试验的数据，记录了 200 名随机受试者，每名受试者随机承受周围中两周的药物剂量，每天的剂量每周变动。每周 7 天距离报告一次裸露程度。数据集还蕴含无关在第 28 天测量的后果和受试者性别的信息。嵌套的第二个数据包含针对 300 个癌症病例和 300 个按年龄匹配的对照的每个记录。前 2 个察看后果是：

> head(nested)
id case age riskset exp15 exp20 exp25 exp30 exp35 exp40 exp45 exp50 exp55
1 1 1 81 240 5 84 34 45 128 81 14 52 11
2 2 1 69 129 11 8 25 6 8 12 19 60 16

exp60
1 16
2 10

变量病例定义病例 / 对照状态，而其余变量报告受试者的年龄和他 / 她所属的危险。随工夫变动的职业裸露档案存储在变量 exp15–exp60 中，对应于 15 至 19 岁，20 至 24 岁等最高 65 岁的均匀年裸露量。

裸露历史矩阵

扩大的 DLNM 框架与规范 DLNM 框架之间的次要区别是裸露历史矩阵的定义，即对 n 个观测值的滞后 ` 经验的一系列裸露。依据钻研设计和随工夫变动的裸露信息，须要以不同的形式将这个 n×（L −’0 + 1）矩阵组合在一起。

在第一个示例中，我为数据框药物中的试验数据建设了裸露历史记录矩阵。

每个受试者的接触曲线用于重建接触历史矩阵。在这种状况下，滞后 0 的裸露量对应于对所有受试者测量终局的第 28 天的裸露量。其余的裸露历史记录可追溯到滞后 27，对应于第一天的裸露。代码，用于将按周存储的裸露材料扩大为每日裸露历史记录的矩阵：每个受试者的接触曲线用于重建接触历史矩阵。

> drug\[1:3,1:14\]
lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13
1 37 37 37 37 37 37 37 40 40 40 40 40 40 40
2 0 0 0 0 0 0 0 55 55 55 55 55 55 55

下面针对前三个主题报告了滞后 0 -13 的接触历史。前七个滞后（0–6）对应于上周的裸露，而滞后 7–13 对应于第三周，依此类推。在第二个示例中，我应用以 5 年为距离的裸露量分布图来嵌套数据框的裸露量历史矩阵。这些数据被扩大为滞后 3–40 的裸露历史矩阵，滞后单位等于一年。然而，在这种状况下，因为每个对象在不同的年龄进行采样，因而计算更加简单。具体地，从受检者的年龄开始沿着裸露曲线向后计算裸露历史。此步骤须要一些额定的计算和数据处理。能够得出给定工夫的裸露曲线的裸露历史，

> nest <- t(apply(nested, 1, function(sub) exphi(repc(0,0,0,sub5:14\]),


> nest\[1:3,1:11\]
lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13
1 0 0 0 0 0 0 0 0 0 0 0
2 0 10 10 10 10 10 16 16 16 16 16

下面针对前三个主题报告了滞后 0–10 的裸露历史。假如第一个对象在 81 岁时进行采样，则经验了在滞后 0 处介于 80 和 81 之间，在滞后 1 处介于 79 和 80 之间的裸露，依此类推。因为他 / 她的上一次裸露年龄为 65 岁，因而将滞后 10 的裸露历史记录设置为 0。在 69 岁时进行采样的第二个对象的滞后 3 的裸露历史记录设置为 0，对应于裸露事件在 66。

这些接触历史与之前显示的接触详情和年龄统一。在这种状况下，应用雷同的裸露情况，在每个受试者奉献不同危险集时计算每个受试者的屡次裸露历史。通常，此矩阵的计算取决于钻研设计，裸露信息，滞后单位和所需的近似程度。

工夫序列以外的利用

一个简略的 DLM

在第一个示例中，我将 dlnm 利用于数据集药物，剖析了药物日剂量与未指定衰弱后果之间的工夫依赖性。第一步是函数的定义：

crossbasis(drug, lag=27, argvar=list("lin")

后果存储在对象 cbdrug 中，即具备非凡属性的已转换变量的矩阵。参数 argvar 和 arglag 别离定义了裸露反馈和滞后反馈函数，此处抉择它们为简略线性函数和三次样条。通过函数 summary 取得摘要：

CROSSBASIS FUNCTIONS
observations: 200
range: 0 to 100
lag period: 0 27
total df: 4
BASIS FOR VAR:
fun: ns
knots: 9 18

cbdrug 能够蕴含在回归模型的公式中，在这种状况下，该模型是假如高斯分布，管制性别影响的简略线性模型。通过函数 crosspred（）预测来解释预计的滞后关联：

 crosspred(cbdrug, mdrug, at=0:20*5)

成果摘要保留在“crosspred”类的对象 pdrug 中

allfit alllow allhigh
30.29 20.12 40.46

下面的代码提取了与 50 次裸露相干的总体累积效应的估算值，能够进行解释：在 28 天滞后工夫内继续一直地裸露于 50 次之后的总体后果减少。还包含 95%的置信区间。

能够生成图：

> plot(drug, zlab="Effect", xlab="Dose, ylab="Lag (days")

代码的第一行产生图 1 中的图形，显示成果在剂量和滞后值的范畴内如何变动。该图表明，在摄入后的头几天，该剂量的药物作用显著，而后在 15-20 天后趋于隐没。从横截面来看，图别离显示了裸露 60 的滞后反馈曲线和滞后 10 的裸露 - 反馈曲线。图中的滞后反馈曲线表明了效应的指数衰减。

更为简单的 DLNM 

在第二个示例中，我应用嵌套的数据集来评估长期裸露于职业病中如何影响癌症产生的危险。剖析步骤与阐明的步骤雷同。最后的假如是，过来三年中继续的裸露（对应于滞后 0–2）不会影响产生癌症的危险。

抉择的基函数是用于预测变量的二次样条和三次样条。通过 clogit（）执行条件逻辑回归。而后预测成果摘要。代码是：

> library(survival)
clogit(case~cbnest+strata(riskset, nested)

图 1 中显示的雷同类型的图可通过以下形式取得：

> plot(est, zlab="OR", xlab="Exposure, ylab=Lag (years")

图中的 3 - D 图再次被解释为职业裸露与癌症危险之间的关联。在此示例中，滞后时间段以年为单位示意。该图表明危险的初始减少（以比值比（OR）掂量），而后升高。从横截面来看，图中预计的滞后反馈曲线显示了裸露后 10 至 15 年的危险峰值，只管置信区间十分宽泛，但危险在裸露后 30 年回到根底程度。

扩大预测

之前取得的预测后果是在间接指定的曝露和滞后值的网格上计算的。

咱们也能够计算新的成果摘要，在给定裸露曲线的状况下生成裸露历史矩阵。例如，咱们能够应用嵌套病例对照剖析来计算，假如受试者裸露于裸露 10 年达 5 年，而后未裸露于 5 年，再裸露于 13 年达 10 年的总体累计 OR。从此裸露量配置中，咱们能够计算出裸露工夫完结时的裸露历史，并预测。

> hist
lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15 lag16
20 13 13 13 13 13 13 13 0 0 0 0 0 10 10

产生工夫 3 到 40 的滞后工夫的裸露历史。通过自变量工夫设置特定工夫，在这种状况下，该工夫对应于裸露工夫的完结工夫（以指数示意）。包含最近的 21 次裸露至 0，以实现长达 40 年的裸露历史。当初，咱们能够应用 hist 作为 crosspred（）的参数来预测总体累积成果。留神，滞后周期必须与预计中应用的统一。

> with(pnestt, allRfit,aRRlow,allRRigh)
20 20 20
3.5031.2409.900

与在整个滞后期间没有裸露的受试者相比，预计的 OR 为 3.5（95%CI：1.2–9.9）。能够应用雷同的办法来获取特定裸露量散布随工夫的动静预测。这个思维是基于假设的裸露 - 滞后 - 反馈关联，在给定随工夫变动的裸露历史的状况下，及时地动静预测危险。实际上，对于每个给定的工夫，随着特定的裸露事件波及不同的滞后工夫，裸露历史会发生变化。举例来说，我展现了如何应用试验数据分析来估算特定药物处方后的动静预测成果。

假如某位患者承受 10 剂量的医治，继续 2 周，而后他 / 她减少至 50，继续 1 周，而后停药 1 周，而后以 20 的剂量从新开始医治 2 周。首先，我创立每日裸露材料：

> expdrug <- rep(c(10,50,0,20),c(2,1,1,2)*7)

当初能够沿裸露曲线程序来创立所有工夫点的裸露历史矩阵：

> nhist <- exphi(expdr, lag=27)

当初能够在 crosspred（）中应用此矩阵来获取动静预测。

当初能够应用该对象绘制动静预测：

> plot(drug,"overall", ylab="Effect xlab="Time (days", ylim=c(-10,27)
> axis(2, at=-1:5*5)
> par(new=TRUE)
> axis(4, at=0:6*10, cex.axis=0.8)

在图中绘制了整体累积关联。此图显示了与下面具体介绍的药物处方相干的基线后果的变动。正如预期的那样，成果会随剂量动态变化，但会呈现滞后。

利用改良函数

对于第一个示例，咱们能够批改先前剖析。图 2 倡议在高裸露量下可能会削弱成果。这个事实和裸露散布的偏斜度能够通过对数变换来解决。首先，让咱们定义一个新的函数：

log <- function(x) log(x+1)

当初能够建模裸露 - 反馈曲线：

nest2 <- crossbasis(est, lag=c(3,40), argvar=list("log"),

CROSSBASIS FUNCTIONS

observations: 600
range: 0 to 1064
lag period: 3 40
total df: 3
BASIS FOR VAR:
fun: mylog
BASIS FOR LAG:
fun: ns
knots: 10 30
intercept: FALSE
Boundary.knots: 3 40

替换新创建的对象：

能够将图中显示的后果与最后显示的后果进行比拟。该比拟表明对数的假如使精度大大提高。

对图的查看表明，滞后反馈曲线遵循指数衰减轨迹。利用衰减函数而不是三次样条曲线可能是正当的。衰减函数能够定义为：

decay <- function(x,scale=5) 
basis <- exp(-x/scale)
attributes(basis)$scale <- scale

参数（默认值为 5）用于管制衰减水平。同样，咱们能够应用此新函数来取得变换：

> cbdrug2 <- crossbasis(Qdrug, lag=27,
arglag=list(fun="fdecay")

CROSSBASIS FUNCTIONS
observations: 200
range: 0 to 100
lag period: 0 27
total df: 1
BASIS FOR VAR:
fun: lin
intercept: FALSE
BASIS FOR LAG:
fun: fdecay
scale: 6

同样，能够重复使用计算步骤以执行批改后的剖析：

> lines(drug, var=60, lty=2)
> lines(drug, lag=10, lty=2)

后果报告在图中。与之前的后果进行比拟（以虚线示意）显示了精度的显着进步。

回归剖析的通用工具

软件包 dlnm 中的性能也能够用作回归剖析的通用工具。第一个示例演示了如何应用带有回归函数 lm（）的回归样条来评估 30-39 岁的女性样本中均匀身高和体重之间的关系。

> library(splines)
> oneheight <- onebasis(women$height, "ns" df=5)
> mwomen <- lm(weight ~ oneheight data=women)

应用一个简略的代码来获取预测和绘图：

> with(pwomen, cbind(fit, low, high)\["70",)
allfit alllow allhigh
18.92287 18.46545 19.38030

能够简略地查看带有置信区间的预计关联，绘制关联。

第二个示例应用惩办样条对平滑关联进行剖析。

> library(mgcv)

> b2 <- gam(y ~ s(x0,bs="cr") + s(x1,bs="cr") + s(x2,bs="cr") + s(x3,bs="cr"),
family=poisson, data=datmethod="REML")
> plot(b2, select=3)

该代码应用通过函数 s（）的回归样条，对带有多个变量的模仿数据执行 GAM 预计平滑关系。也能够应用 dlnm 取得预测和绘图，其中：

allRRfit allRRlow allRRhigh
1.3405415 0.8309798 2.1625694
> plot(gam, ylim=c(0,3)col=2)

参考文献 

A. Gasparrini. Modeling exposure-lag-response associations with distributed lag non-linear models. Statistics in Medicine, 33(5):881–899, 2014.

最受欢迎的见解

1. 用 SPSS 预计 HLM 档次线性模型模型

2.R 语言线性判别分析（LDA），二次判别分析（QDA）和正则判别分析（RDA）

3. 基于 R 语言的 lmer 混合线性回归模型

4.R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真剖析

5. 在 r 语言中应用 GAM（狭义相加模型）进行电力负荷工夫序列剖析

6. 应用 SAS，Stata，HLM，R，SPSS 和 Mplus 的分层线性模型 HLM

7.R 语言中的岭回归、套索回归、主成分回归：线性模型抉择和正则化

8.R 语言用线性回归模型预测空气质量臭氧数据

9.R 语言分层线性模型案例