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

最近我們被客戶要求撰寫關于DLM和DLNM的研究報告，包括一些圖形和統(tǒng)計輸出。

本文說明了R語言中實現(xiàn)分布滯后線性和非線性模型（DLM和DLNM）的建模。首先，本文描述了除時間序列數(shù)據(jù)之外的DLM / DLNM的一般化方法，在Gasparrini [2014]中有更詳細的描述?（?點擊文末“閱讀原文”獲取完整代碼數(shù)據(jù)?）?。

前言

本文中包含的結果并不代表科學發(fā)現(xiàn)，而僅出于說明目的進行報告。

數(shù)據(jù)

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

讓我們看一下數(shù)據(jù)框的前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

數(shù)據(jù)集包含來自一項試驗的數(shù)據(jù)，記錄了200名隨機受試者，每名受試者隨機接受四周中兩周的藥物劑量，每天的劑量每周變化。每周7天間隔報告一次暴露水平。數(shù)據(jù)集還包含有關在第28天測量的結果和受試者性別的信息。嵌套的第二個數(shù)據(jù)包括針對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

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

暴露歷史矩陣

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

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

每個受試者的接觸曲線用于重建接觸歷史矩陣。在這種情況下，滯后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年為間隔的暴露量分布圖來嵌套數(shù)據(jù)框的暴露量歷史矩陣。這些數(shù)據(jù)被擴展為滯后3–40的暴露歷史矩陣，滯后單位等于一年。但是，在這種情況下，由于每個對象在不同的年齡進行采樣，因此計算更加復雜。具體地，從受檢者的年齡開始沿著暴露曲線向后計算暴露歷史。此步驟需要一些額外的計算和數(shù)據(jù)處理。可以得出給定時間的暴露曲線的暴露歷史，

>?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歲時進行采樣，則經(jīng)歷了在滯后0處介于80和81之間，在滯后1處介于79和80之間的暴露，依此類推。由于他/她的上一次暴露年齡為65歲，因此將滯后10的暴露歷史記錄設置為0。在69歲時進行采樣的第二個對象的滯后3的暴露歷史記錄設置為0，對應于暴露事件在66。

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

時間序列以外的應用

一個簡單的DLM

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

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

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

CROSSBASIS?FUNCTIONSobservations:?200range:?0?to?100lag?period:?0?27total?df:?4BASIS?FOR?VAR:fun:?nsknots:?9?18

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

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

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

allfit?alllow?allhigh30.29?20.12?40.46

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

可以生成圖：

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

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

點擊標題查閱往期內(nèi)容

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

左右滑動查看更多

01

02

03

04

更為復雜的DLNM?

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

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

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

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

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

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

擴展預測

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

我們也可以計算新的效果摘要，在給定暴露曲線的情況下生成暴露歷史矩陣。例如，我們可以使用嵌套病例對照分析來計算，假設受試者暴露于暴露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

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

>?with(pnestt,?allRfit,aRRlow,allRRigh)20?20?203.5031.2409.900

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

假設某位患者接受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)建所有時間點的暴露歷史矩陣：

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

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

現(xiàn)在可以使用該對象繪制動態(tài)預測：

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

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

應用改進函數(shù)

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

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

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

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

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

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

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

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

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

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

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

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

結果報告在圖中。與之前的結果進行比較（以虛線表示）顯示了精度的顯著提高。

回歸分析的通用工具

軟件包dlnm中的功能也可以用作回歸分析的通用工具。第一個示例演示了如何使用帶有回歸函數(shù)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

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

第二個示例使用懲罰樣條對平滑關聯(liá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（）的回歸樣條，對帶有多個變量的模擬數(shù)據(jù)執(zhí)行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.

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

點擊標題查閱往期內(nèi)容

分布滯后線性和非線性模型（DLNM）分析空氣污染（臭氧）、溫度對死亡率時間序列數(shù)據(jù)的影響
R語言中的分布滯后非線性模型DLNM與發(fā)病率和空氣污染示例
【視頻】R語言中的分布滯后非線性模型（DLNM）與發(fā)病率，死亡率和空氣污染示例
R語言分布滯后線性和非線性模型（DLNM）分析空氣污染（臭氧）、溫度對死亡率時間序列數(shù)據(jù)的影響
R語言分布滯后線性和非線性模型（DLMs和DLNMs）分析時間序列數(shù)據(jù)
R語言分布滯后非線性模型（DLNM）空氣污染研究溫度對死亡率影響建模應用R語言分布滯后非線性模型（DLNM）研究發(fā)病率，死亡率和空氣污染示例
R語言分布滯后線性和非線性模型（DLM和DLNM）建模
R語言廣義相加模型 (GAMs)分析預測CO2時間序列數(shù)據(jù)
Python | ARIMA時間序列模型預測航空公司的乘客數(shù)量
R語言中生存分析模型的時間依賴性ROC曲線可視化
R語言ARIMA，SARIMA預測道路交通流量時間序列分析:季節(jié)性、周期性
ARIMA模型預測CO2濃度時間序列-python實現(xiàn)
R語言基于遞歸神經(jīng)網(wǎng)絡RNN的溫度時間序列預測
R語言用多元ARMA,GARCH ,EWMA, ETS,隨機波動率SV模型對金融時間序列數(shù)據(jù)建模
R語言神經(jīng)網(wǎng)絡模型預測車輛數(shù)量時間序列
卡爾曼濾波器：用R語言中的KFAS建模時間序列
在Python中使用LSTM和PyTorch進行時間序列預測
R語言從經(jīng)濟時間序列中用HP濾波器，小波濾波和經(jīng)驗模態(tài)分解等提取周期性成分分析
使用PYTHON中KERAS的LSTM遞歸神經(jīng)網(wǎng)絡進行時間序列預測
Python中的ARIMA模型、SARIMA模型和SARIMAX模型對時間序列預測
R語言k-Shape時間序列聚類方法對股票價格時間序列聚類
R語言多元Copula GARCH 模型時間序列預測