單臂試驗連續型數據的貝葉斯Meta分析方法及實現

張天嵩

系統評價和Meta分析是循證醫學重要的證據信息來源[1，2]，通過Meta分析對多個研究結果的統計綜合被廣泛用于評價不同干預的相對效應[3]。在牛津循證醫學中心2011年發布的證據水平分級標準中，將以RCT為基礎的系統評價作為評價某種干預措施受益或危害的1級水平[4]。雖然RCT是公認的最佳治療性研究設計方案，但來自RCT的證據常缺乏或不足，難以回答所有臨床問題，或有時RCT不需要、不可行或不恰當，如新藥1期臨床試驗，則可以采用單臂試驗設計方法[5，6]，國、內外均有依據單臂試驗批準新藥上市的案例[7]。所謂單臂試驗是指將具有某種靶健康狀況的個體樣本給予試驗性治療并觀察某種治療效應的設計[6]，所有受試對象都在同一觀察組，不設立其他實驗組或對照組[5]。近年來，國內外循證醫學工作者開始探索單臂試驗數據的Meta分析方法。單臂試驗的觀察指標按資料性質可分為定量與定性資料。對于定性資料，多以比例、率、比數和發病密度等為效應量[5]，目前有眾多研究文獻介紹此類數據在貝葉斯和頻率學框架下的Meta分析方法，以及基于R、Stata、SAS、RevMan等軟件實現[5, 8-13]；而對連續型數據的Meta分析方法則研究或介紹不多，有少數介紹使用Stata軟件[5]和Meta-Analyst軟件[14]基于頻率學框架下采用倒方差法進行Meta分析的文獻，難以滿足實踐需要。本文介紹基于經典的正態-正態層次模型(NNHM)框架下單臂試驗連續型數據的貝葉斯Meta分析方法，并通過實例來介紹R軟件實現過程。

1 方法

1.1 研究數據 數據來源于2個Meta分析。數據1：來自Dai等的單臂試驗的Meta分析文獻[15]，主要觀察來氟米特治療銀屑病性關節炎的有效性和安全性，本文選取的觀察指標為銀屑病皮損面積和嚴重度指數(PASI)評分，為治療后評分與基線評分的變化值，具體數據信息為每個試驗的樣本量、評分均數及標準差，在表1中分別以“n”、“mean”和“sd”表示；數據2：來自Nagpal等的Meta分析文獻[16]，主要觀察造血干細胞療法對卒中早期臨床干預的有效性和安全性，該文獻納入了單臂試驗和對照試驗的數據，本文選取單臂試驗的相關數據，觀察指標為美國國立衛生研究院卒中量表(NIHSS)評分，為治療后評分與基線評分的變化值，具體數據信息為評分的均數及其95%CI下限和上限，在表1中分別以“mean”、“ll”、“ul”表示。數據見表1。

表1 分析數據

1.2 模型與方法 經典的Meta分析隨機效應模型基于NNHM框架[17]：假設納入Meta分析有i=1,…，N個獨立研究，每個研究的真實值為θi，其估計值及估計值標準誤分別為yi和Si，研究間真實效應變異為τ2，總的合并效應為μ，則模型有兩層。

yi～N(θi，Si2)為抽樣模型，假定yi服從未知均值θi和已知標準誤Si的正態分布；θi～ N(μ,τ2)為參數模型，假定θi服從均數為μ、方差為τ2(標準差為τ)的正態分布，其中未知參數μ為主要感興趣的參數，而異質性參數τ則一般被認為是討厭參數。在頻率學框架下，該模型可以使用Stata、R、SAS、RevMan等軟件來擬合和實現，常用的參數估計方法有最大似然法(ML)和約束最大似然法(REML)等。在貝葉斯框架下，通常需要根據實際情況為參數μ和τ指定合理的先驗。一般情況下，對于參數μ，指定為無信息先驗分布；而對τ，可以指定為無信息先驗分布如均勻分布，或弱信息先驗分布如半正態分布等，該模型可由貝葉斯軟件WinBugs等軟件擬合。

