限制性立方樣條是趨勢逼近中理想的函數模型，也是劑量-反應Meta分析領域常用的一種模型。它是基于參數法下的函數，實質是一條各節點處光滑的分段多項式，并且在定義域每一個子區間中都是三次多項式。限制性立方樣條函數對非線性趨勢的逼近效果比較好，并且可以通過改變節點個數和（或）位置來調整所逼近的曲線。前篇文章我們介紹并歸納了常用參數法線性及非線性劑量-反應Meta分析模型，本篇我們將從模型的構建、參數的合并、節點選取等方面詳細探討限制性立方樣條函數模型的方法學。
引用本文： 黃靜宇, 張超, 李勝, 鄺心穎, 牛玉明, 徐暢. 劑量-反應Meta分析之限制性立方樣條函數的應用. 中國循證醫學雜志, 2015, 15(12): 1471-1474. doi: 10.7507/1672-2531.20150239 復制
劑量-反應Meta分析是一種研究暴露（干預）水平變化與結局指標發生風險之間關系的一種Meta分析方法，在觀察性研究領域應用十分廣泛 [1]。基于比較完善的函數模型及成熟的參數估計方法，優質劑量-反應Meta分析能為循證實踐及決策提供高級別證據 [2]。劑量-反應Meta分析中，自變量與應變量之間通常有兩種關系，其一是線性關系，其二則是非線性關系。當然，亦存在一種介于線性與非線性之間的關系，即邊緣非線性關系。線性關系通過一次回歸函數較易擬合，劑量的核心問題在于非線性趨勢的逼近。目前循證實踐中應用最為廣泛的逼近模型為限制性立方樣條函數，它是基于參數法下的函數，實質是一條各節點處光滑的分段多項式，并且它在定義域每一個“非限制”子區間中都是三次多項式 [3]。限制性立方樣條函數對非線性趨勢的逼近效果比較理想，并且可以通過改變節點個數和（或）位置來調整所逼近的曲線，適用于絕大多數的非線性關系。本文將從函數模型的特點與缺陷等方面對其進行深入介紹及探討。
1 基礎知識
限制性立方樣條函數用于劑量-反應Meta分析模型的構建是基于“兩步法”思想，即先擬合每篇原始研究的立方樣條函數，然后將這些函數的參數進行加權合并 [4]。其關鍵在于參數研究間的協方差的估算，解決此問題通常用到的方法為最小二乘法 [4]。為方便后續說明，我們以一個對數線性模型及表格數據為例（表 1）。
| $Y=bX+\varepsilon $ |
Y表示非參照組的效應量對數矢量值，X表示暴露水平，即表一中的A，ε表示研究內協方差矩陣。該模型即構成最簡單的劑量-反應關系模型，并被強制通過原點。通過計算該函數的參數及方差協方差，即可反映自變量與應變量的關系，也即暴露水平和結局事件發生風險的關系。
表1
劑量-反應Meta分析數據類型
以i代表研究，j代表研究內各層暴露水平。假設bi（i=1，2…）為單篇研究的斜率，kn為節點（n為節點個數），Sj為研究內協方差矩陣。顯然，該矩陣為對稱矩陣，并以各層方差Varj（i=1，2…）的矢量構成對角線上的元素。引入研究內協方差是由于各層暴露水平的方差并不獨立 [4]。根據上述表格，Nj和Mj的期望值可由未校正的四格表的值和校正后的效應量RR進行估算，如下：
| $({{N}_{j}}\times {{M}_{0}}/{{N}_{0}}\times {{M}_{j}})=exp(InR{{R}_{j}})~$ |
| $({{N}_{j}}\times ({{M}_{0}}+{{N}_{0}})\text{ }/{{N}_{0}}\times ({{M}_{j}}+{{N}_{j}}))=exp(InR{{R}_{j}})$ |
其中，公式（1）適用于病例-對照研究或隨機對照研究，公式（2）則適用于隊列研究。需要注意的是，式中的RR為廣義RR，可代表OR、RR及HR [1]。同上順序，參照劑量的協方差計算如下：
| ${{S}_{0}}=(1/{{N}_{0}}+1/{{M}_{0}})~$ |
| ${{S}_{0}}=(1/{{N}_{0}}-1/({{M}_{0}}+{{N}_{0}})\text{ })$ |
Nj和Mj的期望值以及參照劑量協方差計算出來后，可根據他們進行研究內協方差的估算，按照公式（1）和（2）的順序如下：
| ${{S}_{j}}=\sqrt{1/{{N}_{j}}+1/{{M}_{j}}+1/{{N}_{0}}+1/{{M}_{0}}}$ |
| ${{S}_{j}}=\sqrt{\left( 1/{{N}_{0}}-1/\left( {{M}_{0}}+{{N}_{0}} \right) \right)+\left( 1/{{N}_{j}}-1\left( {{M}_{j}}+{{N}_{j}} \right) \right)}$ |
方差之間相關系數rab（a、b∈j）計算公式入下：
| ${{r}_{ab}}={{S}_{0}}/{{S}_{a}}{{S}_{b}}$ |
通過以上值可求出協方差矩陣內任對角元素vab：
| ${{v}_{ab}}={{r}_{ab}}\sqrt{{{v}_{a}}{{v}_{b}}}$ |
矩陣對角線的元素為表一中的SEj。每篇研究內擬合的回歸模型的斜率則可根據最小二乘法 [5]及上述值進行計算：
| $\hat{b}={{S}_{j}}{{{{X}'}}_{i}}\sum{_{j}^{-1}}\ln R{{R}_{j}}$ |
2 限制性立方樣條回歸模型
2.1 模型的建立及意義
非線性關系的擬合思路是在線性關系上加入高次項，如二次項和（或）三次項，這樣構成了高次多項式模型。限制性立方樣條函數屬于多項式中的一種，最大的特點在于它進行了樣條插值并對趨勢首尾兩端進行了線性限制，即在第一個節點前的趨勢和最后一個節點后的趨勢強制限制為線性。如下給出非限制性立方樣條函數的模型 [6]：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}X+{{b}_{2}}{{X}_{2}}+{{b}_{3}}{{X}_{3}}+{{\sum }_{i=1}}^{n}{{b}_{3+i}}~max(X-{{k}_{i}},0){{~}^{3}}+{{\varepsilon }_{i}}$ |
| $f\text{ }X{{k}_{i}},max(X-{{k}_{i}},0){{~}^{3}}={{(x-{{k}_{i}})}^{3}}~$ |
| $f\text{ }X\le {{k}_{i}},max(X-{{k}_{i}},0){{~}^{3}}=0$ |
∑i=1n b3+i max(X-ki,0)3，i=1，…，n，表示ki為取自X的分布節點，若X的值大于或等于ki，此項為∑i=1nβ3+i (X-ki)3，相反，則此項為0。對該函數左端限制為線性，意味著將b2和b3強制設置為0，也即去掉2次和3次項，函數變為：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}X+{{\sum }_{i=1}}^{n}{{b}_{1+i}}~max(X-{{k}_{i}},0){{~}^{3}}+{{\varepsilon }_{i}}~$ |
對原函數進行右端限制跟左端限制處理方法類似，僅需將左端限制中的X換成-X，模型如下：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}\left( -X \right)+{{\sum }_{i=1}}^{n}{{b}_{1+i}}~max(X-{{k}_{i}},0){{~}^{3}}+{{\varepsilon }_{i}}$ |
對原函數進行雙側限制后，產生限制性立方樣條函數，該函數有n-1個回歸參數，將max(-ki,0)3以μi表示：
| $g\left( X \right)={{b}_{0}}+{{\sum }_{i=1}}^{n-1}{{b}_{1+i}}{{X}_{i}}+{{\varepsilon }_{i}}$ |
| ${{X}_{i}}=\left\{ {{\mu }_{i-1}}-{{\mu }_{n-1}}\left( \frac{\left( {{k}_{n}}-{{k}_{i-1}} \right)}{\left( {{k}_{n}}-{{k}_{n-1}} \right)} \right)+{{\mu }_{n}}\left( \frac{\left( {{k}_{n}}-{{k}_{i-1}} \right)}{\left( {{k}_{n}}-{{k}_{n-1}} \right)} \right) \right\}/{{\left( {{k}_{n}}-{{k}_{1}} \right)}^{2}}$ |
樣條中（公式15），X1表示原始暴露變量X，余下Xi為X的函數（公式16）。若X≤ki,公式（16）為0，若ki＜X≤kn-1,則Xi=(X-ki)3，若kn-1＜X≤kn，Xi有兩個立方項，即公式（16）中存在μi-1和μn-1的兩項，若kn＜X，公式（16）有3個樣條項，但3個立方項之和為0 [4]。從公式及其含義中不難看出，第一個節點前和最后一個節點后的趨勢均被限制為了線性。
以X的分布上10th（k1）、50th（k2）、90th（k3）3個節點為例，3個節點將整體趨勢分為了4小段，即0th～10th、10th～50th、50th～90th、90th～100th。第一段內的X≤ k1，由公式（11）和（16）可知，∑i=1n-2b1+i Xi= 0，公式（15）變為：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}{{X}_{1}}+{{\varepsilon }_{i}}~$ |
即整體趨勢的第一段為線性，最后一段同樣為線性，中間兩段仍為立方樣條函數。限制的目的是為了防止四舍五入法平方或立方后產生的誤差。通過第一部分參數估計方法計算出函數中的待估參數，即可反應自變量與應變量之間的因果關系。
2.2 非線性檢驗
在確定模型后，我們需要通過觀察值判斷所擬合的趨勢是線性還是非線性。仍以3個節點為例，產生含兩個X的項，公式（15）則變為：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}X+{{b}_{2}}{{X}_{2}}+{{\varepsilon }_{i}}$ |
判斷整體趨勢是否為非線性，我們只需要考慮函數的非線性部分，即b2X2。由此可以通過兩種方法進行檢驗，其一為似然性檢驗，其二為卡方檢驗（Wald檢驗）。第一種方法分別計算非線性的和線性情況下的似然性，并比較似然性大小，若非線性的似然性大，則認為該趨勢為非線性，反之亦然 [4]。第二種方法則需進行假設，我們不妨假設第二條樣條函數的參數b2為0，即函數非線性部分為0，并計算它取0值的可能性，如果可能性較小（通常取P＜0.05），則可認為該函數為非線性；若取0值的可能性較大，則還不能認為它是非線性 [7]。
3 劑量-反應Meta分析的應用
進行劑量-反應Meta分析，需要將上述步驟中單篇研究的參數求出，然后通過隨機或者固定或混合效應模型進行加權合并即可，以固定效應模型為例，加權平均后的參數的計算公式為：
| $\hat{\beta }=\sum{_{1}^{i}}{{W}_{i}}{{b}_{i}}/\sum{{{W}_{i}}}$ |
其中Wi表示第i個研究中的權重。通過第一部分給出的方差及方差之間的協方差計算公式，我們可以計算出第i個研究中回歸系數bi。值得注意的是，在Meta分析中，方差-協方差矩陣的對角元素（即方差）為單個研究中的方差的集合。如納入10個研究，每個研究有5個不同暴露水平（對應5個方差），則此Meta分析的方差-協方差矩陣對角元素總共有50個。以矩陣公式表示入下：
| $\left[ \begin{matrix} {{s}_{11}} & {{s}_{11-12}} & \cdots & {{s}_{11-ij}} \\ {{s}_{11-12}} & {{s}_{12}} & \cdots & {{s}_{12-ij}} \\ \vdots & \vdots & \vdots & \vdots \\ {{s}_{11-ij}} & {{s}_{12-ij}} & \cdots & {{s}_{ij}} \\ \end{matrix} \right]$ |
s11-ij表示第一個研究的參照劑量下RR的方差與第i個研究內第j層暴露劑量對應的RR的方差的協方差。該矩陣為對稱矩陣，即以對角線（方差）為中軸，左下角與右上角的所有元素（協方差）相等。
3.1 軟件
基于限制性立方樣條函數的劑量-反應Meta分析有多種軟件可以實現其運算，創始人Orsini N提供了至少3種軟件的程序，分別為Stata軟件、SAS軟件和R軟件。軟件的操作較為簡單，只需安裝特定程序包并按照他們提供的例題進行操作即可完成。以Stata軟件為例，所需安裝的程序為glst和xblc兩個，安裝命令如下：
ssc install glst
findit xblc
余下操作（包括數據錄入形式）請參考Orsini N [7]的原文，我們后續文章也會對其進行詳細介紹。
3.2 注意事項
3.2.1 節點問題
由于限制性立方樣條函數的靈活性，對非線性逼近效果非常理想，通常情況下三次樣條函數即可滿足大多數非線性趨勢的逼近。對于實際應用中的節點取值問題，目前仍無強有力的定論。理論上講，節點個數越多，越接近真實趨勢，更能反映局部變動情況，但節點個數多時也可能導致過度逼近；節點個數少時，曲線平滑，更能反映整體的趨勢走向。通常認為，3到7個節點是比較合理的。值得注意的是，如果樣本量較大，建議增加節點個數，而小樣本量則建議盡可能用3個節點以保證結果的統計效能 [4]。
3.2.2 線性關系問題
基于限制性立方樣條函數的劑量-反應Meta分析對非線性的逼近較理想，但對線性關系的反應仍有一定缺陷，如“V”型關系是個典型的分段線性關系，此時則無法使用普通線性回歸進行逼近，而需要使用更靈活的分段線性回歸。對于分段線性回歸在劑量-反應Meta分析中的應用，我們會在后續文章中詳細介紹。
4 總結
本文詳細介紹并分析了限制性立方樣條函數在劑量-反應Meta分析模型中應用的理論及相關問題，而實際上這些理論及問題在循證實踐中也十分常見。在進行一項劑量-反應Meta分析時，我們需要在方法學部分對該函數進行一定的描述，比如，如何進行非線性檢驗，納入研究對暴露水平的分層要求（至少3層），以及選取幾個節點進行趨勢逼近等等 [8, 9]。本文對限制性立方樣條函數給出了較為系統和詳細的說明，對不熟悉或無法理解該函數的循證醫學同行都會有一定幫助。
本文同時也介紹了限制性立方樣條函數在劑量-反應Meta分析中需注意的問題。盡管它適用于大多數分布的趨勢逼近，但在一些特殊情況下仍有一定缺陷，如我們在第三部分討論的線性關系問題中存在的缺陷。盡管如此，限制性立方樣條函數模型仍是目前使用最廣泛的擬合劑量-反應Meta分析的方法 [8]。
當然，能應用于劑量-反應Meta分析的函數模型遠不止限制性立方樣條函數，其他各種分段多項式函數均可替代樣條函數，如靈活分段多項式。甚至已有統計學家開發出了基于半參數法下的逼近函數。但是這些函數也或多或少存在一些缺陷，而導致他們的應用并不廣泛。劑量-反應Meta分析模型仍存在許多共同缺陷，最迫切需要解決的問題之一即原始數據的完整性，許多原始研究并未給出劑量-反應Meta分析需要的數據（表 1），因此多數情況下不得不對這些缺失值進行估算，使結果可信度降低。這些問題仍需要進一步的深入研究來解決。
劑量-反應Meta分析是一種研究暴露（干預）水平變化與結局指標發生風險之間關系的一種Meta分析方法，在觀察性研究領域應用十分廣泛 [1]。基于比較完善的函數模型及成熟的參數估計方法，優質劑量-反應Meta分析能為循證實踐及決策提供高級別證據 [2]。劑量-反應Meta分析中，自變量與應變量之間通常有兩種關系，其一是線性關系，其二則是非線性關系。當然，亦存在一種介于線性與非線性之間的關系，即邊緣非線性關系。線性關系通過一次回歸函數較易擬合，劑量的核心問題在于非線性趨勢的逼近。目前循證實踐中應用最為廣泛的逼近模型為限制性立方樣條函數，它是基于參數法下的函數，實質是一條各節點處光滑的分段多項式，并且它在定義域每一個“非限制”子區間中都是三次多項式 [3]。限制性立方樣條函數對非線性趨勢的逼近效果比較理想，并且可以通過改變節點個數和（或）位置來調整所逼近的曲線，適用于絕大多數的非線性關系。本文將從函數模型的特點與缺陷等方面對其進行深入介紹及探討。
1 基礎知識
限制性立方樣條函數用于劑量-反應Meta分析模型的構建是基于“兩步法”思想，即先擬合每篇原始研究的立方樣條函數，然后將這些函數的參數進行加權合并 [4]。其關鍵在于參數研究間的協方差的估算，解決此問題通常用到的方法為最小二乘法 [4]。為方便后續說明，我們以一個對數線性模型及表格數據為例（表 1）。
| $Y=bX+\varepsilon $ |
Y表示非參照組的效應量對數矢量值，X表示暴露水平，即表一中的A，ε表示研究內協方差矩陣。該模型即構成最簡單的劑量-反應關系模型，并被強制通過原點。通過計算該函數的參數及方差協方差，即可反映自變量與應變量的關系，也即暴露水平和結局事件發生風險的關系。
表1
劑量-反應Meta分析數據類型
以i代表研究，j代表研究內各層暴露水平。假設bi（i=1，2…）為單篇研究的斜率，kn為節點（n為節點個數），Sj為研究內協方差矩陣。顯然，該矩陣為對稱矩陣，并以各層方差Varj（i=1，2…）的矢量構成對角線上的元素。引入研究內協方差是由于各層暴露水平的方差并不獨立 [4]。根據上述表格，Nj和Mj的期望值可由未校正的四格表的值和校正后的效應量RR進行估算，如下：
| $({{N}_{j}}\times {{M}_{0}}/{{N}_{0}}\times {{M}_{j}})=exp(InR{{R}_{j}})~$ |
| $({{N}_{j}}\times ({{M}_{0}}+{{N}_{0}})\text{ }/{{N}_{0}}\times ({{M}_{j}}+{{N}_{j}}))=exp(InR{{R}_{j}})$ |
其中，公式（1）適用于病例-對照研究或隨機對照研究，公式（2）則適用于隊列研究。需要注意的是，式中的RR為廣義RR，可代表OR、RR及HR [1]。同上順序，參照劑量的協方差計算如下：
| ${{S}_{0}}=(1/{{N}_{0}}+1/{{M}_{0}})~$ |
| ${{S}_{0}}=(1/{{N}_{0}}-1/({{M}_{0}}+{{N}_{0}})\text{ })$ |
Nj和Mj的期望值以及參照劑量協方差計算出來后，可根據他們進行研究內協方差的估算，按照公式（1）和（2）的順序如下：
| ${{S}_{j}}=\sqrt{1/{{N}_{j}}+1/{{M}_{j}}+1/{{N}_{0}}+1/{{M}_{0}}}$ |
| ${{S}_{j}}=\sqrt{\left( 1/{{N}_{0}}-1/\left( {{M}_{0}}+{{N}_{0}} \right) \right)+\left( 1/{{N}_{j}}-1\left( {{M}_{j}}+{{N}_{j}} \right) \right)}$ |
方差之間相關系數rab（a、b∈j）計算公式入下：
| ${{r}_{ab}}={{S}_{0}}/{{S}_{a}}{{S}_{b}}$ |
通過以上值可求出協方差矩陣內任對角元素vab：
| ${{v}_{ab}}={{r}_{ab}}\sqrt{{{v}_{a}}{{v}_{b}}}$ |
矩陣對角線的元素為表一中的SEj。每篇研究內擬合的回歸模型的斜率則可根據最小二乘法 [5]及上述值進行計算：
| $\hat{b}={{S}_{j}}{{{{X}'}}_{i}}\sum{_{j}^{-1}}\ln R{{R}_{j}}$ |
2 限制性立方樣條回歸模型
2.1 模型的建立及意義
非線性關系的擬合思路是在線性關系上加入高次項，如二次項和（或）三次項，這樣構成了高次多項式模型。限制性立方樣條函數屬于多項式中的一種，最大的特點在于它進行了樣條插值并對趨勢首尾兩端進行了線性限制，即在第一個節點前的趨勢和最后一個節點后的趨勢強制限制為線性。如下給出非限制性立方樣條函數的模型 [6]：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}X+{{b}_{2}}{{X}_{2}}+{{b}_{3}}{{X}_{3}}+{{\sum }_{i=1}}^{n}{{b}_{3+i}}~max(X-{{k}_{i}},0){{~}^{3}}+{{\varepsilon }_{i}}$ |
| $f\text{ }X{{k}_{i}},max(X-{{k}_{i}},0){{~}^{3}}={{(x-{{k}_{i}})}^{3}}~$ |
| $f\text{ }X\le {{k}_{i}},max(X-{{k}_{i}},0){{~}^{3}}=0$ |
∑i=1n b3+i max(X-ki,0)3，i=1，…，n，表示ki為取自X的分布節點，若X的值大于或等于ki，此項為∑i=1nβ3+i (X-ki)3，相反，則此項為0。對該函數左端限制為線性，意味著將b2和b3強制設置為0，也即去掉2次和3次項，函數變為：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}X+{{\sum }_{i=1}}^{n}{{b}_{1+i}}~max(X-{{k}_{i}},0){{~}^{3}}+{{\varepsilon }_{i}}~$ |
對原函數進行右端限制跟左端限制處理方法類似，僅需將左端限制中的X換成-X，模型如下：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}\left( -X \right)+{{\sum }_{i=1}}^{n}{{b}_{1+i}}~max(X-{{k}_{i}},0){{~}^{3}}+{{\varepsilon }_{i}}$ |
對原函數進行雙側限制后，產生限制性立方樣條函數，該函數有n-1個回歸參數，將max(-ki,0)3以μi表示：
| $g\left( X \right)={{b}_{0}}+{{\sum }_{i=1}}^{n-1}{{b}_{1+i}}{{X}_{i}}+{{\varepsilon }_{i}}$ |
| ${{X}_{i}}=\left\{ {{\mu }_{i-1}}-{{\mu }_{n-1}}\left( \frac{\left( {{k}_{n}}-{{k}_{i-1}} \right)}{\left( {{k}_{n}}-{{k}_{n-1}} \right)} \right)+{{\mu }_{n}}\left( \frac{\left( {{k}_{n}}-{{k}_{i-1}} \right)}{\left( {{k}_{n}}-{{k}_{n-1}} \right)} \right) \right\}/{{\left( {{k}_{n}}-{{k}_{1}} \right)}^{2}}$ |
樣條中（公式15），X1表示原始暴露變量X，余下Xi為X的函數（公式16）。若X≤ki,公式（16）為0，若ki＜X≤kn-1,則Xi=(X-ki)3，若kn-1＜X≤kn，Xi有兩個立方項，即公式（16）中存在μi-1和μn-1的兩項，若kn＜X，公式（16）有3個樣條項，但3個立方項之和為0 [4]。從公式及其含義中不難看出，第一個節點前和最后一個節點后的趨勢均被限制為了線性。
以X的分布上10th（k1）、50th（k2）、90th（k3）3個節點為例，3個節點將整體趨勢分為了4小段，即0th～10th、10th～50th、50th～90th、90th～100th。第一段內的X≤ k1，由公式（11）和（16）可知，∑i=1n-2b1+i Xi= 0，公式（15）變為：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}{{X}_{1}}+{{\varepsilon }_{i}}~$ |
即整體趨勢的第一段為線性，最后一段同樣為線性，中間兩段仍為立方樣條函數。限制的目的是為了防止四舍五入法平方或立方后產生的誤差。通過第一部分參數估計方法計算出函數中的待估參數，即可反應自變量與應變量之間的因果關系。
2.2 非線性檢驗
在確定模型后，我們需要通過觀察值判斷所擬合的趨勢是線性還是非線性。仍以3個節點為例，產生含兩個X的項，公式（15）則變為：
| $g\left( X \right)={{b}_{0}}+{{b}_{1}}X+{{b}_{2}}{{X}_{2}}+{{\varepsilon }_{i}}$ |
判斷整體趨勢是否為非線性，我們只需要考慮函數的非線性部分，即b2X2。由此可以通過兩種方法進行檢驗，其一為似然性檢驗，其二為卡方檢驗（Wald檢驗）。第一種方法分別計算非線性的和線性情況下的似然性，并比較似然性大小，若非線性的似然性大，則認為該趨勢為非線性，反之亦然 [4]。第二種方法則需進行假設，我們不妨假設第二條樣條函數的參數b2為0，即函數非線性部分為0，并計算它取0值的可能性，如果可能性較小（通常取P＜0.05），則可認為該函數為非線性；若取0值的可能性較大，則還不能認為它是非線性 [7]。
3 劑量-反應Meta分析的應用
進行劑量-反應Meta分析，需要將上述步驟中單篇研究的參數求出，然后通過隨機或者固定或混合效應模型進行加權合并即可，以固定效應模型為例，加權平均后的參數的計算公式為：
| $\hat{\beta }=\sum{_{1}^{i}}{{W}_{i}}{{b}_{i}}/\sum{{{W}_{i}}}$ |
其中Wi表示第i個研究中的權重。通過第一部分給出的方差及方差之間的協方差計算公式，我們可以計算出第i個研究中回歸系數bi。值得注意的是，在Meta分析中，方差-協方差矩陣的對角元素（即方差）為單個研究中的方差的集合。如納入10個研究，每個研究有5個不同暴露水平（對應5個方差），則此Meta分析的方差-協方差矩陣對角元素總共有50個。以矩陣公式表示入下：
| $\left[ \begin{matrix} {{s}_{11}} & {{s}_{11-12}} & \cdots & {{s}_{11-ij}} \\ {{s}_{11-12}} & {{s}_{12}} & \cdots & {{s}_{12-ij}} \\ \vdots & \vdots & \vdots & \vdots \\ {{s}_{11-ij}} & {{s}_{12-ij}} & \cdots & {{s}_{ij}} \\ \end{matrix} \right]$ |
s11-ij表示第一個研究的參照劑量下RR的方差與第i個研究內第j層暴露劑量對應的RR的方差的協方差。該矩陣為對稱矩陣，即以對角線（方差）為中軸，左下角與右上角的所有元素（協方差）相等。
3.1 軟件
基于限制性立方樣條函數的劑量-反應Meta分析有多種軟件可以實現其運算，創始人Orsini N提供了至少3種軟件的程序，分別為Stata軟件、SAS軟件和R軟件。軟件的操作較為簡單，只需安裝特定程序包并按照他們提供的例題進行操作即可完成。以Stata軟件為例，所需安裝的程序為glst和xblc兩個，安裝命令如下：
ssc install glst
findit xblc
余下操作（包括數據錄入形式）請參考Orsini N [7]的原文，我們后續文章也會對其進行詳細介紹。
3.2 注意事項
3.2.1 節點問題
由于限制性立方樣條函數的靈活性，對非線性逼近效果非常理想，通常情況下三次樣條函數即可滿足大多數非線性趨勢的逼近。對于實際應用中的節點取值問題，目前仍無強有力的定論。理論上講，節點個數越多，越接近真實趨勢，更能反映局部變動情況，但節點個數多時也可能導致過度逼近；節點個數少時，曲線平滑，更能反映整體的趨勢走向。通常認為，3到7個節點是比較合理的。值得注意的是，如果樣本量較大，建議增加節點個數，而小樣本量則建議盡可能用3個節點以保證結果的統計效能 [4]。
3.2.2 線性關系問題
基于限制性立方樣條函數的劑量-反應Meta分析對非線性的逼近較理想，但對線性關系的反應仍有一定缺陷，如“V”型關系是個典型的分段線性關系，此時則無法使用普通線性回歸進行逼近，而需要使用更靈活的分段線性回歸。對于分段線性回歸在劑量-反應Meta分析中的應用，我們會在后續文章中詳細介紹。
4 總結
本文詳細介紹并分析了限制性立方樣條函數在劑量-反應Meta分析模型中應用的理論及相關問題，而實際上這些理論及問題在循證實踐中也十分常見。在進行一項劑量-反應Meta分析時，我們需要在方法學部分對該函數進行一定的描述，比如，如何進行非線性檢驗，納入研究對暴露水平的分層要求（至少3層），以及選取幾個節點進行趨勢逼近等等 [8, 9]。本文對限制性立方樣條函數給出了較為系統和詳細的說明，對不熟悉或無法理解該函數的循證醫學同行都會有一定幫助。
本文同時也介紹了限制性立方樣條函數在劑量-反應Meta分析中需注意的問題。盡管它適用于大多數分布的趨勢逼近，但在一些特殊情況下仍有一定缺陷，如我們在第三部分討論的線性關系問題中存在的缺陷。盡管如此，限制性立方樣條函數模型仍是目前使用最廣泛的擬合劑量-反應Meta分析的方法 [8]。
當然，能應用于劑量-反應Meta分析的函數模型遠不止限制性立方樣條函數，其他各種分段多項式函數均可替代樣條函數，如靈活分段多項式。甚至已有統計學家開發出了基于半參數法下的逼近函數。但是這些函數也或多或少存在一些缺陷，而導致他們的應用并不廣泛。劑量-反應Meta分析模型仍存在許多共同缺陷，最迫切需要解決的問題之一即原始數據的完整性，許多原始研究并未給出劑量-反應Meta分析需要的數據（表 1），因此多數情況下不得不對這些缺失值進行估算，使結果可信度降低。這些問題仍需要進一步的深入研究來解決。
