天天看點

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

最近我們被客戶要求撰寫關于DLNM的研究報告,包括一些圖形和統計輸出。 本文提供了運作分布滞後非線性模型的示例,同時描述了預測變量和結果之間的非線性和滞後效應,這種互相關系被定義為暴露-滞後-反應關聯 ( 點選文末“閱讀原文”擷取完整代碼資料******** ) 。

資料

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

該研究是關于随時間變化的職業暴露與癌症之間的關系。該研究包括250個風險集,每個風險集都有一個病例和一個對照,并與年齡相比對。暴露資料以15歲至65歲之間的5歲年齡區間收集。

資料集藥物包含模拟資料,來自一個假設的随機對照試驗,對随時間變化劑量的藥物的影響。該研究包括200名随機受試者,每人每天接受藥物劑量,持續28天,每周都有變化。每隔7天報告一次。

DLNM方法

在這裡,我提供了一個簡短的摘要來介紹概念和定義。

暴露-滞後-反應關聯

DLNM的模組化類用于描述關聯,在該關聯中,暴露和結果之間的依賴關系會在時間上滞後。可以使用兩個不同且互補的觀點來描述此過程。我們可以說,在時間t處的暴露事件确定了在時間t +l處的未來風險。使用後向視角,時間t的風險由過去在時間t-l經曆的一系列風險确定。這裡的l是滞後,表示暴露和測得的結果之間的滞後。

DLNM統計模型

DLNM類提供了一個概念和分析架構,用于描述和估計暴露-滞後-反應關聯。DLNM的統計發展基于以下選擇:DLNM類為描述和估計暴露-滞後-反應關聯提供了一個概念和分析架構。DLNM的統計發展基于該選擇。

暴露-滞後-反應關聯的一個簡單情況是,預測變量空間中的關系(即暴露-滞後關系)是線性的。可以通過DLM對這種類型的關系進行模組化。在這種情況下,關聯僅取決于滞後反應函數,該函數模拟線性風險如何随滞後變化。滞後反應函數的不同選擇(樣條曲線,多項式,層次,門檻值等)導緻指定了不同的DLM,并暗示了滞後反應關系的替代假設。

DLNM解釋

DLNM的結果可以通過使用3-D繪圖提供沿兩個次元變化的關聯,通過為每個滞後和預測變量的拟合值建構預測網格來解釋。

第一是與特定暴露值相關聯的滞後反應曲線,定義為預測變量特定性關聯。這被解釋為與時間t風險相關的時間t +l的風險貢獻序列。

第二是與特定滞後值相關聯的暴露-反應曲線,該特定滞後值定義為滞後特定關聯。這被解釋為與在時間t處發生的暴露值相關聯的在時間t +l處的暴露-反應關系。

第三個也是最重要的是與在考慮的滞後期内經曆的整個暴露曆史相關的暴露反應曲線,定義為總體累積關聯。使用正向視角,這被解釋為表示時間t發生的給定暴露期間[t,t+L]期間經曆的淨風險的暴露反應關系。

時間序列之外的應用

分布滞後模型首先是在很久以前的計量經濟時間序列分析中提出的[Almon,1965],然後在環境流行病學Schwartz [2000]的時間序列資料中重新提出。DLNM的擴充是由Armstrong [2006]構想的。Gasparrini等人對時間序列資料的模組化架構進行了重新評估。[2010]。有趣的是,已經在不同的研究領域中提出了這種暴露-滞後-反應關聯的模型。一般的想法是通過特定函數權重過去的暴露,這些函數的參數由資料估算。在癌症流行病學[Hauptmann等,2000;Langholz等,1999;Richardson,2009;Thomas,1983;Vacek,1997]和藥物流行病學[Abrahamowicz等]中,說明了類似于DLM的線性-暴露-反應關系模型。

基本函數

指定标準暴露反應和滞後反應關系的基本函數,例如多項式,分層或門檻值函數。例如,樣條線由推薦的包樣條線中包含的函數ns()和bs()指定。多項式是通過函數poly()獲得的。這是一個簡單向量的轉換示例:

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.000
attr(,"degree")
[1] 3
attr(,"scale")
[1] 5
attr(,"intercept")
[1] FALSE
attr(,"class")
[1] "poly" "matrix"      

