網狀Meta分析（NMA）是一種通過綜合多個臨床研究數據且比較多種干預措施的療效和安全性的統計技術，可為證據網絡中所有干預方案提供優劣排序結果，并為臨床決策提供直接證據支撐。當前，NMA通常基于同一類型數據集的匯總，而對于實現跨研究設計與跨數據格式的數據集合并仍存在方法學與軟件操作上難點。R軟件crossnma程序包基于貝葉斯框架與馬爾可夫鏈蒙特卡羅算法，將三級分層模型擴展到標準NMA數據模型中以實現不同數據類型的差異化合并。crossnma程序包通過引入模型變量充分考慮不同類型數據間合并所帶來的偏倚風險對結果的影響。此外，該程序包還提供結果輸出和簡易圖形繪制等功能，這為實現跨研究設計與跨數據格式證據的NMA提供可能。本研究將通過4個個體參與者數據集與2個聚合數據集的實例對crossnma程序包模型方法和軟件操作進行實例演示和講解。
引用本文： 劉潤本, 張智焮, 黃承楊, 李昊陽, 趙軼凡, 張超, 羅杰. 應用R軟件crossnma程序包實現跨研究設計與跨數據格式的網狀Meta分析. 中國循證醫學雜志, 2023, 23(10): 1197-1203. doi: 10.7507/1672-2531.202303092 復制
網狀Meta分析（network meta-analysis，NMA）是一種通過綜合多個臨床研究數據且比較多種干預措施的療效和安全性的統計技術，其可為證據網絡中所有干預方案提供優劣排序結果，并可為臨床決策提供直接證據支撐。通常，NMA在數據類型的選取時，僅局限于隨機對照試驗（randomized controlled trials，RCT），是為了降低不確定性混雜因素介入，以保證網狀數據結果準確性和穩定性。隨著統計學家在方法學上不斷更新，現已可實現將不同研究設計和數據格式的數據集融入到同一個NMA的框架中，以盡可能達到更大樣本結果的統計學效能。
2022年，Tasnim教授團隊通過網狀Meta回歸（network meta-regression，NMR）模型來實現跨研究設計和跨數據格式的數據集（cross-design evidence and cross-format data）合并，并基于R語言平臺開發了crossnma程序包以便研究者使用。NMR模型不僅可將RCT和非隨機研究（nonrandomized studies，NRS）進行融合，還可將聚合數據（aggregate data，AD）和個體參與者數據（individual participant data，IPD）類型數據融入同一NMA當中。此外，NMR模型還考慮到跨研究設計和跨數據格式的數據集間混雜因素和協變量所帶來的潛在異質性等問題。例如：對在已發表的文獻中提取出的AD數據而言，其協變量需要基于整個研究層面水平進行，而IPD數據，則需要從單個樣本病例水平來予以考量和設計，以避免兩個類型數據偏倚和混雜因素相互干擾。
NMA中的證據通常來自RCT，盡管人們認為RCT與NRS相比具有較低的偏倚風險，但Cochrane的一項系統評價發現RCT與NRS在提供不同干預效果估計上并無實質性差異[1]，此外，實證研究證實RCT中可能存在不同類型的偏倚。Schulz等[2]發現分配隱藏不足或缺乏盲法的RCT更傾向于夸大估計的治療效果，并會提供有偏的主觀結果。Chalmers等[3]在非盲或隨機化不規范的研究中發現干預組與對照組間的療效存在較大差異。與此同時，在分配隱藏較差或缺乏盲法的研究中，Wood等[4]發現相對治療效果往往被評估者主觀判斷或患者自我報告夸大。這些方面結論為RCT與NRS跨研究設計數據類型的合并提供了方法學基礎。
前期研究通常采用匹配調整間接比較[5]和模擬治療比較方法[6]通過二次分配權重法和回歸模型法來調整來自于IPD與AD數據集的效應修正因子。然而，這種效應修正因子的調整方法需要對每個干預措施來優化校正，因此，需要每個干預措施比較中至少有一個IPD數據。Jansen等[7]將三級分層模型擴展到標準NMA數據模型當中，該方法通過整合AD數據集潛在的IPD分布來實現將二分類結局數據擴展到其他類型數據，從而達到了IPD與AD差異化數據合并[8]。
考慮到不同偏倚風險（risk of bias，RoB）差異，crossnma程序包采用了Dias等[9]開發的模型將每個研究的RoB因子考慮到模型當中，以降低跨研究設計間數據合并可能帶來的偏差。本研究將通過4個IPD數據集[10-13]與2個AD數據集[14-15]的實例對crossnma程序包模型方法和軟件操作進行實例演示和講解。
1 R軟件crossnma程序包方法學簡介
R軟件crossnma程序包根據IPD網狀數據與AD網狀數據特征，運用網狀Meta回歸模型分別擬定對應的回歸模型，并將兩個數據模型結果進行整合。
1.1 針對IPD數據的未調整網狀Meta回歸模型
 |
