單臂試驗(yàn)連續(xù)型數(shù)據(jù)的貝葉斯Meta分析方法及實(shí)現(xiàn)

張?zhí)灬?/p>

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

1 方法

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

表1 分析數(shù)據(jù)

1.2 模型與方法 經(jīng)典的Meta分析隨機(jī)效應(yīng)模型基于NNHM框架[17]：假設(shè)納入Meta分析有i=1,…，N個(gè)獨(dú)立研究，每個(gè)研究的真實(shí)值為θi，其估計(jì)值及估計(jì)值標(biāo)準(zhǔn)誤分別為yi和Si，研究間真實(shí)效應(yīng)變異為τ2，總的合并效應(yīng)為μ，則模型有兩層。

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

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

表2 單臂試驗(yàn)連續(xù)型數(shù)據(jù)貝葉斯Meta分析代碼(R軟件)及簡(jiǎn)要解釋

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

2 結(jié)果

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

數(shù)據(jù)1的森林圖如圖1和2所示。森林圖中上半部分顯示了納入Meta分析的每個(gè)研究的效應(yīng)量(即模型中的yi)的點(diǎn)估計(jì)及基于相應(yīng)標(biāo)準(zhǔn)誤獲得的95%CI，以黑色表示，即森林圖中的“quoted estimate”；以及每個(gè)研究真實(shí)效應(yīng)參數(shù)(即模型中的θ_i)的點(diǎn)估計(jì)及95%CI,以灰色表示，即森林圖中的“shrinkage estimate”；森林圖底部顯示合并效應(yīng)量(平均值)點(diǎn)估計(jì)(后驗(yàn)中位數(shù))和95%CI(與圖中“mean”相應(yīng))，以及預(yù)測(cè)值及區(qū)間(與圖中“prediction”相應(yīng))。

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

圖1 數(shù)據(jù)1森林圖[τ為HN(1)先驗(yàn)]圖2 數(shù)據(jù)1森林圖[τ為HC(1)先驗(yàn)]

3 討論

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

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

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

最后，因?yàn)锽ayesmeta包具有以下優(yōu)點(diǎn)：①計(jì)算速度快，精確度和可重復(fù)性強(qiáng)；②操作簡(jiǎn)單，不需要深度掌握MCMC算法等專業(yè)知識(shí)；③允許快速敏感性檢驗(yàn)、便捷地大規(guī)模模擬等，所以值得推廣和使用。該軟件包除了用于單臂試驗(yàn)連續(xù)型數(shù)據(jù)Meta分析外，還可以用于RCT等雙臂設(shè)計(jì)的二分類數(shù)據(jù)、連續(xù)型數(shù)據(jù)、生存數(shù)據(jù)等貝葉斯Meta分析。