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

最近我們被客戶(hù)要求撰寫(xiě)關(guān)于分布滯后非線性模型（DLNM）的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。

本文提供了運(yùn)行分布滯后非線性模型的示例，同時(shí)描述了預(yù)測(cè)變量和結(jié)果之間的非線性和滯后效應(yīng)，這種相互關(guān)系被定義為暴露-滯后-反應(yīng)關(guān)聯(lián)

數(shù)據(jù)

數(shù)據(jù)集包含1987-2000年期間每日死亡率（CVD、呼吸道），天氣（溫度，相對(duì)濕度）和污染數(shù)據(jù)（PM10和臭氧）。數(shù)據(jù)是由健康影響研究所贊助的《國(guó)家發(fā)病率，死亡率和空氣污染研究》（NMMAPS）的一部分[Samet et al.,2000a,b]。

該研究是關(guān)于隨時(shí)間變化的職業(yè)暴露與癌癥之間的關(guān)系。該研究包括250個(gè)風(fēng)險(xiǎn)集，每個(gè)風(fēng)險(xiǎn)集都有一個(gè)病例和一個(gè)對(duì)照，并與年齡相匹配。暴露數(shù)據(jù)以15歲至65歲之間的5歲年齡區(qū)間收集。

數(shù)據(jù)集藥物包含模擬數(shù)據(jù)，來(lái)自一個(gè)假設(shè)的隨機(jī)對(duì)照試驗(yàn)，對(duì)隨時(shí)間變化劑量的藥物的影響。該研究包括200名隨機(jī)受試者，每人每天接受藥物劑量，持續(xù)28天，每周都有變化。每隔7天報(bào)告一次。

DLNM方法

在這里，我提供了一個(gè)簡(jiǎn)短的摘要來(lái)介紹概念和定義。

暴露-滯后-反應(yīng)關(guān)聯(lián)

DLNM的建模類(lèi)用于描述關(guān)聯(lián)，在該關(guān)聯(lián)中，暴露和結(jié)果之間的依賴(lài)關(guān)系會(huì)在時(shí)間上滯后?？梢允褂脙蓚€(gè)不同且互補(bǔ)的觀點(diǎn)來(lái)描述此過(guò)程。我們可以說(shuō)，在時(shí)間t處的暴露事件確定了在時(shí)間t +l處的未來(lái)風(fēng)險(xiǎn)。使用后向視角，時(shí)間t的風(fēng)險(xiǎn)由過(guò)去在時(shí)間t-l經(jīng)歷的一系列風(fēng)險(xiǎn)確定。這里的l是滯后，表示暴露和測(cè)得的結(jié)果之間的滯后。

DLNM統(tǒng)計(jì)模型

DLNM類(lèi)提供了一個(gè)概念和分析框架，用于描述和估計(jì)暴露-滯后-反應(yīng)關(guān)聯(lián)。DLNM的統(tǒng)計(jì)發(fā)展基于以下選擇：DLNM類(lèi)為描述和估計(jì)暴露-滯后-反應(yīng)關(guān)聯(lián)提供了一個(gè)概念和分析框架。DLNM的統(tǒng)計(jì)發(fā)展基于該選擇。

暴露-滯后-反應(yīng)關(guān)聯(lián)的一個(gè)簡(jiǎn)單情況是，預(yù)測(cè)變量空間中的關(guān)系（即暴露-滯后關(guān)系）是線性的。可以通過(guò)DLM對(duì)這種類(lèi)型的關(guān)系進(jìn)行建模。在這種情況下，關(guān)聯(lián)僅取決于滯后反應(yīng)函數(shù)，該函數(shù)模擬線性風(fēng)險(xiǎn)如何隨滯后變化。滯后反應(yīng)函數(shù)的不同選擇（樣條曲線，多項(xiàng)式，層次，閾值等）導(dǎo)致指定了不同的DLM，并暗示了滯后反應(yīng)關(guān)系的替代假設(shè)。

DLNM解釋

DLNM的結(jié)果可以通過(guò)使用3-D繪圖提供沿兩個(gè)維度變化的關(guān)聯(lián)，通過(guò)為每個(gè)滯后和預(yù)測(cè)變量的擬合值構(gòu)建預(yù)測(cè)網(wǎng)格來(lái)解釋。

