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

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

前言

本文說明了R語言中實(shí)現(xiàn)分布滯后線性和非線性模型（DLM和DLNM）的建模。首先，本文描述了除時(shí)間序列數(shù)據(jù)之外的DLM / DLNM的一般化方法，在Gasparrini [2014]中有更詳細(xì)的描述。本文中包含的結(jié)果并不代表科學(xué)發(fā)現(xiàn)，而僅出于說明目的進(jìn)行報(bào)告。

數(shù)據(jù)

主要通過兩個(gè)示例來說明軟件的應(yīng)用，使用藥物數(shù)據(jù)作為數(shù)據(jù)對(duì)象。數(shù)據(jù)集分別包含一項(xiàng)關(guān)于藥物的假設(shè)試驗(yàn)和嵌套病例對(duì)照研究的模擬數(shù)據(jù)，兩者均包括隨時(shí)間變化的暴露量度。

讓我們看一下數(shù)據(jù)框的前2個(gè)觀察樣本：

1. > head(data, 2)

2. id out sex day1 day8. day15. day22.

3. 1 1 46 M 0 0 40 37

4. 2 2 50 F 0 47 55 0

數(shù)據(jù)集包含來自一項(xiàng)試驗(yàn)的數(shù)據(jù)，記錄了200名隨機(jī)受試者，每名受試者隨機(jī)接受四周中兩周的藥物劑量，每天的劑量每周變化。每周7天間隔報(bào)告一次暴露水平。數(shù)據(jù)集還包含有關(guān)在第28天測(cè)量的結(jié)果和受試者性別的信息。嵌套的第二個(gè)數(shù)據(jù)包括針對(duì)300個(gè)癌癥病例和300個(gè)按年齡匹配的對(duì)照的每個(gè)記錄。前2個(gè)觀察結(jié)果是：

1. > head(nested)

2. id case age riskset exp15 exp20 exp25 exp30 exp35 exp40 exp45 exp50 exp55

3. 1 1 1 81 240 5 84 34 45 128 81 14 52 11

4. 2 2 1 69 129 11 8 25 6 8 12 19 60 16

5. 


6. exp60

7. 1 16

8. 2 10

變量病例定義病例/對(duì)照狀態(tài)，而其他變量報(bào)告受試者的年齡和他/她所屬的風(fēng)險(xiǎn)。隨時(shí)間變化的職業(yè)暴露檔案存儲(chǔ)在變量exp15–exp60中，對(duì)應(yīng)于15至19歲，20至24歲等最高65歲的平均年暴露量。

暴露歷史矩陣

