原文链接：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 372 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 exp551 1 1 81 240 5 84 34 45 128 81 14 52 112 2 1 69 129 11 8 25 6 8 12 19 60 16exp601 162 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 lag131 37 37 37 37 37 37 37 40 40 40 40 40 40 402 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 lag131 0 0 0 0 0 0 0 0 0 0 02 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 FUNCTIONSobservations: 200range: 0 to 100lag period: 0 27total df: 4BASIS FOR VAR:fun: nsknots: 9 18

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

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

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

allfit alllow allhigh30.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。从此裸露量配置中，咱们能够计算出裸露工夫完结时的裸露历史，并预测。

> histlag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15 lag1620 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 203.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 FUNCTIONSobservations: 600range: 0 to 1064lag period: 3 40total df: 3BASIS FOR VAR:fun: mylogBASIS FOR LAG:fun: nsknots: 10 30intercept: FALSEBoundary.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 FUNCTIONSobservations: 200range: 0 to 100lag period: 0 27total df: 1BASIS FOR VAR:fun: linintercept: FALSEBASIS FOR LAG:fun: fdecayscale: 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 allhigh18.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 allRRhigh1.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语言分层线性模型案例