第一個未命名的參數x指定要轉換的向量,而參數度設定多項式的度。定義分層函數是通過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      

結果是帶有附加類别“層”的基礎矩陣。轉換是定義對比的虛拟參數化。參數break定義了層的右開放區間的下邊界。

門檻值函數通過thr()指定。一個例子:

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      

結果是具有附加類别“ thr”的基礎矩陣。參數thr.value定義一個帶有一個或兩個門檻值的向量,而side用于指定高(“ h”,預設值),低(“ l”)或雙精度(“ d”)門檻值參數化。

點選标題查閱往期内容

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

​​【視訊】R語言中的分布滞後非線性模型(DLNM)與發病率,死亡率和空氣污染示例​​

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

左右滑動檢視更多

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

01

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

02

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

03

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

04

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

基本轉換

此函數代表以dlnm為機關進行基本轉換的主要函數,适用于指定暴露-反應和滞後-反應關系。它的作用是應用標明的轉換并以适用于其他函數(例如crossbasis()和crosspred())的格式生成基本矩陣。以下示例複制了該部分中顯示的多項式變換:

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.000
attr(,"fun")
[1] "poly"
attr(,"degree")
[1] 3
attr(,"scale")
[1] 5
attr(,"intercept")
[1] FALSE
attr(,"class")
[1] "onebasis" "matrix"
attr(,"range")
[1] 1 5      

結果是帶有附加類“ onebasis”的基礎矩陣。同樣,第一個未命名參數x指定要轉換的向量,而第二個參數fun将字元轉換定義為應用轉換而調用的函數的名稱。具體來說,基本矩陣包括fun和range屬性,以及定義轉換的被調用函數的參數。如前所述,onebasis()還可以根據特定要求調用使用者定義的函數。一個簡單的例子:

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

交叉基

這是dlnm軟體包中的主要函數。它在内部調用onebasis()來生成暴露-反應和滞後-反應關系的基矩陣,并通過特殊的張量積将它們組合起來,以建立交叉基,該交叉基在模型中同時指定了暴露-滞後-反應關聯性。它的第一個參數x的類定義如何解釋資料。可以使用第二個變量lag修改滞後期。

作為一個簡單的示例,我模拟了2-5個滞後期内3個對象的暴露曆史矩陣:它們中的每一個都将傳遞給onebasis()來分别建構暴露-反應和滞後-反應關系的矩陣。僅用于時間序列資料的附加參數組定義了被視為單獨無關序列的觀察組,例如在季節性分析中可能有用。作為一個簡單的示例,我模拟了2-5個滞後期内3個對象的暴露曆史矩陣:它們中的每一個都将傳遞給onebasis()來分别建構暴露-反應和滞後-反應關系的矩陣。作為一個簡單的示例,我模拟了2-5個滞後期内3個對象的暴露曆史矩陣:

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

然後,我應用交叉基參數化,将二次多項式作為暴露反應函數,并将分層函數2-3和4-5定義為滞後反應函數的分層函數:

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.4236111
sub2 2.333333 1.0833333 1.4583333 0.6736111
sub3 2.916667 1.9166667 2.5625000 1.8402778      

該函數傳回“ crossbasis”類的矩陣對象。它首先使用argvar和arglag清單中的參數調用onebasis(),以建立暴露反應空間和滞後反應空間的矩陣基礎。在另一個示例中,我将crossbasis()應用于資料集中的變量temp,該資料集表示1987-2000年期間日平均溫度序列:

> summary(cb)
CROSSBASIS FUNCTIONS
observations: 5114
range: -26.66667 to 33.33333
lag period: 0 30
total df: 10
BASIS FOR VAR:
fun: thr
thr.value: 10 20
side: d
intercept: FALSE
BASIS FOR LAG:
fun: ns
knots: 1 4 12
intercept: TRUE
Boundary.knots: 0 30      

此處,将暴露反應模組化為門檻值為10和20的雙門檻值函數。滞後時間設定為0到30。滞後反應函數留給預設的自然三次樣條(fun =“ ns”),其滞後值為1、4和12。

預測