擴(kuò)展的DLNM框架與標(biāo)準(zhǔn)DLNM框架之間的主要區(qū)別是暴露歷史矩陣的定義，即對(duì)n個(gè)觀測(cè)值的滯后`經(jīng)歷的一系列暴露。根據(jù)研究設(shè)計(jì)和隨時(shí)間變化的暴露信息，需要以不同的方式將這個(gè)n×（L ?'0 + 1）矩陣組合在一起。

在第一個(gè)示例中，我為數(shù)據(jù)框藥物中的試驗(yàn)數(shù)據(jù)建立了暴露歷史記錄矩陣。

每個(gè)受試者的接觸曲線用于重建接觸歷史矩陣。在這種情況下，滯后0的暴露量對(duì)應(yīng)于對(duì)所有受試者測(cè)量結(jié)局的第28天的暴露量。其余的暴露歷史記錄可追溯到滯后27，對(duì)應(yīng)于第一天的暴露。代碼，用于將按周存儲(chǔ)的暴露資料擴(kuò)展為每日暴露歷史記錄的矩陣：每個(gè)受試者的接觸曲線用于重建接觸歷史矩陣。

1. 


2. > drug[1:3,1:14]

3. lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13

4. 1 37 37 37 37 37 37 37 40 40 40 40 40 40 40

5. 2 0 0 0 0 0 0 0 55 55 55 55 55 55 55

上面針對(duì)前三個(gè)主題報(bào)告了滯后0-13的接觸歷史。前七個(gè)滯后（0–6）對(duì)應(yīng)于上周的暴露，而滯后7–13對(duì)應(yīng)于第三周，依此類推。在第二個(gè)示例中，我使用以5年為間隔的暴露量分布圖來嵌套數(shù)據(jù)框的暴露量歷史矩陣。這些數(shù)據(jù)被擴(kuò)展為滯后3–40的暴露歷史矩陣，滯后單位等于一年。但是，在這種情況下，由于每個(gè)對(duì)象在不同的年齡進(jìn)行采樣，因此計(jì)算更加復(fù)雜。具體地，從受檢者的年齡開始沿著暴露曲線向后計(jì)算暴露歷史。此步驟需要一些額外的計(jì)算和數(shù)據(jù)處理?？梢缘贸鼋o定時(shí)間的暴露曲線的暴露歷史，

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

2. 


3. 


4. > nest[1:3,1:11]

5. lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13

6. 1 0 0 0 0 0 0 0 0 0 0 0

7. 2 0 10 10 10 10 10 16 16 16 16 16

上面針對(duì)前三個(gè)主題報(bào)告了滯后0–10的暴露歷史。假設(shè)第一個(gè)對(duì)象在81歲時(shí)進(jìn)行采樣，則經(jīng)歷了在滯后0處介于80和81之間，在滯后1處介于79和80之間的暴露，依此類推。由于他/她的上一次暴露年齡為65歲，因此將滯后10的暴露歷史記錄設(shè)置為0。在69歲時(shí)進(jìn)行采樣的第二個(gè)對(duì)象的滯后3的暴露歷史記錄設(shè)置為0，對(duì)應(yīng)于暴露事件在66。

這些接觸歷史與之前顯示的接觸概況和年齡一致。在這種情況下，使用相同的暴露狀況，在每個(gè)受試者貢獻(xiàn)不同風(fēng)險(xiǎn)集時(shí)計(jì)算每個(gè)受試者的多次暴露歷史。通常，此矩陣的計(jì)算取決于研究設(shè)計(jì)，暴露信息，滯后單位和所需的近似水平。

時(shí)間序列以外的應(yīng)用

一個(gè)簡(jiǎn)單的DLM

在第一個(gè)示例中，我將dlnm應(yīng)用于數(shù)據(jù)集藥物，分析了藥物日劑量與未指定健康結(jié)果之間的時(shí)間依賴性。第一步是函數(shù)的定義：

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

結(jié)果存儲(chǔ)在對(duì)象cbdrug中，即具有特殊屬性的已轉(zhuǎn)換變量的矩陣。參數(shù)argvar和arglag分別定義了暴露反應(yīng)和滯后反應(yīng)函數(shù)，此處選擇它們?yōu)楹?jiǎn)單線性函數(shù)和三次樣條。通過函數(shù)summary獲得摘要：

1. CROSSBASIS FUNCTIONS

2. observations: 200

3. range: 0 to 100

4. lag period: 0 27

5. total df: 4

6. BASIS FOR VAR:

7. fun: ns

8. knots: 9 18

cbdrug可以包含在回歸模型的公式中，在這種情況下，該模型是假設(shè)高斯分布，控制性別影響的簡(jiǎn)單線性模型。通過函數(shù)crosspred（）預(yù)測(cè)來解釋估計(jì)的滯后關(guān)聯(lián)：

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

效果摘要保存在“ crosspred”類的對(duì)象pdrug中

1. allfit alllow allhigh

2. 30.29 20.12 40.46

上面的代碼提取了與50次暴露相關(guān)的總體累積效應(yīng)的估算值，可以進(jìn)行解釋：在28天滯后時(shí)間內(nèi)持續(xù)不斷地暴露于50次之后的總體結(jié)果增加。還包括95％的置信區(qū)間。

可以生成圖：

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

2. 


代碼的第一行產(chǎn)生圖1中的圖形，顯示效果在劑量和滯后值的范圍內(nèi)如何變化。該圖表明，在攝入后的頭幾天，該劑量的藥物作用明顯，然后在15-20天后趨于消失。從橫截面來看，圖分別顯示了暴露60的滯后反應(yīng)曲線和滯后10的暴露-反應(yīng)曲線。圖中的滯后反應(yīng)曲線表明了效應(yīng)的指數(shù)衰減。

更為復(fù)雜的DLNM?

在第二個(gè)示例中，我使用嵌套的數(shù)據(jù)集來評(píng)估長(zhǎng)期暴露于職業(yè)病中如何影響癌癥發(fā)生的風(fēng)險(xiǎn)。分析步驟與說明的步驟相同。最初的假設(shè)是，過去三年中持續(xù)的暴露（對(duì)應(yīng)于滯后0–2）不會(huì)影響發(fā)生癌癥的風(fēng)險(xiǎn)。

選擇的基函數(shù)是用于預(yù)測(cè)變量的二次樣條和三次樣條。通過clogit（）執(zhí)行條件邏輯回歸。然后預(yù)測(cè)效果摘要。代碼是：

1. > library(survival)

2. clogit(case~cbnest+strata(riskset, nested)

圖1中顯示的相同類型的圖可通過以下方式獲得：

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

圖中的3-D圖再次被解釋為職業(yè)暴露與癌癥風(fēng)險(xiǎn)之間的關(guān)聯(lián)。在此示例中，滯后時(shí)間段以年為單位表示。該圖表明風(fēng)險(xiǎn)的初始增加（以比值比（OR）衡量），然后降低。從橫截面來看，圖中估計(jì)的滯后反應(yīng)曲線顯示了暴露后10至15年的風(fēng)險(xiǎn)峰值，盡管置信區(qū)間非常寬泛，但風(fēng)險(xiǎn)在暴露后30年回到基礎(chǔ)水平。

擴(kuò)展預(yù)測(cè)

之前獲得的預(yù)測(cè)結(jié)果是在直接指定的曝露和滯后值的網(wǎng)格上計(jì)算的。

我們也可以計(jì)算新的效果摘要，在給定暴露曲線的情況下生成暴露歷史矩陣。例如，我們可以使用嵌套病例對(duì)照分析來計(jì)算，假設(shè)受試者暴露于暴露10年達(dá)5年，然后未暴露于5年，再暴露于13年達(dá)10年的總體累計(jì)OR。從此暴露量配置中，我們可以計(jì)算出暴露時(shí)間結(jié)束時(shí)的暴露歷史，并預(yù)測(cè)。

1. 


2. > hist

3. lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15 lag16

4. 20 13 13 13 13 13 13 13 0 0 0 0 0 10 10

產(chǎn)生時(shí)間3到40的滯后時(shí)間的暴露歷史。通過自變量時(shí)間設(shè)置特定時(shí)間，在這種情況下，該時(shí)間對(duì)應(yīng)于暴露時(shí)間的結(jié)束時(shí)間（以指數(shù)表示）。包括最近的21次暴露至0，以完成長(zhǎng)達(dá)40年的暴露歷史?，F(xiàn)在，我們可以使用hist作為crosspred（）的參數(shù)來預(yù)測(cè)總體累積效果。注意，滯后周期必須與估計(jì)中使用的一致。

1. > with(pnestt, allRfit,aRRlow,allRRigh)

2. 20 20 20

3. 3.5031.2409.900

與在整個(gè)滯后期間沒有暴露的受試者相比，估計(jì)的OR為3.5（95％CI：1.2–9.9）?？梢允褂孟嗤姆椒▉慝@取特定暴露量分布隨時(shí)間的動(dòng)態(tài)預(yù)測(cè)。這個(gè)思想是基于假定的暴露-滯后-反應(yīng)關(guān)聯(lián)，在給定隨時(shí)間變化的暴露歷史的情況下，及時(shí)地動(dòng)態(tài)預(yù)測(cè)風(fēng)險(xiǎn)。實(shí)際上，對(duì)于每個(gè)給定的時(shí)間，隨著特定的暴露事件涉及不同的滯后時(shí)間，暴露歷史會(huì)發(fā)生變化。舉例來說，我展示了如何使用試驗(yàn)數(shù)據(jù)分析來估算特定藥物處方后的動(dòng)態(tài)預(yù)測(cè)效果。

假設(shè)某位患者接受10劑量的治療，持續(xù)2周，然后他/她增加至50，持續(xù)1周，然后停藥1周，然后以20的劑量重新開始治療2周。首先，我創(chuàng)建每日暴露資料：

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

現(xiàn)在可以沿暴露曲線順序來創(chuàng)建所有時(shí)間點(diǎn)的暴露歷史矩陣：

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

現(xiàn)在可以在crosspred（）中使用此矩陣來獲取動(dòng)態(tài)預(yù)測(cè)。

現(xiàn)在可以使用該對(duì)象繪制動(dòng)態(tài)預(yù)測(cè)：

1. > plot(drug,"overall", ylab="Effect xlab="Time (days", ylim=c(-10,27)

2. > axis(2, at=-1:5*5)

3. > par(new=TRUE)

4. > axis(4, at=0:6*10, cex.axis=0.8)

在圖中繪制了整體累積關(guān)聯(lián)。此圖顯示了與上面詳細(xì)介紹的藥物處方相關(guān)的基線結(jié)果的變化。正如預(yù)期的那樣，效果會(huì)隨劑量動(dòng)態(tài)變化，但會(huì)出現(xiàn)滯后。

應(yīng)用改進(jìn)函數(shù)

對(duì)于第一個(gè)示例，我們可以修改先前分析。圖2建議在高暴露量下可能會(huì)減弱效果。這個(gè)事實(shí)和暴露分布的偏斜度可以通過對(duì)數(shù)變換來解決。首先，讓我們定義一個(gè)新的函數(shù)：

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

現(xiàn)在可以建模暴露-反應(yīng)曲線：

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

2. 


3. CROSSBASIS FUNCTIONS

4. 


5. observations: 600

6. range: 0 to 1064

7. lag period: 3 40

8. total df: 3

9. BASIS FOR VAR:

10. fun: mylog

11. BASIS FOR LAG:

12. fun: ns

13. knots: 10 30

14. intercept: FALSE

15. Boundary.knots: 3 40

替換新創(chuàng)建的對(duì)象：

可以將圖中顯示的結(jié)果與最初顯示的結(jié)果進(jìn)行比較。該比較表明對(duì)數(shù)的假設(shè)使精度大大提高。

對(duì)圖的檢查表明，滯后反應(yīng)曲線遵循指數(shù)衰減軌跡。應(yīng)用衰減函數(shù)而不是三次樣條曲線可能是合理的。衰減函數(shù)可以定義為：

1. decay <- function(x,scale=5)

2. basis <- exp(-x/scale)

3. attributes(basis)$scale <- scale

4. 


參數(shù)（默認(rèn)值為5）用于控制衰減程度。同樣，我們可以使用此新函數(shù)來獲得變換：

1. > cbdrug2 <- crossbasis(Qdrug, lag=27,

2. arglag=list(fun="fdecay")

3. 


4. CROSSBASIS FUNCTIONS

5. observations: 200

6. range: 0 to 100

7. lag period: 0 27

8. total df: 1

9. BASIS FOR VAR:

10. fun: lin

11. intercept: FALSE

12. BASIS FOR LAG:

13. fun: fdecay

14. scale: 6

同樣，可以重復(fù)使用計(jì)算步驟以執(zhí)行修改后的分析：

1. 


2. > lines(drug, var=60, lty=2)

3. > lines(drug, lag=10, lty=2)

結(jié)果報(bào)告在圖中。與之前的結(jié)果進(jìn)行比較（以虛線表示）顯示了精度的顯著提高。

回歸分析的通用工具

軟件包dlnm中的功能也可以用作回歸分析的通用工具。第一個(gè)示例演示了如何使用帶有回歸函數(shù)lm（）的回歸樣條來評(píng)估30-39歲的女性樣本中平均身高和體重之間的關(guān)系。

1. > library(splines)

2. > oneheight <- onebasis(women$height, "ns" df=5)

3. > mwomen <- lm(weight ~ oneheight data=women)

使用一個(gè)簡(jiǎn)單的代碼來獲取預(yù)測(cè)和繪圖：

1. > with(pwomen, cbind(fit, low, high)["70",)

2. allfit alllow allhigh

3. 18.92287 18.46545 19.38030

可以簡(jiǎn)單地查看帶有置信區(qū)間的估計(jì)關(guān)聯(lián)，繪制關(guān)聯(lián)。

第二個(gè)示例使用懲罰樣條對(duì)平滑關(guān)聯(lián)進(jìn)行分析。

1. > library(mgcv)

2. 


3. > b2 <- gam(y ~ s(x0,bs="cr") + s(x1,bs="cr") + s(x2,bs="cr") + s(x3,bs="cr"),

4. family=poisson, data=datmethod="REML")

5. > plot(b2, select=3)

該代碼使用通過函數(shù)s（）的回歸樣條，對(duì)帶有多個(gè)變量的模擬數(shù)據(jù)執(zhí)行GAM估計(jì)平滑關(guān)系。也可以使用dlnm獲得預(yù)測(cè)和繪圖，其中：

1. 


2. allRRfit allRRlow allRRhigh

3. 1.3405415 0.8309798 2.1625694

4. > plot(gam, ylim=c(0,3)col=2)

參考文獻(xiàn)?

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

最受歡迎的見解

1.用SPSS估計(jì)HLM層次線性模型模型

2.R語言線性判別分析（LDA），二次判別分析（QDA）和正則判別分析（RDA）

3.基于R語言的lmer混合線性回歸模型

4.R語言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析

5.在r語言中使用GAM（廣義相加模型）進(jìn)行電力負(fù)荷時(shí)間序列分析

6.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM

7.R語言中的嶺回歸、套索回歸、主成分回歸：線性模型選擇和正則化

8.R語言用線性回歸模型預(yù)測(cè)空氣質(zhì)量臭氧數(shù)據(jù)

9.R語言分層線性模型案例