1.3 模型擬合 在貝葉斯框架下，NNHM可由R軟件(3.4.4版)的bayesmeta包的bayesmeta()函數來擬合[18]。在本研究中，參數μ設定為無信息先驗μ～N(0，104)；參數τ指定為弱信息先驗，分別為半正態(HN)先驗和半柯西(HC)先驗，尺度均設為1。以forestplot()函數繪制森林圖。同時，為了比較與頻率學框架下Meta分析結果，以metafor包的rma.uni()函數擬合該模型，選用REML法估計參數τ。具體代碼及簡要解釋見表2。

表2 單臂試驗連續型數據貝葉斯Meta分析代碼(R軟件)及簡要解釋

續表2 tau.prior=function(x){(dhalfcauchy(x, scale=1))})print(bma22)##繪制森林圖forestplot(bma22)##擬合模型3:頻率學框架,隨機效應模型,采用REML法估計參數τ。meta23<-rma.uni(yi=score,sei=se,method="REML",data=stemcell)summary(meta23)

2 結果

2個研究數據貝葉斯和頻率學Meta分析方法對效應參數μ和研究間異質性τ的估計結果如表3所示，2種方法獲得參數μ點估計值有差異，95%CI寬度有明顯的不同；參數τ估計值差異較大，基于HN(1)先驗的貝葉斯方法似乎可以獲得更為精確的估計值。需要指出的是，bayesmeta包提供了更多的貝葉斯分析結果信息，如在其“marginal posterior summary”結果顯示部分，提供2個主要參數μ和τ后驗眾數、后驗均數、后驗中位數，及相應95%CI的下限和上限，為了與森林圖結果一致，表1中所示選取的是后驗中位數及95% CI的下限和上限結果。

數據1的森林圖如圖1和2所示。森林圖中上半部分顯示了納入Meta分析的每個研究的效應量(即模型中的yi)的點估計及基于相應標準誤獲得的95%CI，以黑色表示，即森林圖中的“quoted estimate”；以及每個研究真實效應參數(即模型中的θ_i)的點估計及95%CI,以灰色表示，即森林圖中的“shrinkage estimate”；森林圖底部顯示合并效應量(平均值)點估計(后驗中位數)和95%CI(與圖中“mean”相應)，以及預測值及區間(與圖中“prediction”相應)。

表3 基于NNHM框架下貝葉斯和頻率學Meta分析結果

圖1 數據1森林圖[τ為HN(1)先驗]圖2 數據1森林圖[τ為HC(1)先驗]

3 討論

單臂試驗在醫學研究中較為常見，如果其測量結局為連續型數據，則具有一個共同特征[19]：針對某一測量結局指標，除了干預實施后進行測量外，可能會對每一個參與者在干預實施前進行基線測量，因此可以獲得最終測量值與基線測量值的差值，稱為基線改變值或改變評分，所以系統評價者在Meta分析提取數據時可能會考慮同時提取基線改變值和最終測量值。對于連續型數據，原始研究對測量結果可能存在多種報告形式，常見的有：①測量指標均數及相應標準誤(或方差)；②受試者樣本量、測量指標均數及相應標準差；③測量指標的均數及95% CI；④測量指標的中位數及上下4分位數間距等；⑤測量指標的中位數、最大值及最小值等。無論何種形式的數據，只要最終能獲得每個研究的效應量及其方差或標準誤，即可進行Meta分析，一般可以通過2種途徑：一是向原始文獻的通訊作者發郵件索要原始數據；二是通過數據轉換，具體方法可以參看國內外相關文獻[19-23]。本文數據1已知樣本量、測量指標均數及相應標準差，可以通過公式“均數標準誤=標準差/樣本量的平方根”獲得均數標準誤[24]；數據2(因文獻[16]未報告每個研究的樣本量，本研究聚焦在貝葉斯Meta分析方法介紹，簡單處理為假定每個研究樣本量很大，95%CI來自正態分布)已知測量指標95%CI，則可以通過公式“均數標準誤=(95%CI上限-95%CI下限)/(2×1.96)”來獲得均數標準誤[19，22]，如果研究樣本量小，或者95%CI由t分布計算所得，則上述公式中的1.96則需要一個更大一點t分布值代替[19，22]，該值可以用自由度(研究樣本量-1)查t分布表或用Stata、R、Microsof Excel等軟件計算獲得。

