引用本文： 付金玉, 秦超英. Meta分析中各研究偏倚的量化分析. 中國循證醫學雜志, 2016, 16(9): 1112-1116. doi: 10.7507/1672-2531.20160169 復制
Meta分析是匯總具有相同目的的多個獨立研究結果，剖析研究間的差異特征，對研究結果進行定量綜合評價的方法[1]。在傳統的Meta分析隨機效應模型中[2]，沒有考慮偏倚的影響，通常假設納入Meta分析的各個研究的效應量服從正態分布，由此得出的總體效應量只受研究間異質性和研究內方差的影響。實際上，異質性由多種因素造成，偏倚可認為是其中一個因素，然而由于偏倚本身的特殊性，如發表偏倚或研究者的偏好，使觀測效應量普遍偏高或偏低，即不是傳統模型中的正態分布，而是呈現非對稱的偏正態分布趨勢。因此，我們希望把偏倚從異質性中分離出來，單獨加以量化表征，使得模型能夠描述各研究觀測效應量非對稱分布的情形，而偏正態分布模型可以描述這類情形。于是，人們提出了對Meta分析中的偏倚進行量化[3-5]，即假設納入Meta分析中的各個研究的效應量服從相同形狀參數的偏正態分布，通過形狀參數對偏倚進行量化，進而可分析偏倚和異質性對結果的影響。但由于各個研究存在不同，如如何處理試驗過程中參與者的流失等，致使各個研究的偏倚不盡相同。因此本文提出將各個研究的偏倚通過不同的形狀參數進行量化，即假設納入Meta分析中的各個研究的效應量服從不同形狀參數的偏正態分布，進而建立基于形狀參數不同的偏正態分布下的隨機效應模型，然后采用極大似然估計法求出每個研究偏倚的估計值，進而可以分析各個研究的偏倚和研究間的異質性對總體效應量的影響。同時，本文還通過模擬計算和實例列舉，驗證本文提出的模型和算法的可行性和有效性。
1 基于偏正態的隨機效應模型
設效應量yi服從形狀參數相同的偏正態分布[6]隨機效應模型[5]為：
| $ {y_i}=\mu+\tau {\xi _i}+{\varepsilon _i}, i=1, \ldots, n $ |
其中yi，i=1, …, n是n個相互獨立的觀測效應量，n表示納入研究數。ξi~SN(0, 1, γ)，εi~N(0, δi2)且ξi和εi相互獨立，γ為形狀參數，表示每個研究的偏倚參數；δi2為研究內方差。μ為總體效應量，τ2為研究間異質性方差。當γ=0時退化為正態分布。
但實際中，由于每個研究的設計方案，研究對象和樣本容量等不同，導致每個研究的偏倚并不一定相同，于是我們對模型（1）進行改進。即假設納入Meta分析的各研究的效應量服從不同形狀參數的偏正態分布，亦即ξi~SN(0, 1, γi)，其中γi表示第i個研究的偏倚參數，這表明每個研究的偏倚均不相同。
在此條件下，觀測效應量yi的期望和方差分別為：
| $ \begin{array}{l} E\left( {{y_i}} \right) = \mu + \tau E\left( {{\xi _i}} \right) + E\left( {{\varepsilon _i}} \right) = \mu + \tau {\gamma _i}\sqrt {\frac{2}{\pi }} \sqrt {\frac{1}{{1 + \gamma _i^2}}} \\ Var\left( {{y_i}} \right) = {\tau ^2}Var\left( {{\xi _i}} \right) + Var\left( {{\varepsilon _i}} \right) = \delta _i^2 + {\tau ^2}\left( {1 - \frac{{2\gamma _i^2}}{{\pi \left( {1 + \gamma _i^2} \right)}}} \right) \end{array} $ |
記
| $ {m_i} = 1 - \frac{2}{\pi }\frac{{\gamma _i^2}}{{\gamma _i^2 + 1}} $ |
則Var(yi)=δi2+miτ2，從而mi也可等價地表示各研究的偏倚的參數。
1.1 觀測效應量yi的密度函數
| $ {\xi _i} = \frac{{{\gamma _i}}}{{1 + \gamma _i^2}}{t_i} + \frac{1}{{\sqrt {1 + \gamma _i^2} }}{x_i} $ |
其中ti~HN(0, 1)，xi~N(0, 1)。于是，觀測效應量yi可表示為：
| $ {y_i} = \mu + \frac{{\tau {\gamma _i}}}{{\sqrt {1 + \gamma _i^2} }}{t_i} + \frac{\tau }{{\sqrt {1 + \gamma _i^2} }}{x_i} + {\varepsilon _i} + {a_i}{t_i} + {s_i} $ |
其中，
ti和si相互獨立。
利用上式第二個等式可求得觀測效應量yi的密度函數為：
| $ \begin{array}{l} f\left( {{y_i}} \right) = \int\limits_0^\infty {\frac{2}{{\sqrt {2\pi } }}} {e^{\frac{{t_{^i}^2}}{2}}} \cdot \frac{1}{{\sqrt {2\pi \sigma _{^i}^2} }}{e^{\frac{{{{\left( {{y_i} - {a_i}{t_i} - \mu } \right)}^2}}}{{2\sigma _{^i}^2}}}}d{t_i}\\ \;\;\;\;\;\;\;\;\; = \frac{2}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}\phi \left( {\frac{{\left( {{y_i} - \mu } \right)}}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}} \right)\Phi \left( {\frac{{\tau {\gamma _i}}}{{\sqrt {{\tau ^2} + \delta _{^i}^2\left( {1 + \gamma _{^i}^2} \right)} }}\frac{{\left( {{y_i} - \mu } \right)}}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}} \right) \end{array} $ |
其中，。
1.2 總體效應量的估計
對改進的偏正態分布隨機效應模型采用馬爾科夫估計法，可得效應量μ的估計式為：
| $ \hat \mu = \frac{{\sum\limits_{i = 1}^n {\omega _i^*{y_i}} }}{{\sum\limits_{i = 1}^n {\omega _i^*} }} $ |
其中，，mi為（2）式，未知參數τ2和γi的估計在第三部分給出。
此外，總體效應量的方差為：
| $ {\rm{var}}\left( {\hat \mu } \right) = \frac{1}{{\sum\limits_{i = 1}^n {\omega _i^*} }} $ |
那么總體效應量的95%CI是
| $ \begin{array}{l} L{L_{\hat \mu }} = \hat \mu - 1.96\sqrt {{\rm{var}}\left( {\hat \mu } \right)} \\ U{L_{\hat \mu }} = \hat \mu + 1.96\sqrt {{\rm{var}}\left( {\hat \mu } \right)} \end{array} $ |
1.3 極大似然估計
為了求出總體效應量的估計值，需要先求參數τ2和γi的估計τ2和γi，i=1, …, n。根據模型（1），可采用極大似然估計法[9]對研究間異質性方差τ2和每個研究的偏倚γi進行估計。由觀測效應量yi的密度函數表達式（3）可得樣本的似然函數為：
| $ \begin{array}{l} L\left( {{\tau ^2},{\gamma _i},\mu } \right) = \prod\limits_{i = 1}^n {f\left( {{y_i}} \right)} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \prod\limits_{i = 1}^n {\frac{2}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}\phi \left( {\frac{{\left( {{y_i} - \mu } \right)}}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}} \right)\Phi \left( {\frac{{\tau {\gamma _i}\left( {{y_i} - \mu } \right)}}{{\sqrt {\left( {{\tau ^2} + \delta _{^i}^2} \right)\left( {{\tau ^2} + \delta _{^i}^2\left( {1 + \gamma _{^i}^2} \right)} \right)} }}} \right)} \end{array} $ |
從而對數似然函數為：
| $ \begin{array}{l} l\left( {{\tau ^2},{\gamma _i},\mu } \right) = {\left( {\frac{2}{\pi }} \right)^{\frac{n}{2}}} - \frac{1}{2}\sum\limits_{i = 1}^n {\ln \left( {{\tau ^2} + \delta _{^i}^2} \right)} - \frac{1}{2}\sum\limits_{i = 1}^n {\frac{{{{\left( {{y_i} - \mu } \right)}^2}}}{{{\tau ^2} + \delta _{^i}^2}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{i = 1}^n {\ln \left( {\Phi \left( {\frac{{\tau {\gamma _i}\left( {{y_i} - \mu } \right)}}{{\sqrt {\left( {{\tau ^2} + \delta _{^i}^2} \right)\left( {{\tau ^2} + \delta _{^i}^2\left( {1 + \gamma _{^i}^2} \right)} \right)} }}} \right)} \right)} \; \end{array} $ |
（5）式關于τ2求偏導，可得相應的似然方程為：
| $ \begin{array}{l} \frac{{\partial l{{\left( {{\tau ^2},{\gamma _i},\mu } \right)}^2}}}{{\partial {\tau ^2}}}\left| {{\tau ^2} - {{\hat \tau }^2}} \right. = - \frac{1}{2}\sum\limits_{i = 1}^n {\frac{1}{{{{\hat \tau }^2} + \delta _{^i}^2}}} + \frac{1}{2}\sum\limits_{i = 1}^n {{{\left( {{y_i} - \hat \mu } \right)}^2}} \frac{1}{{{{\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)}^2}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\sum\limits_{i = 1}^n {\left[ {\frac{{\phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}{{\Phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}} \right]} {\lambda _i} = 0 \end{array} $ |
其中，
| $ \begin{array}{l} {\lambda _i} = \frac{{{{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)\left[ {\frac{1}{\tau } - \hat \tau {{\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)}^{ - 1}} - \hat \tau {{\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)}^{ - 1}}} \right]}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }},\\ i = 1, \ldots ,n. \end{array} $ |
（5）式關于γi求偏導，可得相應的似然方程為：
| $ \begin{array}{l} \frac{{\partial l\left( {{\tau ^2},{\gamma _i},\mu } \right)}}{{\partial {\gamma _i}}}\left| {{\gamma ^2}{\rm{ = }}{{\hat \gamma }^2}} \right.{\rm{ = }}\frac{{\phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}{{\Phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}{\rho _i} = 0\\ i = 1, \ldots ,n \end{array} $ |
其中，。
上述似然方程（6）和（7）式無法求得估計值的解析表達式，可采用牛頓迭代法求解，即似然方程組的解就是研究間異質性方差和各個研究的偏倚的估計值，其中的初值可采用固定效應模型下的總體效應量代替，即用代替后，由（4）式計算。
2 仿真計算
采用蒙特卡羅仿真法，用matlab編程產生偏正態分布隨機數，進而產生四格表數據，根據四格表數據計算出效應量的觀測值和研究內方差，再利用公式（6）、（7）和（4）分別求出異質性方差、每個研究的偏倚和總體效應量的估計值。最后求出相應的OR值。仿真中總體效應量、異質性方差和每個研究的偏倚的真值分別取μ=0.5，τ2=1，γ=[0.28, 0.3, 0.27, 0.35, 0.28, 0.3, 0.32]選取研究個數n=7。具體步驟如下：
第一步，由偏正態分布SN(μ, τ2, γi)產生仿真計算所需要的7個原始數據θi，i=1, …, n用以產生四格表數據。
第二步，由偏正態分布SN(n, n/4, γi)分別產生試驗組的樣本容量niT和對照組的樣本容量niC。
第三步，由二項分布B(niT, piT)生成試驗組的病例數xiT，其中piT是均勻分布U(0.05, 0.65)上的隨機數。同理，由二項分布B(niC, piC)生成對照組的病例數xiC，其中piC由下式求得：
| $ p_i^C = p_i^T\exp \left( {{\theta _i}} \right)/\left[ {1 - p_i^T + p_i^T\exp \left( {{\theta _i}} \right)} \right] $ |
進而計算第i個研究的觀測效應量yi，它的值為四格表數據(xiT, niT-xiT, xiC, niC-xiC)中比值比的對數，即：
| $ {y_i} = \log \left\{ {\left[ {x_i^T/\left( {n_i^T - x_i^T} \right)} \right]/\left[ {x_i^C/\left( {n_i^C - x_i^C} \right)} \right]} \right\} $ |
同時，可計算效應值yi的研究內方差：
| $ \begin{array}{l} \hat \delta _i^2 = 1/\left( {x_i^T + 0.5} \right) + 1/\left( {\left( {n_i^T - x_i^T} \right) + 0.5} \right)\\ \;\;\;\;\;\; + 1/\left( {x_i^C + 0.5} \right) + 1/\left( {\left( {n_i^C - x_i^C} \right) + 0.5} \right) \end{array} $ |
其中每項的分母添加0.5是為了避免分母出現0時，造成計算溢出。
由第一步、第二步和第三步得到模擬數據見表 1，用R畫出效應值yi的核密度圖如圖（1）。
表1
偏正態分布模型下Meta分析的模擬數據
圖1
各個研究的效應量yi核密度圖
第四步，采用牛頓迭代法對公式（6）和（7）進行迭代求解得到估計值，，i=1, …, n。然后利用公式（4）得到總體效應量的估計值，見表 2。
表2
偏正態分布下和正態分布下各參數的真實值和估計值
第五步，利用公式OR=exp()求得表 2中的OR值。
從表 2可以看出：基于效應量yi服從形狀參數不同的偏正態分布的Meta分析隨機效應模型下，研究間的異質性方差比真值偏小，每個研究的偏倚的估計值接近真值，同時總體效應量的估計值和真實值相差不大，誤差為0.548 9。當采用效應量yi服從正態分布的Meta分析隨機效應模型時，總體效應量的估計值和真值相差0.566 3。由此可見，采用效應量服從形狀參數不同的偏正態分布的模型比效應量yi服從正態分布的模型，由于偏倚的剔除使得總體效應量的估計更接近真實值。另外，γi相同的偏正態分布隨機效應模型的異質性方差、每個研究的偏倚和總體效應量和γi不同的偏正態分布隨機效應模型幾近相同，但還是有一點細微的差別，這種較小的差別可能是由于偏倚本身對結果的影響較小所致，但結果的不同在一定程度上能夠說明該模型下各個研究偏倚的不同對總體效應量有一定的影響。最后，從ORi/OR0可以得出剔除基于形狀參數不同的偏正態分布的隨機Meta分析效應模型的偏倚所造成的異質性后的研究間異質性對總體效應量的影響比正態分布下的小。
3 實例分析
實例源于覃金春等的研究[10]，目的是通過對照試驗比較系統性腹膜后淋巴切除（SL）和非系統性腹膜后淋巴切除（USL）治療晚期上皮性卵巢癌的5年生存率。表 3為臨床試驗數據，以及每個研究的觀測效應量yi和研究內方差，圖（2）給出了觀測效應量yi的核密度圖，表 4為采用上述方法得到的計算結果。
表3
SL和USL治療晚期上皮性卵巢癌5年生存率的試驗數據
圖2
各個研究的效應量yi的核密度圖
表4
SL和USL治療晚期上皮性卵巢癌5年生存率的試驗結果
由圖（2）可以看出，觀測效應量yi近似服從偏正態分布。從表 4中可以看出，在偏正態隨機效應模型Meta分析模型下，OR值為0.531 3，比正態隨機效應模型Meta分析下的OR值0.474 8略大。OR值不同說明各個研究的偏倚對總體效應量有一定的影響。總體效應量的置信區間分別是（-0.953 7，-0.311 1）和（-1.079 8，-0.410 0），說明系統性腹膜后淋巴結切除（SL）對治療晚期上皮性卵巢癌有一定的效果，由表 4也可以看出在偏正態隨機效應Meta分析模型下，置信區間的長度是0.642 6，比正態模型下的置信區間長度0.669 8要短，由此表明偏正態隨機效應模型下效應量的區間估計的精度更高。
4 討論
目前，在進行Meta分析過程中，通常假設效應量服從正態分布，然后進行同質性檢驗，若研究間異質性較小，即P值小于0.05，則采用效應量服從正態分布的固定效應模型，否則，采用效應量服從正態分布的隨機效應模型，即只考慮到異質性，而未考慮由偏倚所引起的異質性。但實際上，由于對真實值進行測量時產生的片面或系統偏差，研究設計的缺陷和基于有缺陷資料的采集及其分析都會對結果產生一定的影響。因此，為了剔除偏倚對結果的影響，使結果更接近真實值本文提出了在效應量服從偏態分布的基礎上，每個研究的偏倚具體量化，即使各個研究的效應量服從不同形狀參數的偏正態分布，建立形狀參數不同的偏正態分布的隨機效應模型。在該模型下，通過計算估計每個研究的偏倚和去除所有偏倚之后異質性方差的大小。
本文通過蒙特卡羅仿真和實例分析對效應量yi服從形狀參數不同的偏正態分布隨機效應模型與效應量yi服從形狀參數相同的偏正態分布隨機效應模型和正態模型進行了比較，從理論上講，正態模型和形狀參數相同的偏正態分布隨機效應模型只是形狀參數不同的偏正態模型的特殊情況，即分別是形狀參數為零和形狀參數相同時的情形。從計算結果來看，在效應量yi服從形狀參數不同的偏正態分布的隨機效應模型時，采用本文提出的模型和方法得到的總體效應量的估計較另外兩種模型更接近真值，因為它剔除了各個研究的偏倚對總體效應量的影響。因此，無論在理論上還是實際應用中，本文提出的方法相對來說是可行和有效的。
Meta分析是匯總具有相同目的的多個獨立研究結果，剖析研究間的差異特征，對研究結果進行定量綜合評價的方法[1]。在傳統的Meta分析隨機效應模型中[2]，沒有考慮偏倚的影響，通常假設納入Meta分析的各個研究的效應量服從正態分布，由此得出的總體效應量只受研究間異質性和研究內方差的影響。實際上，異質性由多種因素造成，偏倚可認為是其中一個因素，然而由于偏倚本身的特殊性，如發表偏倚或研究者的偏好，使觀測效應量普遍偏高或偏低，即不是傳統模型中的正態分布，而是呈現非對稱的偏正態分布趨勢。因此，我們希望把偏倚從異質性中分離出來，單獨加以量化表征，使得模型能夠描述各研究觀測效應量非對稱分布的情形，而偏正態分布模型可以描述這類情形。于是，人們提出了對Meta分析中的偏倚進行量化[3-5]，即假設納入Meta分析中的各個研究的效應量服從相同形狀參數的偏正態分布，通過形狀參數對偏倚進行量化，進而可分析偏倚和異質性對結果的影響。但由于各個研究存在不同，如如何處理試驗過程中參與者的流失等，致使各個研究的偏倚不盡相同。因此本文提出將各個研究的偏倚通過不同的形狀參數進行量化，即假設納入Meta分析中的各個研究的效應量服從不同形狀參數的偏正態分布，進而建立基于形狀參數不同的偏正態分布下的隨機效應模型，然后采用極大似然估計法求出每個研究偏倚的估計值，進而可以分析各個研究的偏倚和研究間的異質性對總體效應量的影響。同時，本文還通過模擬計算和實例列舉，驗證本文提出的模型和算法的可行性和有效性。
1 基于偏正態的隨機效應模型
設效應量yi服從形狀參數相同的偏正態分布[6]隨機效應模型[5]為：
| $ {y_i}=\mu+\tau {\xi _i}+{\varepsilon _i}, i=1, \ldots, n $ |
其中yi，i=1, …, n是n個相互獨立的觀測效應量，n表示納入研究數。ξi~SN(0, 1, γ)，εi~N(0, δi2)且ξi和εi相互獨立，γ為形狀參數，表示每個研究的偏倚參數；δi2為研究內方差。μ為總體效應量，τ2為研究間異質性方差。當γ=0時退化為正態分布。
但實際中，由于每個研究的設計方案，研究對象和樣本容量等不同，導致每個研究的偏倚并不一定相同，于是我們對模型（1）進行改進。即假設納入Meta分析的各研究的效應量服從不同形狀參數的偏正態分布，亦即ξi~SN(0, 1, γi)，其中γi表示第i個研究的偏倚參數，這表明每個研究的偏倚均不相同。
在此條件下，觀測效應量yi的期望和方差分別為：
| $ \begin{array}{l} E\left( {{y_i}} \right) = \mu + \tau E\left( {{\xi _i}} \right) + E\left( {{\varepsilon _i}} \right) = \mu + \tau {\gamma _i}\sqrt {\frac{2}{\pi }} \sqrt {\frac{1}{{1 + \gamma _i^2}}} \\ Var\left( {{y_i}} \right) = {\tau ^2}Var\left( {{\xi _i}} \right) + Var\left( {{\varepsilon _i}} \right) = \delta _i^2 + {\tau ^2}\left( {1 - \frac{{2\gamma _i^2}}{{\pi \left( {1 + \gamma _i^2} \right)}}} \right) \end{array} $ |
記
| $ {m_i} = 1 - \frac{2}{\pi }\frac{{\gamma _i^2}}{{\gamma _i^2 + 1}} $ |
則Var(yi)=δi2+miτ2，從而mi也可等價地表示各研究的偏倚的參數。
1.1 觀測效應量yi的密度函數
| $ {\xi _i} = \frac{{{\gamma _i}}}{{1 + \gamma _i^2}}{t_i} + \frac{1}{{\sqrt {1 + \gamma _i^2} }}{x_i} $ |
其中ti~HN(0, 1)，xi~N(0, 1)。于是，觀測效應量yi可表示為：
| $ {y_i} = \mu + \frac{{\tau {\gamma _i}}}{{\sqrt {1 + \gamma _i^2} }}{t_i} + \frac{\tau }{{\sqrt {1 + \gamma _i^2} }}{x_i} + {\varepsilon _i} + {a_i}{t_i} + {s_i} $ |
其中，
ti和si相互獨立。
利用上式第二個等式可求得觀測效應量yi的密度函數為：
| $ \begin{array}{l} f\left( {{y_i}} \right) = \int\limits_0^\infty {\frac{2}{{\sqrt {2\pi } }}} {e^{\frac{{t_{^i}^2}}{2}}} \cdot \frac{1}{{\sqrt {2\pi \sigma _{^i}^2} }}{e^{\frac{{{{\left( {{y_i} - {a_i}{t_i} - \mu } \right)}^2}}}{{2\sigma _{^i}^2}}}}d{t_i}\\ \;\;\;\;\;\;\;\;\; = \frac{2}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}\phi \left( {\frac{{\left( {{y_i} - \mu } \right)}}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}} \right)\Phi \left( {\frac{{\tau {\gamma _i}}}{{\sqrt {{\tau ^2} + \delta _{^i}^2\left( {1 + \gamma _{^i}^2} \right)} }}\frac{{\left( {{y_i} - \mu } \right)}}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}} \right) \end{array} $ |
其中，。
1.2 總體效應量的估計
對改進的偏正態分布隨機效應模型采用馬爾科夫估計法，可得效應量μ的估計式為：
| $ \hat \mu = \frac{{\sum\limits_{i = 1}^n {\omega _i^*{y_i}} }}{{\sum\limits_{i = 1}^n {\omega _i^*} }} $ |
其中，，mi為（2）式，未知參數τ2和γi的估計在第三部分給出。
此外，總體效應量的方差為：
| $ {\rm{var}}\left( {\hat \mu } \right) = \frac{1}{{\sum\limits_{i = 1}^n {\omega _i^*} }} $ |
那么總體效應量的95%CI是
| $ \begin{array}{l} L{L_{\hat \mu }} = \hat \mu - 1.96\sqrt {{\rm{var}}\left( {\hat \mu } \right)} \\ U{L_{\hat \mu }} = \hat \mu + 1.96\sqrt {{\rm{var}}\left( {\hat \mu } \right)} \end{array} $ |
1.3 極大似然估計
為了求出總體效應量的估計值，需要先求參數τ2和γi的估計τ2和γi，i=1, …, n。根據模型（1），可采用極大似然估計法[9]對研究間異質性方差τ2和每個研究的偏倚γi進行估計。由觀測效應量yi的密度函數表達式（3）可得樣本的似然函數為：
| $ \begin{array}{l} L\left( {{\tau ^2},{\gamma _i},\mu } \right) = \prod\limits_{i = 1}^n {f\left( {{y_i}} \right)} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \prod\limits_{i = 1}^n {\frac{2}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}\phi \left( {\frac{{\left( {{y_i} - \mu } \right)}}{{\sqrt {{\tau ^2} + \delta _{^i}^2} }}} \right)\Phi \left( {\frac{{\tau {\gamma _i}\left( {{y_i} - \mu } \right)}}{{\sqrt {\left( {{\tau ^2} + \delta _{^i}^2} \right)\left( {{\tau ^2} + \delta _{^i}^2\left( {1 + \gamma _{^i}^2} \right)} \right)} }}} \right)} \end{array} $ |
從而對數似然函數為：
| $ \begin{array}{l} l\left( {{\tau ^2},{\gamma _i},\mu } \right) = {\left( {\frac{2}{\pi }} \right)^{\frac{n}{2}}} - \frac{1}{2}\sum\limits_{i = 1}^n {\ln \left( {{\tau ^2} + \delta _{^i}^2} \right)} - \frac{1}{2}\sum\limits_{i = 1}^n {\frac{{{{\left( {{y_i} - \mu } \right)}^2}}}{{{\tau ^2} + \delta _{^i}^2}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{i = 1}^n {\ln \left( {\Phi \left( {\frac{{\tau {\gamma _i}\left( {{y_i} - \mu } \right)}}{{\sqrt {\left( {{\tau ^2} + \delta _{^i}^2} \right)\left( {{\tau ^2} + \delta _{^i}^2\left( {1 + \gamma _{^i}^2} \right)} \right)} }}} \right)} \right)} \; \end{array} $ |
（5）式關于τ2求偏導，可得相應的似然方程為：
| $ \begin{array}{l} \frac{{\partial l{{\left( {{\tau ^2},{\gamma _i},\mu } \right)}^2}}}{{\partial {\tau ^2}}}\left| {{\tau ^2} - {{\hat \tau }^2}} \right. = - \frac{1}{2}\sum\limits_{i = 1}^n {\frac{1}{{{{\hat \tau }^2} + \delta _{^i}^2}}} + \frac{1}{2}\sum\limits_{i = 1}^n {{{\left( {{y_i} - \hat \mu } \right)}^2}} \frac{1}{{{{\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)}^2}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\sum\limits_{i = 1}^n {\left[ {\frac{{\phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}{{\Phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}} \right]} {\lambda _i} = 0 \end{array} $ |
其中，
| $ \begin{array}{l} {\lambda _i} = \frac{{{{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)\left[ {\frac{1}{\tau } - \hat \tau {{\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)}^{ - 1}} - \hat \tau {{\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)}^{ - 1}}} \right]}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }},\\ i = 1, \ldots ,n. \end{array} $ |
（5）式關于γi求偏導，可得相應的似然方程為：
| $ \begin{array}{l} \frac{{\partial l\left( {{\tau ^2},{\gamma _i},\mu } \right)}}{{\partial {\gamma _i}}}\left| {{\gamma ^2}{\rm{ = }}{{\hat \gamma }^2}} \right.{\rm{ = }}\frac{{\phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}{{\Phi \left( {\frac{{\hat \tau {{\hat \gamma }_i}\left( {{y_i} - \hat \mu } \right)}}{{\sqrt {\left( {{{\hat \tau }^2} + \delta _{^i}^2} \right)\left( {{{\hat \tau }^2} + \delta _{^i}^2\left( {1 + \hat \gamma _i^2} \right)} \right)} }}} \right)}}{\rho _i} = 0\\ i = 1, \ldots ,n \end{array} $ |
其中，。
上述似然方程（6）和（7）式無法求得估計值的解析表達式，可采用牛頓迭代法求解，即似然方程組的解就是研究間異質性方差和各個研究的偏倚的估計值，其中的初值可采用固定效應模型下的總體效應量代替，即用代替后，由（4）式計算。
2 仿真計算
采用蒙特卡羅仿真法，用matlab編程產生偏正態分布隨機數，進而產生四格表數據，根據四格表數據計算出效應量的觀測值和研究內方差，再利用公式（6）、（7）和（4）分別求出異質性方差、每個研究的偏倚和總體效應量的估計值。最后求出相應的OR值。仿真中總體效應量、異質性方差和每個研究的偏倚的真值分別取μ=0.5，τ2=1，γ=[0.28, 0.3, 0.27, 0.35, 0.28, 0.3, 0.32]選取研究個數n=7。具體步驟如下：
第一步，由偏正態分布SN(μ, τ2, γi)產生仿真計算所需要的7個原始數據θi，i=1, …, n用以產生四格表數據。
第二步，由偏正態分布SN(n, n/4, γi)分別產生試驗組的樣本容量niT和對照組的樣本容量niC。
第三步，由二項分布B(niT, piT)生成試驗組的病例數xiT，其中piT是均勻分布U(0.05, 0.65)上的隨機數。同理，由二項分布B(niC, piC)生成對照組的病例數xiC，其中piC由下式求得：
| $ p_i^C = p_i^T\exp \left( {{\theta _i}} \right)/\left[ {1 - p_i^T + p_i^T\exp \left( {{\theta _i}} \right)} \right] $ |
進而計算第i個研究的觀測效應量yi，它的值為四格表數據(xiT, niT-xiT, xiC, niC-xiC)中比值比的對數，即：
| $ {y_i} = \log \left\{ {\left[ {x_i^T/\left( {n_i^T - x_i^T} \right)} \right]/\left[ {x_i^C/\left( {n_i^C - x_i^C} \right)} \right]} \right\} $ |
同時，可計算效應值yi的研究內方差：
| $ \begin{array}{l} \hat \delta _i^2 = 1/\left( {x_i^T + 0.5} \right) + 1/\left( {\left( {n_i^T - x_i^T} \right) + 0.5} \right)\\ \;\;\;\;\;\; + 1/\left( {x_i^C + 0.5} \right) + 1/\left( {\left( {n_i^C - x_i^C} \right) + 0.5} \right) \end{array} $ |
其中每項的分母添加0.5是為了避免分母出現0時，造成計算溢出。
由第一步、第二步和第三步得到模擬數據見表 1，用R畫出效應值yi的核密度圖如圖（1）。
表1
偏正態分布模型下Meta分析的模擬數據
圖1
各個研究的效應量yi核密度圖
第四步，采用牛頓迭代法對公式（6）和（7）進行迭代求解得到估計值，，i=1, …, n。然后利用公式（4）得到總體效應量的估計值，見表 2。
表2
偏正態分布下和正態分布下各參數的真實值和估計值
第五步，利用公式OR=exp()求得表 2中的OR值。
從表 2可以看出：基于效應量yi服從形狀參數不同的偏正態分布的Meta分析隨機效應模型下，研究間的異質性方差比真值偏小，每個研究的偏倚的估計值接近真值，同時總體效應量的估計值和真實值相差不大，誤差為0.548 9。當采用效應量yi服從正態分布的Meta分析隨機效應模型時，總體效應量的估計值和真值相差0.566 3。由此可見，采用效應量服從形狀參數不同的偏正態分布的模型比效應量yi服從正態分布的模型，由于偏倚的剔除使得總體效應量的估計更接近真實值。另外，γi相同的偏正態分布隨機效應模型的異質性方差、每個研究的偏倚和總體效應量和γi不同的偏正態分布隨機效應模型幾近相同，但還是有一點細微的差別，這種較小的差別可能是由于偏倚本身對結果的影響較小所致，但結果的不同在一定程度上能夠說明該模型下各個研究偏倚的不同對總體效應量有一定的影響。最后，從ORi/OR0可以得出剔除基于形狀參數不同的偏正態分布的隨機Meta分析效應模型的偏倚所造成的異質性后的研究間異質性對總體效應量的影響比正態分布下的小。
3 實例分析
實例源于覃金春等的研究[10]，目的是通過對照試驗比較系統性腹膜后淋巴切除（SL）和非系統性腹膜后淋巴切除（USL）治療晚期上皮性卵巢癌的5年生存率。表 3為臨床試驗數據，以及每個研究的觀測效應量yi和研究內方差，圖（2）給出了觀測效應量yi的核密度圖，表 4為采用上述方法得到的計算結果。
表3
SL和USL治療晚期上皮性卵巢癌5年生存率的試驗數據
圖2
各個研究的效應量yi的核密度圖
表4
SL和USL治療晚期上皮性卵巢癌5年生存率的試驗結果
由圖（2）可以看出，觀測效應量yi近似服從偏正態分布。從表 4中可以看出，在偏正態隨機效應模型Meta分析模型下，OR值為0.531 3，比正態隨機效應模型Meta分析下的OR值0.474 8略大。OR值不同說明各個研究的偏倚對總體效應量有一定的影響。總體效應量的置信區間分別是（-0.953 7，-0.311 1）和（-1.079 8，-0.410 0），說明系統性腹膜后淋巴結切除（SL）對治療晚期上皮性卵巢癌有一定的效果，由表 4也可以看出在偏正態隨機效應Meta分析模型下，置信區間的長度是0.642 6，比正態模型下的置信區間長度0.669 8要短，由此表明偏正態隨機效應模型下效應量的區間估計的精度更高。
4 討論
目前，在進行Meta分析過程中，通常假設效應量服從正態分布，然后進行同質性檢驗，若研究間異質性較小，即P值小于0.05，則采用效應量服從正態分布的固定效應模型，否則，采用效應量服從正態分布的隨機效應模型，即只考慮到異質性，而未考慮由偏倚所引起的異質性。但實際上，由于對真實值進行測量時產生的片面或系統偏差，研究設計的缺陷和基于有缺陷資料的采集及其分析都會對結果產生一定的影響。因此，為了剔除偏倚對結果的影響，使結果更接近真實值本文提出了在效應量服從偏態分布的基礎上，每個研究的偏倚具體量化，即使各個研究的效應量服從不同形狀參數的偏正態分布，建立形狀參數不同的偏正態分布的隨機效應模型。在該模型下，通過計算估計每個研究的偏倚和去除所有偏倚之后異質性方差的大小。
本文通過蒙特卡羅仿真和實例分析對效應量yi服從形狀參數不同的偏正態分布隨機效應模型與效應量yi服從形狀參數相同的偏正態分布隨機效應模型和正態模型進行了比較，從理論上講，正態模型和形狀參數相同的偏正態分布隨機效應模型只是形狀參數不同的偏正態模型的特殊情況，即分別是形狀參數為零和形狀參數相同時的情形。從計算結果來看，在效應量yi服從形狀參數不同的偏正態分布的隨機效應模型時，采用本文提出的模型和方法得到的總體效應量的估計較另外兩種模型更接近真值，因為它剔除了各個研究的偏倚對總體效應量的影響。因此，無論在理論上還是實際應用中，本文提出的方法相對來說是可行和有效的。