第一是與特定暴露值相關(guān)聯(lián)的滯后反應(yīng)曲線，定義為預(yù)測(cè)變量特定性關(guān)聯(lián)。這被解釋為與時(shí)間t風(fēng)險(xiǎn)相關(guān)的時(shí)間t +l的風(fēng)險(xiǎn)貢獻(xiàn)序列。

第二是與特定滯后值相關(guān)聯(lián)的暴露-反應(yīng)曲線，該特定滯后值定義為滯后特定關(guān)聯(lián)。這被解釋為與在時(shí)間t處發(fā)生的暴露值相關(guān)聯(lián)的在時(shí)間t +l處的暴露-反應(yīng)關(guān)系。

第三個(gè)也是最重要的是與在考慮的滯后期內(nèi)經(jīng)歷的整個(gè)暴露歷史相關(guān)的暴露反應(yīng)曲線，定義為總體累積關(guān)聯(lián)。使用正向視角，這被解釋為表示時(shí)間t發(fā)生的給定暴露期間[t，t+L]期間經(jīng)歷的凈風(fēng)險(xiǎn)的暴露反應(yīng)關(guān)系。

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

分布滯后模型首先是在很久以前的計(jì)量經(jīng)濟(jì)時(shí)間序列分析中提出的[Almon，1965]，然后在環(huán)境流行病學(xué)Schwartz [2000]的時(shí)間序列數(shù)據(jù)中重新提出。DLNM的擴(kuò)展是由Armstrong [2006]構(gòu)想的。Gasparrini等人對(duì)時(shí)間序列數(shù)據(jù)的建模框架進(jìn)行了重新評(píng)估。[2010]。有趣的是，已經(jīng)在不同的研究領(lǐng)域中提出了這種暴露-滯后-反應(yīng)關(guān)聯(lián)的模型。一般的想法是通過(guò)特定函數(shù)加權(quán)過(guò)去的暴露，這些函數(shù)的參數(shù)由數(shù)據(jù)估算。在癌癥流行病學(xué)[Hauptmann等，2000；Langholz等，1999；Richardson，2009；Thomas，1983；Vacek，1997]和藥物流行病學(xué)[Abrahamowicz等]中，說(shuō)明了類(lèi)似于DLM的線性-暴露-反應(yīng)關(guān)系模型。

基本函數(shù)

指定標(biāo)準(zhǔn)暴露反應(yīng)和滯后反應(yīng)關(guān)系的基本函數(shù)，例如多項(xiàng)式，分層或閾值函數(shù)。例如，樣條線由推薦的包樣條線中包含的函數(shù)ns（）和bs（）指定。多項(xiàng)式是通過(guò)函數(shù)poly（）獲得的。這是一個(gè)簡(jiǎn)單向量的轉(zhuǎn)換示例：

poly(1:5,degree=3)1 2 3[1,] 0.2 0.04 0.008[2,] 0.4 0.16 0.064[3,] 0.6 0.36 0.216[4,] 0.8 0.64 0.512[5,] 1.0 1.00 1.000attr(,"degree")[1] 3attr(,"scale")[1] 5attr(,"intercept")[1] FALSEattr(,"class")[1] "poly" "matrix"

第一個(gè)未命名的參數(shù)x指定要轉(zhuǎn)換的向量，而參數(shù)度設(shè)置多項(xiàng)式的度。定義分層函數(shù)是通過(guò)strata（）指定的。

strata(1:5,breaks=c(2,4))[,]1 2[1,] 0 0[2,] 1 0[3,] 1 0[4,] 0 1[5,] 0 1

結(jié)果是帶有附加類(lèi)別“層”的基礎(chǔ)矩陣。轉(zhuǎn)換是定義對(duì)比的虛擬參數(shù)化。參數(shù)break定義了層的右開(kāi)放區(qū)間的下邊界。

閾值函數(shù)通過(guò)thr（）指定。一個(gè)例子：

thr(1:5,thr.value=3,side="d")[,] 1 2[1,] 2 0[2,] 1 0[3,] 0 0[4,] 0 1[5,] 0 2

