劑量-反應Meta分析中，探討自變量與因變量之間的因果關系通常的做法是建立回歸模型。經典的模型中需同時考慮線性及非線性趨勢。傳統的線性模型是基于一次函數的簡單模型，而這種模型適用性較差。本文在線性模型基礎上，通過引入分段函數，建立分段線性劑量-反應Meta分析模型，很好地彌補了傳統線性模型的不足。本文將詳細探討劑量-反應Meta分析中線性及分段線性模型的方法學及應用。
引用本文： 于云祥, 張超, 翁鴻, 鄺心穎, 龔泰芳, 盧云, 徐暢. 劑量-反應Meta分析之線性關系及分段線性回歸模型的應用. 中國循證醫學雜志, 2016, 16(1): 111-114. doi: 10.7507/1672-2531.20160020 復制
在探討自變量與因變量之間因果關系時，通常可以建立一個以自變量為X，應變量為Y的回歸方程 [1]。通過建立的回歸模型，我們可以用觀察值估計線性或非線性趨勢。劑量-反應Meta分析模型也是一種回歸方程，既包含線性模型也包含非線性模型。線性模型即一次函數模型，反應自變量與因變量之間的線性關系；通過對線性模型引入二次項和（或）三次項甚至高次項，可用于擬合自變量與因變量之間的非線性關系 [2]。在循證實踐中，這些非線性模型在擬合非線性趨勢時往往能取得較好的逼近效果，并且可以通過改變模型中的函數項的次數或者通過插值等方法對非線性趨勢進行調整，能適用于大多數情況。目前劑量-反應Meta分析中的線性模型對線性趨勢預測效果則欠佳，特別是當自變量與因變量呈J-型、U-型、S-型、V-型等關系時，其趨勢并非單一走向，傳統線性模型無法準確逼近。本文引入一類新的線性模型——分段線性模型，可以很好地解決這個問題。以下從基礎知識、模型、方法學三方面進行探討，并輔以實例演示，與傳統線性Meta分析模型進行比較。
1 傳統線性模型
1.1 模型表達式
我們在《劑量-反應Meta分析方法學概論》一文 [3]中已對線性模型進行了描述，即引入一個線性回歸函數來反應暴露水平（自變量）與結局指標（因變量）發生風險的線性關系，表達式如下：
| ${{y}_{j}}={{\beta }_{1}}{{x}_{j}}+{{\varepsilon }_{j~}}$ |
yj表示在結局指標發生風險在第j層暴露水平的取值，期望值為β1xj；xj則對應第j層暴露水平的取值。β1是待估參數，也即回歸直線的斜率。εj表示yj的隨機誤差項，即標準誤，滿足均值E(ε)=0，值得注意的是，其方差不等且協方差不為0 [4]。該函數被強制通過原點，也就是在參考劑量暴露水平時對應零風險。繼續以i表示納入研究，Meta分析模型則可表示為：
| ${{y}_{ij}}={{\beta }_{i1}}{{x}_{ij}}+{{\varepsilon }_{ij}}$ |
因為該模型并不完全滿足回歸模型的基本假定，因此參數估計不能使用普通最小二乘法，而需要通過廣義最小二乘法 [4]方可計算單個研究中的回歸斜率及其置信區間。因變量的隨機誤差可通過《劑量-反應Meta分析之限制性立方樣條函數在模型中的應用》一文 [5]中給出的公式求出。在此基礎之上，對單個研究的斜率進行加權合并，得到合并后的回歸斜率?β，如下：
| $\hat{\beta }\text{= }\frac{{{\sum{_{\text{ij}}{{\text{x}}_{\text{ij}}}{{\text{y}}_{\text{ij}}}\text{/}{{\varepsilon }_{\text{ij}}}}}^{2}}}{\sum{_{\text{ij}}{{\text{x}}_{\text{ij}}}^{2}\text{/}{{\varepsilon }_{\text{ij}}}^{2}}}$ |
| ${{\text{S}}_{\beta }}^{\text{2}}\text{= }\frac{\text{1}}{\sum{_{\text{ij}}{{\text{x}}_{\text{ij}}}^{2}\text{/}{{\varepsilon }_{\text{ij}}}^{2}}}$ |
求出斜率后，就可以知道yij的期望值，實際上，我們最終感興趣的結果也即這個期望值。該公式也可以簡化為矩陣形式寫出：
| $\hat{\beta }={{\left( X'{{V}^{-1}}X \right)}^{-1}}X'{{V}^{-1}}Y$ |
| $S{{~}_{\beta }}^{2}={{\left( X'V-1X \right)}^{-1}}$ |
1.2 模型的意義
在循證實踐中，單個研究的暴露水平分布往往不同，如最低暴露劑量不一致，導致函數中暴露基線差異大。為避免此種情況的發生，在進行線性劑量-反應Meta分析時，通常以合并后的斜率來表示自變量與因變量之間的關系。以公式（2）為例，x每變化一個單位，對應y則變化個單位。當然，這種變化可以顯著也可以不顯著，可根據假設檢驗來判斷：H0：=0，計算該零假設成立的情況下，觀察結果在樣本分布中可能發生的概率 [1]。
1.3 操作代碼
多種軟件可以進行線性劑量反應Meta分析，我們以Stata/SE12.0軟件為例，實現代碼如下：
encode study,gen(studynum)
gen double logrr=log(adjrr)
gen double logub=log(ub)
gen double loglb=log(lb)
gen double se=(logub-loglb)/(2*invnormal(0.975) )
glstlogrr dose,se(se) cov(n case) pfirst(id studynum) ts(r)
lincom dose*1,eform
關于代碼的意義及解釋請參考《應用STATA做Meta分析》一書 [6]相關章節及《應用 R軟件dosresmeta程序包和mvmeta程序包實現劑量-反應關系Meta分析》一文 [2]相關內容。
2 分段線性模型
分段線性模型主要是針對存在分段線性趨勢的分布的情況。模型可以用多種函數表示，為方便讀者理解，我們與系列文章 [5]保持一致，使用分段線性樣條函數 [7]。也是對經典限制性立方樣條函數法下的劑量-反應Meta分析模型缺陷的一種補充。值得注意的是，該函數并非本團隊原創，我們的工作只是將該函數引入劑量-反應Meta分析，并研發出了相關代碼。
2.1 模型的表達式
跟限制性立方樣條函數類似，直線樣條函數也是通過節點連接的分段直線函數，需滿足每個節點處連續。假設以ki表示節點（i=1，2，…，n-1），以單個研究為例，模型表達式如下：
| $y=g\left( x \right)={{b}_{1i}}{{X}_{i}}+\varepsilon $ |
式中，b1k為回歸斜率，并且隨著節點不同而改變，Xi是X的函數，意義如下：
| ${{X}_{1}}=X$ |
| ${{X}_{i+1}}=max\left( X-{{k}_{i-1}} \right)$ |
X的值大于或等于ki-1，公式（9）為X-ki-1，反之，則為0。為便于理解，我們將公式進行變換，形成一個類似多元線性回歸的模型：
| $y=g\left( x \right)={{b}_{1}}{{X}_{1}}+{{b}_{1i+1}}{{X}_{i+1}}+\varepsilon ~$ |
以j（j=1，2，3…）表示單個研究，Meta分析模型則如下：
| ${{y}_{j}}={{\beta }_{1ij}}{{x}_{ij}}+{{\varepsilon }_{j}}$ |
2.2 模型的意義
分段線性函數模型的意義與普通線性函數類似，不同的是通過不同節點將該函數分成了多段連續的線性函數，這些由節點分成的子區間的線性函數，斜率可同可不同，因此可更靈活地擬合線性趨勢。
仍然以X的分布上10th（k1）、50th（k2）、90th（k3）這3個節點為例，3個節點將整個趨勢分為了4個小段，即0th～10th、10th～50th、50th～90th、90th～100th。函數表達式如下：
y=g(x)=b1X1+ε；0th≤X1＜10th
y=g(x)=b12X2+ε；10th≤X1＜50th
y=g(x)=b13X3+ε；50th≤X1＜90th
y=g(x)=b14X4+ε；90th≤X1＜100th
根據廣義最小二乘法逐一計算出不同分段的斜率，即可反應自變量與因變量之間的分段線性關系。即在每一個分段內，x每變化一個單位，對應y則變化b1i個單位。
2.3 操作代碼
本文給出主要操作代碼，值得注意的是，該代碼由本研究團隊根據簡單回歸模型、多重回歸模型及非線性Meta分析模型的代碼改編而成，讀者在使用時請注明版權。我們使用一個節點（即X取7）為例，代碼如下：
encode study,gen(studynum)
gen double logrr=log(adjrr)
gen double logub=log(ub)
gen double loglb=log(lb)
gen double se=(logub-loglb)/(2*invnormal(0.975) )
mkspline linsp_dose1 7 linsp_dose2 = dose
reglogrrlinsp*
predicty_linsp
lincom linsp_dose1*-1,eform
lincom linsp_dose2*1,eform
其中，does表示自變量。由上述描述可知，1個節點可產生兩個亞分段（dose1和dose2）。命令的前5行是對效應量進行對數轉換，第6行設置節點，第7和8行進行回歸分析及線性預測，最后兩行即表示不同分段下計算的斜率的對數返回值。該命令的關鍵也在于最后4行。
3 實例演示及比較
我們以Stata軟件中自帶數據為例，比較兩種線性模型方法擬合的線性效果。該數據中，自變量為路程（英里，mpg），因變量為價格（price），數據調出命令如下：
sysuse auto,clear
進行普通線性回歸分析，并繪制線性趨勢圖：
reg price mpg
predicty_lin
twoway scatter price mpg || line y_lin mpg
再以分段線性樣條函數再次分析該數據，并繪制趨勢圖，選取一個節點，即以15為分界點：
mkspline linsp_mpg1 15 linsp_mpg2= mpg
reg price linsp*
predicty_linsp
twoway scatter price mpg || line y_linsp mpg,sort clstyle(solid)
圖 1給出兩種方法的結果及線性預測趨勢圖，可以看出，分段樣條線性函數法擬合的線性趨勢更為接近觀察值的分布。
圖1
普通線性回歸模型與分段樣條線性回歸模型結果及趨勢圖的比較
4 總結
本文介紹并比較了傳統線性劑量-反應Meta分析模型及改進的分段樣條線性模型，通過實例演示可知道分段樣條線性回歸模型在擬合線性趨勢時更為靈活和準確。這種改進方法對劑量-反應Meta分析模型是一種完善。并且通過這種方法，可以拓寬線性關系的適用范圍。
在實踐中，在檢測出自變量與因變量之間為線性相關的可能性較小時，往往直接使用非線性模型。這種處理方式忽略了潛在的分段線性分布的情況。在此種情況下，是否需要先使用簡單的分段線性模型目前尚無定論。非線性模型在逼近局部變動上有較大優勢。但這種優勢在特定情況下可能導致過度擬合。我們建議可先使用分段線性模型，若該模型下的結果與非線性模型結果類似，則優先考慮結構和方法學更為簡單的分段線性模型，反之，則使用非線性模型。但此方法是否有效，仍需進一步研究驗證。
總之，在擬合線性關系時，分段樣條線性回歸模型能更好地逼近自變量與因變量之間的因果關系。
在探討自變量與因變量之間因果關系時，通常可以建立一個以自變量為X，應變量為Y的回歸方程 [1]。通過建立的回歸模型，我們可以用觀察值估計線性或非線性趨勢。劑量-反應Meta分析模型也是一種回歸方程，既包含線性模型也包含非線性模型。線性模型即一次函數模型，反應自變量與因變量之間的線性關系；通過對線性模型引入二次項和（或）三次項甚至高次項，可用于擬合自變量與因變量之間的非線性關系 [2]。在循證實踐中，這些非線性模型在擬合非線性趨勢時往往能取得較好的逼近效果，并且可以通過改變模型中的函數項的次數或者通過插值等方法對非線性趨勢進行調整，能適用于大多數情況。目前劑量-反應Meta分析中的線性模型對線性趨勢預測效果則欠佳，特別是當自變量與因變量呈J-型、U-型、S-型、V-型等關系時，其趨勢并非單一走向，傳統線性模型無法準確逼近。本文引入一類新的線性模型——分段線性模型，可以很好地解決這個問題。以下從基礎知識、模型、方法學三方面進行探討，并輔以實例演示，與傳統線性Meta分析模型進行比較。
1 傳統線性模型
1.1 模型表達式
我們在《劑量-反應Meta分析方法學概論》一文 [3]中已對線性模型進行了描述，即引入一個線性回歸函數來反應暴露水平（自變量）與結局指標（因變量）發生風險的線性關系，表達式如下：
| ${{y}_{j}}={{\beta }_{1}}{{x}_{j}}+{{\varepsilon }_{j~}}$ |
yj表示在結局指標發生風險在第j層暴露水平的取值，期望值為β1xj；xj則對應第j層暴露水平的取值。β1是待估參數，也即回歸直線的斜率。εj表示yj的隨機誤差項，即標準誤，滿足均值E(ε)=0，值得注意的是，其方差不等且協方差不為0 [4]。該函數被強制通過原點，也就是在參考劑量暴露水平時對應零風險。繼續以i表示納入研究，Meta分析模型則可表示為：
| ${{y}_{ij}}={{\beta }_{i1}}{{x}_{ij}}+{{\varepsilon }_{ij}}$ |
因為該模型并不完全滿足回歸模型的基本假定，因此參數估計不能使用普通最小二乘法，而需要通過廣義最小二乘法 [4]方可計算單個研究中的回歸斜率及其置信區間。因變量的隨機誤差可通過《劑量-反應Meta分析之限制性立方樣條函數在模型中的應用》一文 [5]中給出的公式求出。在此基礎之上，對單個研究的斜率進行加權合并，得到合并后的回歸斜率?β，如下：
| $\hat{\beta }\text{= }\frac{{{\sum{_{\text{ij}}{{\text{x}}_{\text{ij}}}{{\text{y}}_{\text{ij}}}\text{/}{{\varepsilon }_{\text{ij}}}}}^{2}}}{\sum{_{\text{ij}}{{\text{x}}_{\text{ij}}}^{2}\text{/}{{\varepsilon }_{\text{ij}}}^{2}}}$ |
| ${{\text{S}}_{\beta }}^{\text{2}}\text{= }\frac{\text{1}}{\sum{_{\text{ij}}{{\text{x}}_{\text{ij}}}^{2}\text{/}{{\varepsilon }_{\text{ij}}}^{2}}}$ |
求出斜率后，就可以知道yij的期望值，實際上，我們最終感興趣的結果也即這個期望值。該公式也可以簡化為矩陣形式寫出：
| $\hat{\beta }={{\left( X'{{V}^{-1}}X \right)}^{-1}}X'{{V}^{-1}}Y$ |
| $S{{~}_{\beta }}^{2}={{\left( X'V-1X \right)}^{-1}}$ |
1.2 模型的意義
在循證實踐中，單個研究的暴露水平分布往往不同，如最低暴露劑量不一致，導致函數中暴露基線差異大。為避免此種情況的發生，在進行線性劑量-反應Meta分析時，通常以合并后的斜率來表示自變量與因變量之間的關系。以公式（2）為例，x每變化一個單位，對應y則變化個單位。當然，這種變化可以顯著也可以不顯著，可根據假設檢驗來判斷：H0：=0，計算該零假設成立的情況下，觀察結果在樣本分布中可能發生的概率 [1]。
1.3 操作代碼
多種軟件可以進行線性劑量反應Meta分析，我們以Stata/SE12.0軟件為例，實現代碼如下：
encode study,gen(studynum)
gen double logrr=log(adjrr)
gen double logub=log(ub)
gen double loglb=log(lb)
gen double se=(logub-loglb)/(2*invnormal(0.975) )
glstlogrr dose,se(se) cov(n case) pfirst(id studynum) ts(r)
lincom dose*1,eform
關于代碼的意義及解釋請參考《應用STATA做Meta分析》一書 [6]相關章節及《應用 R軟件dosresmeta程序包和mvmeta程序包實現劑量-反應關系Meta分析》一文 [2]相關內容。
2 分段線性模型
分段線性模型主要是針對存在分段線性趨勢的分布的情況。模型可以用多種函數表示，為方便讀者理解，我們與系列文章 [5]保持一致，使用分段線性樣條函數 [7]。也是對經典限制性立方樣條函數法下的劑量-反應Meta分析模型缺陷的一種補充。值得注意的是，該函數并非本團隊原創，我們的工作只是將該函數引入劑量-反應Meta分析，并研發出了相關代碼。
2.1 模型的表達式
跟限制性立方樣條函數類似，直線樣條函數也是通過節點連接的分段直線函數，需滿足每個節點處連續。假設以ki表示節點（i=1，2，…，n-1），以單個研究為例，模型表達式如下：
| $y=g\left( x \right)={{b}_{1i}}{{X}_{i}}+\varepsilon $ |
式中，b1k為回歸斜率，并且隨著節點不同而改變，Xi是X的函數，意義如下：
| ${{X}_{1}}=X$ |
| ${{X}_{i+1}}=max\left( X-{{k}_{i-1}} \right)$ |
X的值大于或等于ki-1，公式（9）為X-ki-1，反之，則為0。為便于理解，我們將公式進行變換，形成一個類似多元線性回歸的模型：
| $y=g\left( x \right)={{b}_{1}}{{X}_{1}}+{{b}_{1i+1}}{{X}_{i+1}}+\varepsilon ~$ |
以j（j=1，2，3…）表示單個研究，Meta分析模型則如下：
| ${{y}_{j}}={{\beta }_{1ij}}{{x}_{ij}}+{{\varepsilon }_{j}}$ |
2.2 模型的意義
分段線性函數模型的意義與普通線性函數類似，不同的是通過不同節點將該函數分成了多段連續的線性函數，這些由節點分成的子區間的線性函數，斜率可同可不同，因此可更靈活地擬合線性趨勢。
仍然以X的分布上10th（k1）、50th（k2）、90th（k3）這3個節點為例，3個節點將整個趨勢分為了4個小段，即0th～10th、10th～50th、50th～90th、90th～100th。函數表達式如下：
y=g(x)=b1X1+ε；0th≤X1＜10th
y=g(x)=b12X2+ε；10th≤X1＜50th
y=g(x)=b13X3+ε；50th≤X1＜90th
y=g(x)=b14X4+ε；90th≤X1＜100th
根據廣義最小二乘法逐一計算出不同分段的斜率，即可反應自變量與因變量之間的分段線性關系。即在每一個分段內，x每變化一個單位，對應y則變化b1i個單位。
2.3 操作代碼
本文給出主要操作代碼，值得注意的是，該代碼由本研究團隊根據簡單回歸模型、多重回歸模型及非線性Meta分析模型的代碼改編而成，讀者在使用時請注明版權。我們使用一個節點（即X取7）為例，代碼如下：
encode study,gen(studynum)
gen double logrr=log(adjrr)
gen double logub=log(ub)
gen double loglb=log(lb)
gen double se=(logub-loglb)/(2*invnormal(0.975) )
mkspline linsp_dose1 7 linsp_dose2 = dose
reglogrrlinsp*
predicty_linsp
lincom linsp_dose1*-1,eform
lincom linsp_dose2*1,eform
其中，does表示自變量。由上述描述可知，1個節點可產生兩個亞分段（dose1和dose2）。命令的前5行是對效應量進行對數轉換，第6行設置節點，第7和8行進行回歸分析及線性預測，最后兩行即表示不同分段下計算的斜率的對數返回值。該命令的關鍵也在于最后4行。
3 實例演示及比較
我們以Stata軟件中自帶數據為例，比較兩種線性模型方法擬合的線性效果。該數據中，自變量為路程（英里，mpg），因變量為價格（price），數據調出命令如下：
sysuse auto,clear
進行普通線性回歸分析，并繪制線性趨勢圖：
reg price mpg
predicty_lin
twoway scatter price mpg || line y_lin mpg
再以分段線性樣條函數再次分析該數據，并繪制趨勢圖，選取一個節點，即以15為分界點：
mkspline linsp_mpg1 15 linsp_mpg2= mpg
reg price linsp*
predicty_linsp
twoway scatter price mpg || line y_linsp mpg,sort clstyle(solid)
圖 1給出兩種方法的結果及線性預測趨勢圖，可以看出，分段樣條線性函數法擬合的線性趨勢更為接近觀察值的分布。
圖1
普通線性回歸模型與分段樣條線性回歸模型結果及趨勢圖的比較
4 總結
本文介紹并比較了傳統線性劑量-反應Meta分析模型及改進的分段樣條線性模型，通過實例演示可知道分段樣條線性回歸模型在擬合線性趨勢時更為靈活和準確。這種改進方法對劑量-反應Meta分析模型是一種完善。并且通過這種方法，可以拓寬線性關系的適用范圍。
在實踐中，在檢測出自變量與因變量之間為線性相關的可能性較小時，往往直接使用非線性模型。這種處理方式忽略了潛在的分段線性分布的情況。在此種情況下，是否需要先使用簡單的分段線性模型目前尚無定論。非線性模型在逼近局部變動上有較大優勢。但這種優勢在特定情況下可能導致過度擬合。我們建議可先使用分段線性模型，若該模型下的結果與非線性模型結果類似，則優先考慮結構和方法學更為簡單的分段線性模型，反之，則使用非線性模型。但此方法是否有效，仍需進一步研究驗證。
總之，在擬合線性關系時，分段樣條線性回歸模型能更好地逼近自變量與因變量之間的因果關系。