crossbasis()生成的交叉基矩陣需要包含在回歸模型公式中才能拟合模型。關聯通過函數crosspred()進行彙總,該函數針對預設值或使用者直接選擇的預測值和滞後值的組合的網格進行預測。例如,我使用建立的交叉基矩陣cb,使用資料集時間序列資料來研究溫度與心血管疾病死亡率之間的關聯。首先,我将一個簡單的線性模型與模型公式中包含的交叉基矩陣拟合。然後,我通過使用cross-basis和回歸模型對象作為前兩個參數調用crosspred()來獲得預測:

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

結果是“ crosspred”類的清單對象,其中的存儲預測和有關模型的其他資訊,例如系數和與交叉基參數有關的關聯(協)方差矩陣的一部分。可以為特定的預測器-滞後組合選擇預測的網格。例如,我提取溫度為-10°C且滞後5的預測和置信區間,然後提取25°C的整體累積預測:

> pred$allfit["25"]
25
1.108262      

第一個結果表明,在給定的一天中,-20°C的溫度會在五天後導緻0.95例心血管死亡的增加,或者在給定的一天中,溫度為-6攝氏度時,心血管死亡的數目增加0.95。其他類型的預測可以通過crosspred()獲得。特别是,如果模型連結等于log或logit,則将自動傳回取幂的預測。如果參數cum設定為TRUE,則是累積預測的矩陣cum。

crosspred()的另一種用法是預測特定的暴露曆史記錄集的影響。這可以通過輸入暴露曆史矩陣作為參數來實作。例如,我們可以從拟合模型中預測出,在過去10天暴露于30°C和在滞後期的其餘時間暴露于22°C之後,心血管死亡的總體累積增加:如果參數cum設定為TRUE,則包括增量累積預測的矩陣cum,并将其存儲在元件cum中。crosspred()的另一種用法是預測特定的暴露曆史記錄集的影響。這可以通過輸入暴露曆史矩陣作為參數來實作。例如,我們可以從拟合模型中預測出,在過去10天暴露于30°C和在滞後期的其餘時間暴露于22°C之後,心血管死亡的總體累積增加:

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

dlnm軟體包的主要優點之一是,使用者可以使用标準回歸函數執行DLNM,隻需在模型公式中包括交叉基矩陣即可。函數crosspred()自動處理來自回歸函數lm()和glm(),gam()(軟體包mgcv),coxph()的模型。

降維

DLNM的拟合度可以降低到預測變量或滞後的一個次元,僅以諸如總累積暴露反應表達。該計算通過函數crossreduce()進行,該函數的工作原理與crosspred()非常相似。前兩個自變量base和model指定交叉基矩陣和需要對其執行計算的模型對象。減少的類型由類型定義,帶有選項“ overall”-“ lag”-“ var”,用于彙總總體累積暴露反應,滞後特異性暴露反應或預測變量特異性滞後反應。

繪圖

一維或二維關聯的解釋通過圖形表示來輔助。通過方法函數plot(),lines()和points()為類“ crosspred”和“ crossreduce”提供進階和低級繪圖功能。例如,我使用對象pred中的預測。plot()方法可以通過參數ptype為“ crosspred”對象生成不同類型的圖。具體來說,它會生成整個二維暴露-滞後-反應關聯的圖形。二維關聯可以繪制為3-D或等高線圖,例如:

> plot(pred,ptype="3d",main="3D plot"      
R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料
R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

可以通過選擇不同的ptype獲得定義的關聯的摘要。

> plot(pred,"overall"      

在這種情況下,方法函數plot()在内部調用函數plot.default(),如上面的示例所示,可以将其他特定參數添加到函數調用中。通過設定ptype =“ slices”,可以将滞後特異性和預測因子特異性關聯分别繪制為暴露-反應和滞後-反應關系,因為它們是在3-D曲面中沿特定次元切割的切片。例如:

> plot(pred,"slices",lag=5      
R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料
R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料
R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

這兩個圖分别代表了滞後5的暴露反應和特定于25°C溫度的滞後反應。參數lag和var指定必須分别繪制lag和特定于預測變量的關聯的值。

R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例|附代碼資料

擷取全文完整代碼資料資料。

本文選自《R語言分布滞後非線性模型(DLNM)研究發病率,死亡率和空氣污染示例》。

繼續閱讀