結(jié)果是具有附加類(lèi)別“ thr”的基礎(chǔ)矩陣。參數(shù)thr.value定義一個(gè)帶有一個(gè)或兩個(gè)閾值的向量，而side用于指定高（“ h”，默認(rèn)值），低（“ l”）或雙精度（“ d”）閾值參數(shù)化。

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

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

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

01

02

03

04

基本轉(zhuǎn)換

此函數(shù)代表以dlnm為單位進(jìn)行基本轉(zhuǎn)換的主要函數(shù)，適用于指定暴露-反應(yīng)和滯后-反應(yīng)關(guān)系。它的作用是應(yīng)用選定的轉(zhuǎn)換并以適用于其他函數(shù)（例如crossbasis（）和crosspred（））的格式生成基本矩陣。以下示例復(fù)制了該部分中顯示的多項(xiàng)式變換：

onebasis(1:5,fun="poly",degree=3) b1 b2 b3[1,] 0.2 0.04 0.008[2,] 0.4 0.16 0.064[3,] 0.6 0.36 0.216[4,] 0.8 0.64 0.512[5,] 1.0 1.00 1.000attr(,"fun")[1] "poly"attr(,"degree")[1] 3attr(,"scale")[1] 5attr(,"intercept")[1] FALSEattr(,"class")[1] "onebasis" "matrix"attr(,"range")[1] 1 5

結(jié)果是帶有附加類(lèi)“ onebasis”的基礎(chǔ)矩陣。同樣，第一個(gè)未命名參數(shù)x指定要轉(zhuǎn)換的向量，而第二個(gè)參數(shù)fun將字符轉(zhuǎn)換定義為應(yīng)用轉(zhuǎn)換而調(diào)用的函數(shù)的名稱(chēng)。具體來(lái)說(shuō)，基本矩陣包括fun和range屬性，以及定義轉(zhuǎn)換的被調(diào)用函數(shù)的參數(shù)。如前所述，onebasis（）還可以根據(jù)特定要求調(diào)用用戶(hù)定義的函數(shù)。一個(gè)簡(jiǎn)單的例子：

> mylog <- function(x) log(x) > onebasis(1:5,"mylog") b1[1,] 0.0000000[2,] 0.6931472[3,] 1.0986123[4,] 1.3862944[5,] 1.6094379attr(,"fun")[1] "mylog"attr(,"range")[1] 1 5attr(,"class")[1] "onebasis" "matrix"

交叉基

這是dlnm軟件包中的主要函數(shù)。它在內(nèi)部調(diào)用onebasis（）來(lái)生成暴露-反應(yīng)和滯后-反應(yīng)關(guān)系的基矩陣，并通過(guò)特殊的張量積將它們組合起來(lái)，以創(chuàng)建交叉基，該交叉基在模型中同時(shí)指定了暴露-滯后-反應(yīng)關(guān)聯(lián)性。它的第一個(gè)參數(shù)x的類(lèi)定義如何解釋數(shù)據(jù)?？梢允褂玫诙€(gè)變量lag修改滯后期。

作為一個(gè)簡(jiǎn)單的示例，我模擬了2-5個(gè)滯后期內(nèi)3個(gè)對(duì)象的暴露歷史矩陣：它們中的每一個(gè)都將傳遞給onebasis（）來(lái)分別構(gòu)建暴露-反應(yīng)和滯后-反應(yīng)關(guān)系的矩陣。僅用于時(shí)間序列數(shù)據(jù)的附加參數(shù)組定義了被視為單獨(dú)無(wú)關(guān)序列的觀察組，例如在季節(jié)性分析中可能有用。作為一個(gè)簡(jiǎn)單的示例，我模擬了2-5個(gè)滯后期內(nèi)3個(gè)對(duì)象的暴露歷史矩陣：它們中的每一個(gè)都將傳遞給onebasis（）來(lái)分別構(gòu)建暴露-反應(yīng)和滯后-反應(yīng)關(guān)系的矩陣。作為一個(gè)簡(jiǎn)單的示例，我模擬了2-5個(gè)滯后期內(nèi)3個(gè)對(duì)象的暴露歷史矩陣：