公式（1）中，
服從伯努利分布；
為在第j個IPD數據中，參考干預b的效應量；
為獨立效應，
為協變量，
與
可在兩式中相抵消，此時，式中剩下相對效應
、組內協變量校正因子
（組內協變量系數與干預j的平均協變量的乘積）和組間與組內協變量校正因子差值
。
1.2 針對AD數據的未調整網狀Meta回歸模型
 |
公式（2）中，
為研究中每臂的事件數服從二項式分布，此時模型較為簡單，對照組的概率僅為隨機擾動項
，它可在二者中相抵消，則關于AD數據相對效應的Meta回歸中自變量只有平均協變量
。
此時
作為似然函數是隨機擾動項，服從
，
則先驗服從
，對應的貝葉斯公式為公式（3）。
 |
那么先驗
從何得來，本研究僅使用NRS研究得到均值為
，方差為
的后驗分布，再將這個后驗分布作為RCT的先驗，將均值加上
，方差除以一個0到1的膨脹因子
，經過如此調整就可以得到先驗
。本研究并沒有使用來自NRS研究的AD，而是使用了RCT研究的AD，在程序包中可設置參數將RCT研究的AD視為NRS數據來參與計算。
1.3 針對IPD數據的風險調整網狀Meta回歸模型
在偏差調整模型中作者使用了兩種不同的方法。偏差調整模型I［公式（4）］是將偏差使用加法效應和乘法效應合并到每一個研究，并用二元參數
（0為沒有偏差，1為存在偏差）來衡量是否存在偏差，參數
服從零一分布即
，其中
，
與
控制著分布的偏斜程度。當
，說明分布沒有偏斜其RoB為未知；當
，此時分布有著高偏倚或低偏倚，在本程序包中面對不同偏倚風險的證據均可設置不同的參數來相匹配。
為乘法偏差效應系數，
為加法偏差效應系數。
 |
而偏差調整模型Ⅱ［公式（5）］則是使用一個服從雙峰正態分布的參數
來代替
，
服從隨機效應的混合正態分布
，并將偏倚效應納入其中，此時偏差調整變為兩個分布之和。
 |
1.4 針對AD數據的風險調整網狀Meta回歸模型
 |
相較于IPD數據的風險調整模型，AD數據的風險調整模型［公式（6）］略顯簡單。在估算
方面還有另一種計算方法，使用對數變換后的
再加上隨機擾動項
，即
。對于偏差效應
有著兩種假設，當在隨機效應中
，而在固定效應中
。
為平均偏差效應，當b為非活性藥物干預措施時
，當b和k為活性藥物干預措施時
；
取0時則為零均值偏差效應，取
為共同偏差效應；其中
為偏差方向，當
取值為0時則結果偏向于b，當
取值為1時結果偏向于k。
可由參數bias.group控制，當其取0時為無偏差調整，取1時
為非活性藥物干預措施，取2時
為活性藥物干預措施。
 |
還有另一種加法偏差效應的運用［公式（7）］，此時原來的OR指數
已變為加入偏差效應的
即
，此時
。偏差調整模型Ⅱ［公式（8）］更是大同小異，由此假設一個共同模型來總結相對效應
。
 |
