分布滞後非線性模型（DLNM）表示一個模組化架構，可以靈活地描述在時間序列資料中顯示潛在非線性和滞後影響的關聯。該方法論基于交叉基的定義，交叉基是由兩組基礎函數的組合表示的二維函數空間，它們分别指定了預測變量和滞後變量的關系。本文在R軟體實作DLNM，然後幫助解釋結果，并着重于圖形表示。本文提供指定和解釋DLNM的概念和實踐步驟，并舉例說明了對實際資料的應用
關鍵字：分布滞後模型，時間序列，平滑，滞後效應，R。
1.簡介
統計回歸模型的主要目的是定義一組預測變量與結果之間的關系，然後估計相關影響。當依賴項顯示某些滞後影響時，會進一步增加複雜性：在這種情況下，預測變量的發生（我們稱其為暴露事件）會在遠遠超出事件周期的時間範圍内影響結果。此步驟需要定義更複雜的模型以表征關聯，并指定依賴項的時間結構。
1.1 概念架構
對滞後效應的适當統計模型的說明及其結果的解釋，有助于建立适當的概念架構。這個架構的主要特點是定義了一個額外的次元來描述關聯，它指定了暴露和結果之間在滞後次元上的時間依賴性。這個術語，借用了時間序列分析的文獻，代表了評估影響滞後時暴露事件和結果之間的時間間隔。在長時間暴露的情況下，資料可以通過等距時間段的劃分來構造，定義一系列暴露事件和結果實作。這種劃分也定義了滞後機關。在這個時間結構中，暴露-反應關系可以用兩種相反的觀點中的任何一種來描述：我們可以說一個特定的暴露事件對未來的多個結果産生影響，或者說一個特定的結果可以用過去多個暴露事件的貢獻來解釋。然後，可以使用滞後的概念來描述向前（從固定結果到未來結果）或向後（從固定結果到過去的結果）的關系。
最終，滞後效應統計模型的主要特征是它們的二維結構：該關系同時在預測變量的通常空間和滞後的次元上進行描述。
1.2 分布滞後模型
最近，在評估環境壓力因素的短期影響的研究中已經解決了滞後影響的問題：一些時間序列研究報告說，暴露于高水準的污染或極端溫度會在其發生後的幾天内持續影響健康（ Braga等，2001；Goodman等，2004；Samoli等，2009；Zanobetti和Schwartz，2008）。
給定定義的資料時間結構和簡單的滞後次元定義，時間序列研究設計可提供多種優勢來處理滞後影響，其中時間劃分是由等間隔和有序的時間點直接指定的。在這種情況下，滞後效應可以用分布滞後模型（DLM）來優雅地描述，該模型最初是在計量經濟學中開發的（Almon 1965），最近在環境因素研究中用于量化健康效應（Schwartz 2000; Zanobetti et al。2000; 2007）。Muggeo和Hajat，2009年）。通過這種方法，可以使用多個參數來解釋在不同時滞下的影響，進而将單個暴露事件的影響分布在特定的時間段内，
1.3 本文目的
統計環境R提供了一組用于指定和解釋DLNM結果的工具。本文的目的是提供該程式包函數的全面概述，包括函數的詳細摘要以及以實際資料為例的示例。該示例涉及1987-2000年期間兩個環境因素（空氣污染（臭氧）和溫度）對死亡率的影響。在本文中，我重新考慮了定義DLNM，預測效果并借助圖形函數解釋結果的主要概念和實踐步驟。
2.非線性和滞後效應
在本節中，我介紹了時間序列模型的基本公式，然後介紹了描述非線性效應和滞後效應的方法，後者通過簡單DLM的模型來描述。
2.1 基本模型
時間序列資料的模型通常可以表示為：
其中µt≡E（Yt），Yt是t = 1時的一系列結果...，n，假設來自指數族的分布。函數sj指定變量xj和線性預測變量之間的關系，該變量由參數向量βj定義。變量uk包含具有由相關系數γk指定的線性效應的其他預測變量
之前描述的資料說明性示例中，結果Yt是每日死亡計數，假定是泊松分布，其中E（Y）= µ，V（Y）= φµ。
臭氧和溫度的非線性和滞後影響通過函數sj模組化，該函數定義了預測變量和滞後變量兩個次元之間的關系
2.2 非線性暴露-反應關系
DLNM開發的第一步是定義預測變量空間中的關系。通常，非線性暴露-反應依賴性通過适當的函數s在回歸模型中表示。在完全參數化的方法中，提出了幾種不同的函數，每個函數都具有不同的假設和靈活性。主要選擇通常依賴于描述光滑曲線的函數，例如多項式或樣條函數（Braga等，2001；Dominici等，2004）。關于線性門檻值參數化的使用（Muggeo 2010; Daniels et al。2000）; 或通過虛拟參數化進行簡單分層。
所有這些函數都對原始預測變量進行了轉換，以生成包含在模型中作為線性項的一組轉換變量。相關的基礎函數包括原始變量x的一組完全已知的轉換，這些轉換生成一組稱為基礎變量的新變量。代數表示可以通過以下方式給出：
定義DLNM的第一步是在函數mkbasis（）中執行的，該函數用于建立基礎矩陣Z。此函數的目的是提供一種通用的方式來包含x的非線性效應。舉例來說，我建立了一個将所選基函數應用于向量
的基矩陣：
R> mkais(1:5, tpe = "s", df = 4, egree = 2, cenvlue = 3) 結果是一個清單對象，存儲基礎矩陣和定義該矩陣的自變量。在這種情況下，所選基準是具有4個自由度的二次樣條，由參數類型df和度定義。
可以通過第二個參數類型選擇不同類型的基礎。可用的選項是自然三次方或簡單的B樣條（類型=“ ns”或“ bs”）；虛拟變量層；多項式（“ poly”）；門檻值類型的函數和簡單的線性（“ lin”）。參數df定義了基礎的維數（基礎的列數，基本上是轉換後的變量的數目）。該值可能取決于參數“結點”。如果未定義，則預設情況下将結放置在等距的分位數上。自變量度數選擇“ bs”和“ poly”的多項式度數。
參數cen和cenvalue用于使連續函數（類型“ ns”，“ bs”，“ poly”和“ lin”）的基準居中，如果未提供cenvalue，則預設為原始變量的均值。
2.3滞後效應
定義DLNM的第二步是指定函數，以對附加滞後次元中的關系進行模組化，以實作滞後效果。在這種情況下，給定時間t的結果Yt可以用過去的暴露量xt-L來解釋。給定最大滞後L時，附加滞後次元可以由n×（L +1）矩陣Q表示，例如：
簡單的DLM使用描述結果與滞後風險之間的依賴關系的函數來允許線性關系的滞後效應。
第二步通過函數mklagbasis（）進行，該函數調用mkbasis（）來建構基礎矩陣C。例如：
R> mkgbais(mxlag =5,type ="strta", kots = c(2, 4)) 在此示例中，在通過第一個參數maxlag将最大滞後固定為5之後，滞後向量0：maxlag對應于
，将自動建立并應用所選函數。
3.定義DLNM
DLNM規範的最後一步涉及同時定義預測器和滞後兩個次元中的關系。盡管非線性和滞後效應的術語不同，但這兩個過程在概念上是相似的：定義表示相關空間中關系的基礎。
然後，通過交叉基的定義來指定DLNM，交叉基是二維函數空間，同時描述了沿預測變量範圍及其滞後次元的依存關系。首先，選擇x的基函數得出Z，然後為x的每個基變量建立附加的滞後次元，進而生成一個
數組R˙。通過定義的C，DLNM可以表示為：
選擇交叉基等于如上所述選擇兩組基函數，将其組合以生成交叉基函數。這是通過函數crossbasis（）執行的，該函數調用函數mkbasis（）和mklagbasis（）分别生成兩個基本矩陣Z和C，而不是通過張量積将它們組合起來以産生W。可以使用此函數指定臭氧和溫度的兩個交叉基。相關代碼為：
basi.o3 <- crossbasis(o3 varype= "hthr"
+ vnots = 40, laty = "sata", lanot = c(2,6), mag= 10)
bai.te <- crossbasis(tmp varype = "bs",
+ vrgre 3, vad = 6 cevalu = 25 ladf = 5, malag = 30) 在此示例中，臭氧的交叉基包括一個預測空間的門檻值函數，線性關系超過40.3 µgr / m3，并且虛拟參數化假設沿滞後0-1、2-5和6-10的層具有恒定的分布滞後效應。相比之下，溫度的選項是：以25攝氏度為中心的6 自由度的立方樣條（預設為等距的結點），以及以5自由度的立方樣條（預設為lagtype =“ ns”）（結為25℃）。預設情況下，最多30個滞後。
如果未設定中心值，則預設的中心點是預測變量的平均值（例如，對于上述溫度的交叉基，溫度為25℃）。該值代表來自DLNM的預期效果的參考。參考值的選擇不影響模型的拟合，并且可以根據解釋問題選擇不同的值。
這些選擇可以通過函數summary（）進行檢查。例如：
R> summary(basis.temp) 為了估計相應參數η，可以在通用回歸函數的模型公式中包括交叉基矩陣。在該示例中，最終模型還包括一個自然立方樣條，以模拟季節性趨勢和長期趨勢分量，代碼是：
odel <- glmdeath ~ bais.temp+ basis.o +ns(tim 7 * 14) dw,
+ fmily = quasiposson()) 4.根據DLNM進行預測
如第3節所示，DLNM的規範涉及暴露序列的複雜參數化，但是參數η的估算是使用常見的回歸指令進行的。但是，定義沿兩個次元的關系的此類參數的含義并不簡單。可以通過預測在具有适當暴露值和L + 1滞後的網格上的滞後特定效果來輔助解釋。此外，可以通過将滞後特定貢獻相加來計算從滞後L到0持續暴露所預測的總體效果。預測的效果通過函數crosspred（）在dlnm中計算。以下代碼在示例中計算了對臭氧和溫度的預測：
pre.o <- crosspred(basis, odel at = c(0:6,0., .3)) 傳遞給crosspred（）的前兩個參數是“ crossbasis”類的對象和用于估計的模型對象。像上面的第一個示例一樣，可以通過at參數直接指定必須為其預測效果的暴露值向量。在這裡，我選擇了臭氧中從0到65 µgr / m3的整數，再加上所選門檻值的值和10個機關以上的值（分别為40.3和50.3 µgr / m3）。然後，該函數調用crossbasis（）來建構預測基準，并根據模型中的參數生成預測效果和标準誤差。結果是“ crosspred”類的清單對象，該對象存儲了預測的效果。它包括滞後效應矩陣和總體效應向量，以及相應的标準誤差矩陣和向量。如第5節所示。例如，臭氧增加10個機關的總體效果表示為RR和95％置信區間，可以通過以下公式得出：
R> pred.o3$allRRfit["50.3"] R> cbind(lRlow,alRigh)["50.3",] 5.描述DLNM
由DLNM估算的二維暴露-反應關系可能難以概括。關聯的圖形表示提供了一般描述。調用進階函數plot.default（），persp（）和filled.contour（）來生成散點圖，3-D和等高線圖。例如，臭氧和死亡率之間的關系可以通過RR進行總結，即每次滞後會比門檻值高出10 µgr / m3。該圖如圖1（左）所示，可通過以下方式獲得：
圖1：在門檻值（40.3 µgr / m3）以上的臭氧增加10個機關時，滞後效應（左）和總體效應（右）對死亡率的影響。
R> plot(re.o3) 參數ptype =“ slices”指定圖的類型，在這種情況下，沿着滞後空間在預測值var = 50.3處的預測效果矩陣的切片，對應于在40.3 µgr / m3的門檻值之上增加了10個機關。自變量ci表示置信區間的圖類型。如果使用cumul = TRUE，則繪制累積效果。
根據概念定義，可以使用兩種不同的觀點來讀取圖1中的左圖：它表示在第t天以50.3 µgr / m3的臭氧進行單次暴露後，未來每一天的風險增加。01
02
03
04
或者，可以繪制總體效果，該總體效果是通過使用參數ptype =“ overall”将滞後效應相加得出的：
R> plot(pred ) 圖2：溫度和全因死亡率之間的暴露-反應關系的三維圖，以25°C為參考。
一種更詳細的方法來表示溫度與死亡率之間的平滑關系，其中樣條函數已用于定義這兩個次元的相關性。可以使用3-D和等高線圖對這種複雜的依賴關系進行一般描述，該圖說明了由預測效果的整個網格給出的效果表面。所示的圖是通過以下方式獲得的：
R> plot(pred.temp, "contour") 參考點（此處為25℃）是crossbasis函數在crossbasis（）中中心的值。
三維圖或等高線圖提供了關系的全面摘要，但在表示特定預測值或滞後值的影響方面的能力有限。下面給出了更全面的圖，該圖檔通過以下方式獲得：
R> plot(pred.temp, "slices
+ ci.g , ltensity =20 colr(0))) 圖3（左）顯示了由plot（）和lines（）中的參數var選擇的溫度值的預測滞後效應影響。另外，圖3（右）顯示了針對特定滞後的沿溫度的預測效應的多重曲線圖（左），以及圖3（右）中繪制的相同滞後效應，以及99％的置信區間。
這些圖表顯示了高溫和低溫影響的不同模式，高溫的影響非常強烈且迅速，低溫影響更為延遲，在最初的滞後中為負。
6.模組化政策
DLNM架構提供了機會，可以通過為預測變量和滞後變量兩個次元中的每個次元選擇基本函數來指定廣泛的模型選擇。前面各節中說明的示例代表了一種潛在的模組化替代方法。為了讨論該方法的靈活性以及模型選擇的相關問題，下面顯示了與不同模型的比較，以估計與溫度的關聯。具體來說，為預測變量的空間選擇多項式和層次函數，同時保持相同的自然三次樣條，以模拟長達30天的滞後分布的滞後曲線。指定交叉基礎，運作模型并預測效果的代碼為：
R> basis.temp2 <- crossbasis(emp, vrtpe = "poly",
R> model2 <- update(mdel, .~. - bsis.emp + baiste2)
R> model3 <- updat(model .~. -bais.tmp + bass.mp3) 對于預測變量，第一種方法建議使用與第5節中的原始三次樣條相同的自由度的多項式函數。第二種模型基于一個更簡單的雙門檻值函數，将單個門檻值置于25°C，之前确定為最低死亡率。此選擇還便于模型比較，因為這是其他兩個連續函數的中心點。這三個模型估計的總體效果顯示在由代碼産生的圖4（左）中：
R> plot(pre.temp, "overall", ylim = c(0.5, 2.5), ci = "n", lwd = 1.5,
+ main = "Overall effect")
R> lines(pretemp2, "overall", col = 3, lty = 2, lwd = 2)
R> lines(pretemp3, "overall", col = 4, lty = 4, lwd = 2
+ p, c("natural spline", "polynomial", "double threshold",
+ col = 2:4, lty = c(1:2, 4), lwd = 1.5, inset = 0.1, cex = 0.8) 正如預期的那樣，替代模型會産生不同的結果。特别是，如果與具有等距結點的三次樣條進行比較，則多項式模型會估計出低溫的“擺動”關系。取而代之的是，這兩個函數提供了非常接近的高溫影響估算值。相反，雖然雙門檻值模型的線性假設似乎足以模拟低溫的依賴性，但有一些證據表明，這種方法往往會低估熱的影響。估計的分布滞後曲線的第二次比較如圖4所示（右），如下所示：
R> plot(pred, slices", va =32, im =95 .2="n" 盡管在所有三個模型中都為滞後空間選擇了完全相同的函數，但對預測變量的不同選擇提供了分布滞後曲線的不同估計值，與32°C的參考點相比，代表了32°C的影響。
圖4：溫度為32°C時的總體效應（左）和滞後特異性效應（右）對3種替代模型的全因死亡率的影響（以25°C為參考）。芝加哥1987-2000。
特别是，樣條曲線和多項式模型會産生非常相似的效果（正如預期的那樣，考慮到高溫度尾部曲線在其他次元上的拟合幾乎相同），而雙門檻值模型的曲線顯示出截然不同的形狀。具體而言，由于缺乏此模型的靈活性，是以暗示收獲效果（較長滞後的負估計）可能表示僞像。
缺乏通用标準，無法在可用的選擇中選擇總結關聯的最佳模型，進而減輕了對各種替代産品的規格要求的這種豐富性。在上面的示例中，我對樣條線模型表現出了明顯的偏愛。這種選擇既基于對函數屬性的了解，例如靈活性和穩定性，又基于給出圖4所示結果的合理論據。但是，該結論是有問題的，而不是基于可靠的和一般的統計選擇标準。此外，結論是基于幾個先驗的選擇，就像門檻值位置或結數或多項式次數一樣。
通常，在DLNM中，可以描述兩個不同的選擇級别。第一個涉及不同函數的規範。如上所示，該選擇應既基于假設的暴露反應形狀的合理性，又基于複雜性，可概括性和易于解釋之間的折衷。第二級重點關注特定函數内的不同選擇，例如用于定義樣條曲線基的結的數量和位置。後者更難解決，盡管不是DLNM開發所固有的。一些研究人員在時間序列分析中研究了這個問題，提出了基于資訊準則（Akaike，Bayesian和其他變體），偏自相關或（廣義）交叉驗證的方法（Peng等，2006；Baccini等，2006）。2007）。使用者可以在DLNM中應用相同的方法，但是他應該記住，這些模型的二維性質帶來了額外的複雜性，例如最大滞後的定義。此外，關于執行不同準則的依據還不是結論性的（Dominici等人，2008年）。需要進一步研究以提供有關DLNM中模型選擇的一些指導。
可以建議使用其他方法。Muggeo（2008）提出了一個模型，該模型具有對預測變量空間進行限制的分段參數化，以及基于懲罰性樣條的雙重懲罰基于分布滞後的參數化。此方法包括自動選擇門檻值和分布滞後曲線的平滑度，并且已在R（Muggeo 2010）中完全實作。這種方法與靈活的DLNM的比較可以放寬對預測變量次元上形狀的假設，進而可以提供有關此關系的其他一些見解。
7.資料要求
本文介紹的DLNMs架構是為時間序列資料開發的。（1）中基本模型的一般表達式允許将此方法應用于（廣義）線性模型（GLM）中的任何族分布和連結函數，并擴充到廣義加法模型（GAM）或基于廣義估計方程的模型（GEE）。但是，DLNM的目前實作需要一系列等距，完整和有序的資料。
還使用標明滞後時間段中包含的先前觀察值來計算一系列轉換變量中的每個值。是以，将轉換變量中的第一個最大滞後觀測值設定為NA。允許在x中缺少值，但是由于相同的原因，将相同且下一個maxlag轉換後的值設定為NA。盡管正确，但對于零散的缺失觀測值存在的較長滞後時間的DLNM，這可能會産生計算問題。在這種情況下，可以考慮一些插補方法。