> histlag2 lag3 lag4 lag5 sub1 1 3 5 6 sub2 7 8 9 4 sub3 10 2 11 12

然后，我應(yīng)用交叉基參數(shù)化，將二次多項(xiàng)式作為暴露反應(yīng)函數(shù)，并將分層函數(shù)2-3和4-5定義為滯后反應(yīng)函數(shù)的分層函數(shù)：

lag=c(2,5),argvar=list(fun="poly",degree=2),arglag=list(fun="strata",breaks=4))[,]v1.l1 v1.l2 v2.l1 v2.l2 sub1 1.250000 0.9166667 0.4930556 0.4236111sub2 2.333333 1.0833333 1.4583333 0.6736111sub3 2.916667 1.9166667 2.5625000 1.8402778

該函數(shù)返回“ crossbasis”類(lèi)的矩陣對(duì)象。它首先使用argvar和arglag列表中的參數(shù)調(diào)用onebasis（），以建立暴露反應(yīng)空間和滯后反應(yīng)空間的矩陣基礎(chǔ)。在另一個(gè)示例中，我將crossbasis（）應(yīng)用于數(shù)據(jù)集中的變量temp，該數(shù)據(jù)集表示1987-2000年期間日平均溫度序列：

> summary(cb)CROSSBASIS FUNCTIONSobservations: 5114range: -26.66667 to 33.33333lag period: 0 30total df: 10BASIS FOR VAR:fun: thrthr.value: 10 20side: dintercept: FALSEBASIS FOR LAG:fun: nsknots: 1 4 12intercept: TRUEBoundary.knots: 0 30

此處，將暴露反應(yīng)建模為閾值為10和20的雙閾值函數(shù)。滯后時(shí)間設(shè)置為0到30。滯后反應(yīng)函數(shù)留給默認(rèn)的自然三次樣條（fun =“ ns”），其滯后值為1、4和12。

預(yù)測(cè)

crossbasis（）生成的交叉基矩陣需要包含在回歸模型公式中才能擬合模型。關(guān)聯(lián)通過(guò)函數(shù)crosspred（）進(jìn)行匯總，該函數(shù)針對(duì)默認(rèn)值或用戶(hù)直接選擇的預(yù)測(cè)值和滯后值的組合的網(wǎng)格進(jìn)行預(yù)測(cè)。例如，我使用創(chuàng)建的交叉基矩陣cb，使用數(shù)據(jù)集時(shí)間序列數(shù)據(jù)來(lái)研究溫度與心血管疾病死亡率之間的關(guān)聯(lián)。首先，我將一個(gè)簡(jiǎn)單的線性模型與模型公式中包含的交叉基矩陣擬合。然后，我通過(guò)使用cross-basis和回歸模型對(duì)象作為前兩個(gè)參數(shù)調(diào)用crosspred（）來(lái)獲得預(yù)測(cè)：

crosspred(cb,model,at=-20:30)

結(jié)果是“ crosspred”類(lèi)的列表對(duì)象，其中的存儲(chǔ)預(yù)測(cè)和有關(guān)模型的其他信息，例如系數(shù)和與交叉基參數(shù)有關(guān)的關(guān)聯(lián)（協(xié)）方差矩陣的一部分?？梢詾樘囟ǖ念A(yù)測(cè)器-滯后組合選擇預(yù)測(cè)的網(wǎng)格。例如，我提取溫度為-10°C且滯后5的預(yù)測(cè)和置信區(qū)間，然后提取25°C的整體累積預(yù)測(cè)：

> pred$allfit["25"]25 1.108262

第一個(gè)結(jié)果表明，在給定的一天中，-20°C的溫度會(huì)在五天后導(dǎo)致0.95例心血管死亡的增加，或者在給定的一天中，溫度為-6攝氏度時(shí)，心血管死亡的數(shù)目增加0.95。其他類(lèi)型的預(yù)測(cè)可以通過(guò)crosspred（）獲得。特別是，如果模型鏈接等于log或logit，則將自動(dòng)返回取冪的預(yù)測(cè)。如果參數(shù)cum設(shè)置為T(mén)RUE，則是累積預(yù)測(cè)的矩陣cum。