2 R軟件crossnma程序包操作簡介
2.1 程序包下載與安裝
本研究使用的R 4.2.2，具體安裝方式可參閱《R軟件Metafor程序包在Meta分析中的應用》[16]。值得注意的是，crossnma程序包的運算是基于AJGS軟件來進行數據處理的，因此，需要事先安裝JAGS軟件。目前，JAGS軟件可免費從官方網站下載，具體下載方式請參考《R軟件調用JAGS軟件實現網狀Meta分析》[17]一文。
在使用R軟件crossnma程序包來執行IPD和AD數據NMA時，需要安裝和加載其他輔助性程序包，如：netmeta與meta。具體安裝和加載代碼，見表1。
表1
程序安裝與加載代碼表
2.2 數據介紹與預處理
2.2.1 數據介紹
復發緩解型多發性硬化癥（relapsing remitting multiple sclerosis，RRMS）為多發性硬化癥（multiple sclerosis，MS）患病人數最多的一種臨床亞型，其臨床表現為神經系統退行性慢性炎癥，由于其預后較差、人數較多故而倍受關注。該程序包研發者將以不同的干預措施治療復發緩解型多發性硬化癥后2年內是否復發的IPD和AD兩種不同數據格式臨床數據集為實例。其中，IPD數據集為個體受試者數據，含1 945例數據，4種干預措施（A、B、C和D），而AD數據集僅有4例數據，兩種干預措施（A和C）。本研究IPD數據集其中包含一個NRS研究（Swiss Multiple Sclerosis Cohort[10]）和三個RCT研究（分別來自CONFIRM[11]、DEFINE[12]和AFFIRM[13]），詳見表2。AD數據集包含兩個RCT研究，分別來自Bornstein[14]和Johnson[15]，詳見表3。
表2
IPD部分數據表
表3
AD數據表
2.2.2 數據錄入
本文數據無需錄入僅加載crossnma程序包后就可使用，可使用代碼head()查看數據集：head(ipddata)或head(stddata)。
2.3 數據運算
2.3.1 建立JAGS模型
在完成數據載入后，將可以采用crossnma程序包中的crossnma.model()函數來構建JAGS軟件所需的模型。本處展示未調整與調整模型構建，后續結果和圖形展示均以未調整結果演示為主。具體構建模型代碼，如下：
未調整模型：mod0 <- crossnma.model(treat, id, relapse, n, design, prt.data = ipddata, std.data = stddata, reference = "A", trt.effect = "random", method.bias = "naive")
調整模型Ⅰ：mod1 <- crossnma.model(treat, id, relapse, n, design, prt.data = ipddata, std.data = stddata, reference="A", trt.effect="random", method.bias='adjust1', bias.type = 'add', bias.effect = 'common', bias = rob, unfav = unfavored, bias.group = bias.group)
調整模型Ⅱ：mod2 <- crossnma.model(treat, id, relapse, n, design, prt.data = ipddata, std.data = stddata, reference = "A", trt.effect = "random", method.bias = 'adjust2', bias.type = 'add', bias.effect = 'common', bias = rob, unfav = unfavored, bias.group = bias.group)。
調整模型Ⅰ是將偏差合并入每個研究中，而調整模型Ⅱ則是將偏差調整為兩個分布之和，可詳見本文1.3與1.4方法學部分。在代碼中，需要對method.bias參數進行賦值，研究者只需設置參數為“adjust1”或“adjust2”，即可實現調整模型Ⅰ或調整模型Ⅱ的使用。
上述代碼中，treat、id、relapse及design為IPD與AD數據集所共有的變量。prt.data為IPD數據；std.data為AD數據；reference指示引用何種處理的名稱，如果沒有指定則默認為第一個字母；trt.effect是定義特異性治療效果模型的變量可設置為“random”或“common”，本研究中將未調整模型該參數設置為“random”意為將每個相對處理效果分配為正態分布以便在研究中結合。method.bias則是定義了兩種研究設計RCT與NRS的結合方法，有多個參數可供選擇：參數“naive”可理解為樸素的結合，二者設計上的差異被忽略；“prior”表示使用NRS證據構造RCT研究中相對治療效果的先驗，當證據為高偏差RCT時為“pi.high.rct=’dbeta(10,1)’”，低偏差RCT時為“pi.low.rct=’dbeta(1,10)’”，高偏差NRS時為“pi.high.rct=’dbeta(30,1)’”，低偏差NRS時為“pi.low.rct=’dbeta(1,30)’”；“adjust1”和“adjust2”則是兩個偏差調整方案；當只有一種研究設計類型（RCT或NRS）時，需要指定此參數來指示對研究數據進行樸素結合還是偏差調整的分析；本研究中將未調整模型中該參數設置為“naive”意為不對模型進行調整僅為樸素結合。crossnma.model()函數還可以加入協變量進行分析，協變量可以是數字變量也可以是二分類變量，這對研究者找到異質性有著較大作用。
2.4 模型運算
在數據運算時，crossnma程序包需要使用crossnma()函數來調取JAGS軟件來完成統計運算工作。其工作原理是crossnma()函數將上一環節中已處理好的對象（含模型、數據等信息）傳遞給JAGS軟件進行數據統計。具體代碼為：fit <- crossnma(mod0, n.adapt = 10000, n.iter = 100000, thin = 3, n.chains = 3)。
其中“mod0”為crossnma.model產生的對象，“n.adapt”是MCMC鏈的退火次數默認值為1 000，“n.iter”是MCMC鏈的迭代次數，“thin”是MCMC鏈的間隔取值數，默認值為1，“n.chains”是擬合MCMC鏈數。為后續處理簡便，此處將運算結果命名為“fit”。
2.5 結果輸出
在NMA中，獲取兩兩干預配對比較的網狀結果至關重要。在crossnma程序包中，計算和排列好NMA結果需要使用到league()函數，具體代碼為：league(fit, exp = TRUE, direction = "long")，結果見表4。
表4
網狀Meta分析結果表
2.6 圖形繪制
R軟件crossnma程序包除了能夠計算NMA數據，還可繪制相關功能性圖形，這使得該程序包功能更加全面。當前，crossnma程序包僅可繪制兩種圖形，即網狀關系圖（network plot）與熱力圖（heat plot），這也導致該程序包缺乏部分圖形繪制功能。
2.6.1 網狀關系圖繪制
網狀關系圖是NMA中較為關鍵的圖形之一。在crossnma程序包中，繪制網狀關系圖需要使用netgraph()函數，具體代碼為：netgraph(fit)。
2.6.2 熱力圖繪制
此外，crossnma程序包還提供了熱力圖繪制功能，在crossnma程序包中，熱力圖需要使用heatplot()函數，具體代碼為：heatplot(fit)。
3 討論
本研究從統計學方法與程序包實例操作，分別從模型構建、數據結構分析、軟件參數設置及圖形輸出等方面，對crossnma程序包作全面講解。在基于貝葉斯框架和馬爾可夫鏈蒙特卡羅算法下，crossnma程序包實現了對來自于RCT或NRS跨研究設計的IPD和AD不同數據格式的證據整合，并構建了相應的統計模型，還可通過調整相關參數來將盡可能降低來自于不同研究設計研究間的偏倚風險混入所帶來結果偏差，該功能可在crossnma.model()函數中的method.bias變量中予以實現。基于IPD數據集與AD數據集的合并方法要比單純分析AD數據的NMA具有兩個關鍵性優勢。首先，AD數據研究僅提供相對效應修飾因子與干預措施間的相互作用，而IPD數據研究則考慮到了個體患者水平的相互作用，從而避免了生態偏倚。其次，無論受試者干預分配的情況，IPD數據都會根據預后因素和協變量來進行調整，這些因素和協變量都可幫助預測該疾病結局和病程轉歸。在擬合NMA模型中調整預后因素是不可或缺的，這樣不僅可以提高外部有效性和估計治療效果的精度，并有利于研究結果合理化解釋及糾正隨機化后基線中潛在的不平衡性等問題。
crossnma程序包中調整模型Ⅰ與調整模型Ⅱ對于潛在的偏倚調整的模型構建不同，可借助其他程序包提供的信息偏差準則（deviance information criterion，DIC）來輔助診斷模型與數據鍥合度問題，以幫助操作者做出不同模型選擇決策。此外，若調整模型與未調整模型結果存在較大差距，則說明該類型數據的合并可能受到不同數據格式的偏倚影響較大，建議研究者考慮采用亞組分析或其他合并方法實現合并以減小證據合并帶來的偏倚。
與R軟件gemtc程序包相比，兩個程序包存在一定異同點并各有優劣。盡管，二者均是基于R語言平臺研發，其貝葉斯統計框架均通過JAGS軟件予以實現NMA，然而，它們也存在諸多不同之處。首先，gemtc程序包可應用于連續型變量與二分類變量數據的NMA，而crossnma程序包目前只適用于二分類變量數據NMA。其次，在同時處理不同格式的數據方面，gemtc程序包僅應用于AD數據的合并，而crossnma程序包可同時將AD數據集與IPD數據集進行合并，極大可能增加了NMA樣本量。另外在數據統計與圖形繪制方面，二者在功能上存在較大差距，其中，gemtc程序包具有不一致性檢測、異質性分析、對干預措施進行排序等功能，并能夠提供諸多圖形繪制功能，如網狀關系圖、森林圖、收斂診斷圖、軌跡密度圖等；而crossnma程序包僅能繪制網狀關系圖、軌跡圖等部分圖形，其他圖形的繪制功能仍需進一步研發。
盡管，crossnma程序包中在統計學方法層面上解決了證據定量合成階段的偏差問題，但仍有幾點值得使用者關注。首先，根據以往的經驗證據表明：在高偏倚風險研究中，其干預效果往往會被夸大[18]。一旦出現這種情況時，研究者可采用非正態隨機效應分布降低高偏倚風險研究所帶來結果偏差的影響[19]。其次，在研究設計和解釋結果方面，也應當降低高風險研究（如，NRS）所帶來的偏倚影響。Sarri等[20]提出七步法框架描述了如何將RCT和NRS數據集進行合并，這一框架對NRS進行一定批判性評估，這一觀點或許可通過偏差調整模型予以實現。最后，crossnma程序包在NMA數據統計和圖形繪制等方面仍存在不足，如缺失不一致性檢測、干預排序計算、森林圖及漏斗圖繪制等。
綜上所述，R語言crossnma程序包為NMA的跨研究設計和不同格式數據集的合并提供可能，同時，該程序包通過模型校正因子盡可能降低不同研究設計數據集偏倚對結果所帶來的影響。目前，該程序包在NMA干預排序和常規圖形繪制等方面功能仍存在不足和缺失，需進一步研發和完善。
網狀Meta分析（network meta-analysis，NMA）是一種通過綜合多個臨床研究數據且比較多種干預措施的療效和安全性的統計技術，其可為證據網絡中所有干預方案提供優劣排序結果，并可為臨床決策提供直接證據支撐。通常，NMA在數據類型的選取時，僅局限于隨機對照試驗（randomized controlled trials，RCT），是為了降低不確定性混雜因素介入，以保證網狀數據結果準確性和穩定性。隨著統計學家在方法學上不斷更新，現已可實現將不同研究設計和數據格式的數據集融入到同一個NMA的框架中，以盡可能達到更大樣本結果的統計學效能。
2022年，Tasnim教授團隊通過網狀Meta回歸（network meta-regression，NMR）模型來實現跨研究設計和跨數據格式的數據集（cross-design evidence and cross-format data）合并，并基于R語言平臺開發了crossnma程序包以便研究者使用。NMR模型不僅可將RCT和非隨機研究（nonrandomized studies，NRS）進行融合，還可將聚合數據（aggregate data，AD）和個體參與者數據（individual participant data，IPD）類型數據融入同一NMA當中。此外，NMR模型還考慮到跨研究設計和跨數據格式的數據集間混雜因素和協變量所帶來的潛在異質性等問題。例如：對在已發表的文獻中提取出的AD數據而言，其協變量需要基于整個研究層面水平進行，而IPD數據，則需要從單個樣本病例水平來予以考量和設計，以避免兩個類型數據偏倚和混雜因素相互干擾。
NMA中的證據通常來自RCT，盡管人們認為RCT與NRS相比具有較低的偏倚風險，但Cochrane的一項系統評價發現RCT與NRS在提供不同干預效果估計上并無實質性差異[1]，此外，實證研究證實RCT中可能存在不同類型的偏倚。Schulz等[2]發現分配隱藏不足或缺乏盲法的RCT更傾向于夸大估計的治療效果，并會提供有偏的主觀結果。Chalmers等[3]在非盲或隨機化不規范的研究中發現干預組與對照組間的療效存在較大差異。與此同時，在分配隱藏較差或缺乏盲法的研究中，Wood等[4]發現相對治療效果往往被評估者主觀判斷或患者自我報告夸大。這些方面結論為RCT與NRS跨研究設計數據類型的合并提供了方法學基礎。
前期研究通常采用匹配調整間接比較[5]和模擬治療比較方法[6]通過二次分配權重法和回歸模型法來調整來自于IPD與AD數據集的效應修正因子。然而，這種效應修正因子的調整方法需要對每個干預措施來優化校正，因此，需要每個干預措施比較中至少有一個IPD數據。Jansen等[7]將三級分層模型擴展到標準NMA數據模型當中，該方法通過整合AD數據集潛在的IPD分布來實現將二分類結局數據擴展到其他類型數據，從而達到了IPD與AD差異化數據合并[8]。
考慮到不同偏倚風險（risk of bias，RoB）差異，crossnma程序包采用了Dias等[9]開發的模型將每個研究的RoB因子考慮到模型當中，以降低跨研究設計間數據合并可能帶來的偏差。本研究將通過4個IPD數據集[10-13]與2個AD數據集[14-15]的實例對crossnma程序包模型方法和軟件操作進行實例演示和講解。
1 R軟件crossnma程序包方法學簡介
R軟件crossnma程序包根據IPD網狀數據與AD網狀數據特征，運用網狀Meta回歸模型分別擬定對應的回歸模型，并將兩個數據模型結果進行整合。
1.1 針對IPD數據的未調整網狀Meta回歸模型
|
公式（1）中，服從伯努利分布；為在第j個IPD數據中，參考干預b的效應量；為獨立效應，為協變量，與可在兩式中相抵消，此時，式中剩下相對效應、組內協變量校正因子（組內協變量系數與干預j的平均協變量的乘積）和組間與組內協變量校正因子差值。
1.2 針對AD數據的未調整網狀Meta回歸模型
|
公式（2）中，為研究中每臂的事件數服從二項式分布，此時模型較為簡單，對照組的概率僅為隨機擾動項，它可在二者中相抵消，則關于AD數據相對效應的Meta回歸中自變量只有平均協變量。
此時作為似然函數是隨機擾動項，服從，則先驗服從，對應的貝葉斯公式為公式（3）。
|
那么先驗從何得來，本研究僅使用NRS研究得到均值為，方差為的后驗分布，再將這個后驗分布作為RCT的先驗，將均值加上，方差除以一個0到1的膨脹因子，經過如此調整就可以得到先驗。本研究并沒有使用來自NRS研究的AD，而是使用了RCT研究的AD，在程序包中可設置參數將RCT研究的AD視為NRS數據來參與計算。
1.3 針對IPD數據的風險調整網狀Meta回歸模型
在偏差調整模型中作者使用了兩種不同的方法。偏差調整模型I［公式（4）］是將偏差使用加法效應和乘法效應合并到每一個研究，并用二元參數（0為沒有偏差，1為存在偏差）來衡量是否存在偏差，參數服從零一分布即，其中，與控制著分布的偏斜程度。當，說明分布沒有偏斜其RoB為未知；當，此時分布有著高偏倚或低偏倚，在本程序包中面對不同偏倚風險的證據均可設置不同的參數來相匹配。為乘法偏差效應系數，為加法偏差效應系數。
|
而偏差調整模型Ⅱ［公式（5）］則是使用一個服從雙峰正態分布的參數來代替，服從隨機效應的混合正態分布，并將偏倚效應納入其中，此時偏差調整變為兩個分布之和。
|
1.4 針對AD數據的風險調整網狀Meta回歸模型
|
相較于IPD數據的風險調整模型，AD數據的風險調整模型［公式（6）］略顯簡單。在估算方面還有另一種計算方法，使用對數變換后的再加上隨機擾動項，即。對于偏差效應有著兩種假設，當在隨機效應中，而在固定效應中。為平均偏差效應，當b為非活性藥物干預措施時，當b和k為活性藥物干預措施時；取0時則為零均值偏差效應，取為共同偏差效應；其中為偏差方向，當取值為0時則結果偏向于b，當取值為1時結果偏向于k。可由參數bias.group控制，當其取0時為無偏差調整，取1時為非活性藥物干預措施，取2時為活性藥物干預措施。
|
還有另一種加法偏差效應的運用［公式（7）］，此時原來的OR指數已變為加入偏差效應的即，此時。偏差調整模型Ⅱ［公式（8）］更是大同小異，由此假設一個共同模型來總結相對效應。
|
2 R軟件crossnma程序包操作簡介
2.1 程序包下載與安裝
本研究使用的R 4.2.2，具體安裝方式可參閱《R軟件Metafor程序包在Meta分析中的應用》[16]。值得注意的是，crossnma程序包的運算是基于AJGS軟件來進行數據處理的，因此，需要事先安裝JAGS軟件。目前，JAGS軟件可免費從官方網站下載，具體下載方式請參考《R軟件調用JAGS軟件實現網狀Meta分析》[17]一文。
在使用R軟件crossnma程序包來執行IPD和AD數據NMA時，需要安裝和加載其他輔助性程序包，如：netmeta與meta。具體安裝和加載代碼，見表1。
表1
程序安裝與加載代碼表
2.2 數據介紹與預處理
2.2.1 數據介紹
復發緩解型多發性硬化癥（relapsing remitting multiple sclerosis，RRMS）為多發性硬化癥（multiple sclerosis，MS）患病人數最多的一種臨床亞型，其臨床表現為神經系統退行性慢性炎癥，由于其預后較差、人數較多故而倍受關注。該程序包研發者將以不同的干預措施治療復發緩解型多發性硬化癥后2年內是否復發的IPD和AD兩種不同數據格式臨床數據集為實例。其中，IPD數據集為個體受試者數據，含1 945例數據，4種干預措施（A、B、C和D），而AD數據集僅有4例數據，兩種干預措施（A和C）。本研究IPD數據集其中包含一個NRS研究（Swiss Multiple Sclerosis Cohort[10]）和三個RCT研究（分別來自CONFIRM[11]、DEFINE[12]和AFFIRM[13]），詳見表2。AD數據集包含兩個RCT研究，分別來自Bornstein[14]和Johnson[15]，詳見表3。
表2
IPD部分數據表
表3
AD數據表
2.2.2 數據錄入
本文數據無需錄入僅加載crossnma程序包后就可使用，可使用代碼head()查看數據集：head(ipddata)或head(stddata)。
2.3 數據運算
2.3.1 建立JAGS模型
在完成數據載入后，將可以采用crossnma程序包中的crossnma.model()函數來構建JAGS軟件所需的模型。本處展示未調整與調整模型構建，后續結果和圖形展示均以未調整結果演示為主。具體構建模型代碼，如下：
未調整模型：mod0 <- crossnma.model(treat, id, relapse, n, design, prt.data = ipddata, std.data = stddata, reference = "A", trt.effect = "random", method.bias = "naive")
調整模型Ⅰ：mod1 <- crossnma.model(treat, id, relapse, n, design, prt.data = ipddata, std.data = stddata, reference="A", trt.effect="random", method.bias='adjust1', bias.type = 'add', bias.effect = 'common', bias = rob, unfav = unfavored, bias.group = bias.group)
調整模型Ⅱ：mod2 <- crossnma.model(treat, id, relapse, n, design, prt.data = ipddata, std.data = stddata, reference = "A", trt.effect = "random", method.bias = 'adjust2', bias.type = 'add', bias.effect = 'common', bias = rob, unfav = unfavored, bias.group = bias.group)。
調整模型Ⅰ是將偏差合并入每個研究中，而調整模型Ⅱ則是將偏差調整為兩個分布之和，可詳見本文1.3與1.4方法學部分。在代碼中，需要對method.bias參數進行賦值，研究者只需設置參數為“adjust1”或“adjust2”，即可實現調整模型Ⅰ或調整模型Ⅱ的使用。
上述代碼中，treat、id、relapse及design為IPD與AD數據集所共有的變量。prt.data為IPD數據；std.data為AD數據；reference指示引用何種處理的名稱，如果沒有指定則默認為第一個字母；trt.effect是定義特異性治療效果模型的變量可設置為“random”或“common”，本研究中將未調整模型該參數設置為“random”意為將每個相對處理效果分配為正態分布以便在研究中結合。method.bias則是定義了兩種研究設計RCT與NRS的結合方法，有多個參數可供選擇：參數“naive”可理解為樸素的結合，二者設計上的差異被忽略；“prior”表示使用NRS證據構造RCT研究中相對治療效果的先驗，當證據為高偏差RCT時為“pi.high.rct=’dbeta(10,1)’”，低偏差RCT時為“pi.low.rct=’dbeta(1,10)’”，高偏差NRS時為“pi.high.rct=’dbeta(30,1)’”，低偏差NRS時為“pi.low.rct=’dbeta(1,30)’”；“adjust1”和“adjust2”則是兩個偏差調整方案；當只有一種研究設計類型（RCT或NRS）時，需要指定此參數來指示對研究數據進行樸素結合還是偏差調整的分析；本研究中將未調整模型中該參數設置為“naive”意為不對模型進行調整僅為樸素結合。crossnma.model()函數還可以加入協變量進行分析，協變量可以是數字變量也可以是二分類變量，這對研究者找到異質性有著較大作用。
2.4 模型運算
在數據運算時，crossnma程序包需要使用crossnma()函數來調取JAGS軟件來完成統計運算工作。其工作原理是crossnma()函數將上一環節中已處理好的對象（含模型、數據等信息）傳遞給JAGS軟件進行數據統計。具體代碼為：fit <- crossnma(mod0, n.adapt = 10000, n.iter = 100000, thin = 3, n.chains = 3)。
其中“mod0”為crossnma.model產生的對象，“n.adapt”是MCMC鏈的退火次數默認值為1 000，“n.iter”是MCMC鏈的迭代次數，“thin”是MCMC鏈的間隔取值數，默認值為1，“n.chains”是擬合MCMC鏈數。為后續處理簡便，此處將運算結果命名為“fit”。
2.5 結果輸出
在NMA中，獲取兩兩干預配對比較的網狀結果至關重要。在crossnma程序包中，計算和排列好NMA結果需要使用到league()函數，具體代碼為：league(fit, exp = TRUE, direction = "long")，結果見表4。
表4
網狀Meta分析結果表
2.6 圖形繪制
R軟件crossnma程序包除了能夠計算NMA數據，還可繪制相關功能性圖形，這使得該程序包功能更加全面。當前，crossnma程序包僅可繪制兩種圖形，即網狀關系圖（network plot）與熱力圖（heat plot），這也導致該程序包缺乏部分圖形繪制功能。
2.6.1 網狀關系圖繪制
網狀關系圖是NMA中較為關鍵的圖形之一。在crossnma程序包中，繪制網狀關系圖需要使用netgraph()函數，具體代碼為：netgraph(fit)。
2.6.2 熱力圖繪制
此外，crossnma程序包還提供了熱力圖繪制功能，在crossnma程序包中，熱力圖需要使用heatplot()函數，具體代碼為：heatplot(fit)。
3 討論
本研究從統計學方法與程序包實例操作，分別從模型構建、數據結構分析、軟件參數設置及圖形輸出等方面，對crossnma程序包作全面講解。在基于貝葉斯框架和馬爾可夫鏈蒙特卡羅算法下，crossnma程序包實現了對來自于RCT或NRS跨研究設計的IPD和AD不同數據格式的證據整合，并構建了相應的統計模型，還可通過調整相關參數來將盡可能降低來自于不同研究設計研究間的偏倚風險混入所帶來結果偏差，該功能可在crossnma.model()函數中的method.bias變量中予以實現。基于IPD數據集與AD數據集的合并方法要比單純分析AD數據的NMA具有兩個關鍵性優勢。首先，AD數據研究僅提供相對效應修飾因子與干預措施間的相互作用，而IPD數據研究則考慮到了個體患者水平的相互作用，從而避免了生態偏倚。其次，無論受試者干預分配的情況，IPD數據都會根據預后因素和協變量來進行調整，這些因素和協變量都可幫助預測該疾病結局和病程轉歸。在擬合NMA模型中調整預后因素是不可或缺的，這樣不僅可以提高外部有效性和估計治療效果的精度，并有利于研究結果合理化解釋及糾正隨機化后基線中潛在的不平衡性等問題。
crossnma程序包中調整模型Ⅰ與調整模型Ⅱ對于潛在的偏倚調整的模型構建不同，可借助其他程序包提供的信息偏差準則（deviance information criterion，DIC）來輔助診斷模型與數據鍥合度問題，以幫助操作者做出不同模型選擇決策。此外，若調整模型與未調整模型結果存在較大差距，則說明該類型數據的合并可能受到不同數據格式的偏倚影響較大，建議研究者考慮采用亞組分析或其他合并方法實現合并以減小證據合并帶來的偏倚。
與R軟件gemtc程序包相比，兩個程序包存在一定異同點并各有優劣。盡管，二者均是基于R語言平臺研發，其貝葉斯統計框架均通過JAGS軟件予以實現NMA，然而，它們也存在諸多不同之處。首先，gemtc程序包可應用于連續型變量與二分類變量數據的NMA，而crossnma程序包目前只適用于二分類變量數據NMA。其次，在同時處理不同格式的數據方面，gemtc程序包僅應用于AD數據的合并，而crossnma程序包可同時將AD數據集與IPD數據集進行合并，極大可能增加了NMA樣本量。另外在數據統計與圖形繪制方面，二者在功能上存在較大差距，其中，gemtc程序包具有不一致性檢測、異質性分析、對干預措施進行排序等功能，并能夠提供諸多圖形繪制功能，如網狀關系圖、森林圖、收斂診斷圖、軌跡密度圖等；而crossnma程序包僅能繪制網狀關系圖、軌跡圖等部分圖形，其他圖形的繪制功能仍需進一步研發。
盡管，crossnma程序包中在統計學方法層面上解決了證據定量合成階段的偏差問題，但仍有幾點值得使用者關注。首先，根據以往的經驗證據表明：在高偏倚風險研究中，其干預效果往往會被夸大[18]。一旦出現這種情況時，研究者可采用非正態隨機效應分布降低高偏倚風險研究所帶來結果偏差的影響[19]。其次，在研究設計和解釋結果方面，也應當降低高風險研究（如，NRS）所帶來的偏倚影響。Sarri等[20]提出七步法框架描述了如何將RCT和NRS數據集進行合并，這一框架對NRS進行一定批判性評估，這一觀點或許可通過偏差調整模型予以實現。最后，crossnma程序包在NMA數據統計和圖形繪制等方面仍存在不足，如缺失不一致性檢測、干預排序計算、森林圖及漏斗圖繪制等。
綜上所述，R語言crossnma程序包為NMA的跨研究設計和不同格式數據集的合并提供可能，同時，該程序包通過模型校正因子盡可能降低不同研究設計數據集偏倚對結果所帶來的影響。目前，該程序包在NMA干預排序和常規圖形繪制等方面功能仍存在不足和缺失，需進一步研發和完善。
