全文下載鏈接：http://tecdat.cn/?p=18700

本文說明了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è)觀察樣本：

>?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

數(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é)果是：

>?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

變量病例定義病例/對(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è)受試者的接觸曲線用于重建接觸歷史矩陣。

>?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

上面針對(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í)間的暴露曲線的暴露歷史，

>?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

上面針對(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獲得摘要：

CROSSBASIS?FUNCTIONS observations:?200 range:?0?to?100 lag?period:?0?27 total?df:?4 BASIS?FOR?VAR: fun:?ns 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中

allfit?alllow?allhigh 30.29?20.12?40.46

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

可以生成圖：

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

正在上傳…重新上傳取消

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

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語言里的非線性模型：多項(xiàng)式回歸、局部樣條、平滑樣條、 廣義相加模型GAM分析

左右滑動(dòng)查看更多

01

02

03

04

更為復(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è)效果摘要。代碼是：

>?library(survival) 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è)。

>?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

產(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ì)中使用的一致。

>?with(pnestt,?allRfit,aRRlow,allRRigh) 20?20?20 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è)：

>?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)

在圖中繪制了整體累積關(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)曲線：

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

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

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

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

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

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

>?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

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

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

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

回歸分析的通用工具

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

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

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

>?with(pwomen,?cbind(fit,?low,?high)\["70",) allfit?alllow?allhigh 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)行分析。

>?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)

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

allRRfit?allRRlow?allRRhigh 1.3405415?0.8309798?2.1625694 >?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.

本文摘選《R語言分布滯后線性和非線性模型（DLM和DLNM）建模》，點(diǎn)擊“閱讀原文”獲取全文完整資料。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

分布滯后線性和非線性模型（DLNM）分析空氣污染（臭氧）、溫度對(duì)死亡率時(shí)間序列數(shù)據(jù)的影響

R語言中的分布滯后非線性模型DLNM與發(fā)病率和空氣污染示例

【視頻】R語言中的分布滯后非線性模型（DLNM）與發(fā)病率，死亡率和空氣污染示例

R語言分布滯后線性和非線性模型（DLNM）分析空氣污染（臭氧）、溫度對(duì)死亡率時(shí)間序列數(shù)據(jù)的影響

R語言分布滯后線性和非線性模型（DLMs和DLNMs）分析時(shí)間序列數(shù)據(jù)

R語言分布滯后非線性模型（DLNM）空氣污染研究溫度對(duì)死亡率影響建模應(yīng)用

R語言分布滯后非線性模型（DLNM）研究發(fā)病率，死亡率和空氣污染示例

R語言分布滯后線性和非線性模型（DLM和DLNM）建模

R語言廣義相加模型 (GAMs)分析預(yù)測(cè)CO2時(shí)間序列數(shù)據(jù)

Python | ARIMA時(shí)間序列模型預(yù)測(cè)航空公司的乘客數(shù)量

R語言中生存分析模型的時(shí)間依賴性ROC曲線可視化

R語言ARIMA，SARIMA預(yù)測(cè)道路交通流量時(shí)間序列分析:季節(jié)性、周期性

ARIMA模型預(yù)測(cè)CO2濃度時(shí)間序列-python實(shí)現(xiàn)

R語言基于遞歸神經(jīng)網(wǎng)絡(luò)RNN的溫度時(shí)間序列預(yù)測(cè)

R語言用多元ARMA,GARCH ,EWMA, ETS,隨機(jī)波動(dòng)率SV模型對(duì)金融時(shí)間序列數(shù)據(jù)建模

R語言神經(jīng)網(wǎng)絡(luò)模型預(yù)測(cè)車輛數(shù)量時(shí)間序列

卡爾曼濾波器：用R語言中的KFAS建模時(shí)間序列

在Python中使用LSTM和PyTorch進(jìn)行時(shí)間序列預(yù)測(cè)

R語言從經(jīng)濟(jì)時(shí)間序列中用HP濾波器，小波濾波和經(jīng)驗(yàn)?zāi)B(tài)分解等提取周期性成分分析

使用PYTHON中KERAS的LSTM遞歸神經(jīng)網(wǎng)絡(luò)進(jìn)行時(shí)間序列預(yù)測(cè)

Python中的ARIMA模型、SARIMA模型和SARIMAX模型對(duì)時(shí)間序列預(yù)測(cè)

R語言k-Shape時(shí)間序列聚類方法對(duì)股票價(jià)格時(shí)間序列聚類

R語言多元Copula GARCH 模型時(shí)間序列預(yù)測(cè)