crosspred（）的另一種用法是預(yù)測(cè)特定的暴露歷史記錄集的影響。這可以通過(guò)輸入暴露歷史矩陣作為參數(shù)來(lái)實(shí)現(xiàn)。例如，我們可以從擬合模型中預(yù)測(cè)出，在過(guò)去10天暴露于30°C和在滯后期的其余時(shí)間暴露于22°C之后，心血管死亡的總體累積增加：如果參數(shù)cum設(shè)置為T(mén)RUE，則包括增量累積預(yù)測(cè)的矩陣cum，并將其存儲(chǔ)在組件cum中。crosspred（）的另一種用法是預(yù)測(cè)特定的暴露歷史記錄集的影響。這可以通過(guò)輸入暴露歷史矩陣作為參數(shù)來(lái)實(shí)現(xiàn)。例如，我們可以從擬合模型中預(yù)測(cè)出，在過(guò)去10天暴露于30°C和在滯后期的其余時(shí)間暴露于22°C之后，心血管死亡的總體累積增加：

> crosspred(cb,model,at=histpred)$allfit1 5.934992

dlnm軟件包的主要優(yōu)點(diǎn)之一是，用戶(hù)可以使用標(biāo)準(zhǔn)回歸函數(shù)執(zhí)行DLNM，只需在模型公式中包括交叉基矩陣即可。函數(shù)crosspred（）自動(dòng)處理來(lái)自回歸函數(shù)lm（）和glm（），gam（）（軟件包mgcv），coxph（）的模型。

降維

DLNM的擬合度可以降低到預(yù)測(cè)變量或滯后的一個(gè)維度，僅以諸如總累積暴露反應(yīng)表達(dá)。該計(jì)算通過(guò)函數(shù)crossreduce（）進(jìn)行，該函數(shù)的工作原理與crosspred（）非常相似。前兩個(gè)自變量base和model指定交叉基矩陣和需要對(duì)其執(zhí)行計(jì)算的模型對(duì)象。減少的類(lèi)型由類(lèi)型定義，帶有選項(xiàng)“ overall”-“ lag”-“ var”，用于匯總總體累積暴露反應(yīng)，滯后特異性暴露反應(yīng)或預(yù)測(cè)變量特異性滯后反應(yīng)。

繪圖

一維或二維關(guān)聯(lián)的解釋通過(guò)圖形表示來(lái)輔助。通過(guò)方法函數(shù)plot（），lines（）和points（）為類(lèi)“ crosspred”和“ crossreduce”提供高級(jí)和低級(jí)繪圖功能。例如，我使用對(duì)象pred中的預(yù)測(cè)。plot（）方法可以通過(guò)參數(shù)ptype為“ crosspred”對(duì)象生成不同類(lèi)型的圖。具體來(lái)說(shuō)，它會(huì)生成整個(gè)二維暴露-滯后-反應(yīng)關(guān)聯(lián)的圖形。二維關(guān)聯(lián)可以繪制為3-D或等高線圖，例如：

>?plot(pred,ptype="3d",main="3D plot"

可以通過(guò)選擇不同的ptype獲得定義的關(guān)聯(lián)的摘要。

>?plot(pred,"overall"

在這種情況下，方法函數(shù)plot（）在內(nèi)部調(diào)用函數(shù)plot.default（），如上面的示例所示，可以將其他特定參數(shù)添加到函數(shù)調(diào)用中。通過(guò)設(shè)置ptype =“ slices”，可以將滯后特異性和預(yù)測(cè)因子特異性關(guān)聯(lián)分別繪制為暴露-反應(yīng)和滯后-反應(yīng)關(guān)系，因?yàn)樗鼈兪窃?-D曲面中沿特定維度切割的切片。例如：

>?plot(pred,"slices",lag=5

這兩個(gè)圖分別代表了滯后5的暴露反應(yīng)和特定于25°C溫度的滯后反應(yīng)。參數(shù)lag和var指定必須分別繪制lag和特定于預(yù)測(cè)變量的關(guān)聯(lián)的值。

點(diǎn)擊文末?“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語(yǔ)言分布滯后非線性模型（DLNM）研究發(fā)病率，死亡率和空氣污染示例》。

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

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