一旦獲得每個研究觀察指標測量值及其相應標準誤，則可以在NNHM框架下采用隨機效應模型來進行Meta分析，主要有兩個參數真實效應(real-valued effect)μ和異質性(heterogeneity)τ，如果異質性為0，則退變為固定效應模型。在NNHM框架下的分析方法有頻率學和貝葉斯兩大策略，目前更為常用的頻率學策略一般分2步，第1步估計異質性，第2步基于異質性估計派生出效應估計，但目前有很多估計異質性算法，共同的問題是對于異質性方差分量估計易趨向于0，導致選擇固定效應模型。貝葉斯策略近年來也已逐步應用于Meta分析中，綜合了未知參數的部體信息、樣本信息及先驗信息，特別是可以允許(實際上也需要)指定先驗信息用于數據計算。在貝葉斯隨機效應模型中，根據需要很容易對效應參數μ指定為無信息先驗或弱信息先驗；對研究間異質性參數τ，則常需要根據效應量、研究數量等因素綜合考慮指定為信息先驗，特別是僅有幾個研究納入Meta分析時。有學者認為，除非在研究數量少或期望信息先驗等情況下，一般均宜將τ指定為均勻分布；研究數量少時，推薦τ指定為重尾分布，如半學生t(half-Student-t)分布、HC分布、HN分布等弱信息先驗[25,26]。本研究中所用數據納入Meta分析的研究數量均較少，因此將τ的先驗指定為HN(1),另指定HC(1)作為敏感性分析，以比較不同先驗對參數估計的影響。本研究發現，采用貝葉斯分析方法，納入Meta分析的研究數量較少(如3個研究時)，給定HN(1)和HC(1)先驗，未知參數μ和τ的后驗分布差異較大；而研究數量為中等(如7個時)，兩個未知參數后驗分布則較為接近，提示，Meta分析納入研究數量少時，對于研究間異質性估計要謹慎，值得進一步研究和探索。

貝葉斯方法通常需要高維積分方法等，在計算方面較難實現，但隨著現代計算機技術的出現和發展，特別是馬爾可夫鏈蒙特卡洛(MCMC)方法的建立，使得計算方法制約問題得以解決，但使用MCMC方法需要一定量的專業知識如收斂性診斷和實踐經驗等。因在NNHM框架貝葉斯隨機效應模型中，只有μ和τ兩個未知參數需要估計和推斷，可以利用數值積分方法或重要性抽樣使計算簡單化[27]，研發bayesmeta包的目的就是為了簡單易行地實現基于NNHM框架下完全貝葉斯分析策略，采用稱為散度限制性條件鑲嵌(DIRECT)策略[28]的新算法獲得未知參數的后驗分布，使用者無需擔心MCMC方法相關設置、診斷和后加工處理等貝葉斯方法專業問題。使用bayesmeta包可以快速獲得精確的后驗概括性結果，并能提供類似于經典頻率學Meta分析森林圖的圖示化結果；后驗分布可以準解析形式或預測估計和收縮估計(shrinkage estimation)等高級方法如存取。bayesmeta包還提供了其他結果信息，如預測未來一個新的研究真實效應參數θk+1值，即森林圖中“prediction”相應部分，也可以通過summary()函數獲得具體數字化結果。

最后，因為Bayesmeta包具有以下優點：①計算速度快，精確度和可重復性強；②操作簡單，不需要深度掌握MCMC算法等專業知識；③允許快速敏感性檢驗、便捷地大規模模擬等，所以值得推廣和使用。該軟件包除了用于單臂試驗連續型數據Meta分析外，還可以用于RCT等雙臂設計的二分類數據、連續型數據、生存數據等貝葉斯Meta分析。