隨著單臂臨床試驗日益增多，缺乏頭對頭臨床研究，未調整的間接比較和網狀Meta分析應用局限性日益突出。當個體數據數量有限時，匹配調整間接比較能夠充分地利用單項研究個體數據并通過傾向得分匹配其他研究的匯總數據，有效地平衡不同試驗中人群差異帶來的潛在偏倚，并完成目標干預措施間的療效比較。本研究一方面圍繞匹配調整間接比較有關概念和原理進行介紹，另一方面，基于當前腫瘤藥物評價的廣泛應用，系統展示如何基于R語言使用錨定匹配調整間接比較法針對生存數據進行分析，以期為科學循證決策提供參考。
在特定人群中，頭對頭隨機對照試驗（randomized controlled trial，RCT）是無偏估計干預措施間療效差異的金標準。然而在醫療實踐中，由于創新治療技術快速上市，頭對頭試驗較為缺乏，導致決策所需直接比較證據不足。雖然間接比較和網狀Meta分析（network meta-analysis，NMA）可在一定程度上解決上述問題，但是方法本身未考慮人群基線特征差異，簡單地合并匯總數據（aggregated data，AgD），例如風險比（hazard ratio，HR）等效應量，可能產生偏倚較大的研究結果[1-2]。此外，間接比較和NMA應用的局限性明顯，當缺少共同對照時，難以構建合適“橋梁”完成對不同干預措施的比較[3]。
理想情況下，當能夠完全獲得擬對比干預措施RCT個體數據（individual patient data，IPD）時，可以通過傾向得分匹配調整（propensity score matching，PSM）法對人群間的基線特征進行調整。然而PSM法會受限于IPD的可獲得性，如多數情況下研究者僅能獲得某一研究的IPD和另一研究的AgD，無法開展PSM。為解決該問題，匹配調整間接比較（matching-adjusted indirect comparisons，MAIC）方法得以開發和使用[4]。
MAIC是一種基于對人群傾向得分加權匹配的間接比較方法，目前已廣泛應用于英國國家衛生與臨床優化研究所（National Institute for Health and Care Excellence，NICE）的衛生技術評估研究中[5]。MAIC能夠依據AgD的人群基本特征對IPD校準，然后完成比較分析。根據研究間有無共同對照，可將MAIC分為有共同對照組的錨定MAIC和無共同對照的非錨定MAIC[6]。由于錨定MAIC法有共同對照組作為橋梁，統計方法上認為錨定MAIC可能給出更穩健的結果，因此更推薦錨定MAIC法[7]。
因為錨定MAIC和非錨定MAIC法的原理和應用類似，為方便說明，本研究主要以錨定MAIC方法為例，基于Phillipo等[6]開發的理論和代碼，詳細闡述MAIC的原理，同時考慮到目前腫瘤藥物評價的廣泛開展，本研究將具體介紹如何基于R語言使用錨定MAIC法對生存數據進行分析。
1 方法簡介
錨定MAIC法是使得試驗AB人群（IPD資料，接受干預措施A或干預措施B）加權匹配調整后與試驗AC人群（AgD資料，接受干預措施A或干預措施C）可比。相應的，調整后試驗AB中干預措施t（干預措施t=A或B）的平均結果產出
可用如下公式表達：
 |
公式1中，
代表試驗AB中接受干預措施t的人數；
代表試驗AB中接受干預措施t的第i個患者被納入試驗AC人群結構中的權重。
1.1 明確效應修飾因子或預后變量
在加權匹配調整前，為降低匹配工作量并提高匹配精確度，首先需要根據結局變量特征建立相應的多變量回歸模型，然后依據顯著性水平對相關協變量進行判斷，以明確其是否為重要的潛在效應修飾因子或預后變量。
1.2 估計權重
權重的估計是按照試驗AC人群基本特征使用試驗AB人群數據估計傾向得分Logistic回歸模型，進而得到匹配調整權重。（公式2）
 |
公式2中，
是常數項，
是干預措施t的回歸參數，
是試驗AB中接受干預措施t的第i個患者的協變量。因此將公式2左右兩邊轉化為指數形式并代入公式1中，可得到公式3。
 |
一般情況下，如果有足夠的IPD資料，可以使用類似于PSM中應用的標準方法估計回歸參數。在只有試驗AB的IPD資料，且無法充分獲取試驗AC中聯合協變量分布的情況下，可采用Signorovitch等[8]提出的矩估計法。在矩估計法中，為確保權重能平衡好匹配調整后的試驗AB總體與試驗AC總體之間的平均協變量值，這相當于要求實現最小化公式4，且此時試驗AC的平均協變量值
。
 |
將公式4作為目標函數，為求其最小值時對應的
，可以通過構建對應的梯度（導數）函數（公式5），然后使用BFGS算法進行迭代獲得，進而通過公式4計算得到匹配調整權重。然后按照公式1對IPD校準并完成比較分析。
 |
2 軟件安裝與程序加載
2.1 R語言安裝
截至筆者完稿時，R軟件最新版本為4.2.2，讀者可在其官網下載與計算機操作系統相匹配的版本。
2.2 程序包加載
本研究中需要安裝的程序包包括：dplyr、tidyr、boot、survival、survminer、epiDisplay、ggplot2和ggpubr。通過install.packages函數安裝成功后，可通過library函數加載調用。以dplyr程序包的安裝和加載為例，具體代碼如下：（注：代碼中“#”后面的文字是對代碼的解釋，不參與代碼運行）
install.packages("dplyr") #安裝程序包
library("dplyr") #加載程序包
3 數據加載與預處理
3.1 數據的讀取與錄入
本研究中IPD資料采用survival程序包中內置的myeloid數據，作為試驗AB數據（干預措施A vs. 干預措施B）；AgD資料是筆者擬定的虛擬數據，當作試驗AC數據（干預措施A vs. 干預措施C）。基于該數據集進行演示，目的是通過錨定MAIC方法，使試驗AB數據經調整后與試驗AC數據基線可比，從而獲得匹配調整后干預措施B相較于干預措施C的HR值。myeloid數據是一項基于急性髓細胞白血病試驗的模擬數據集，該數據集一共包含8個變量，646條數據。變量分別是受試者編號（id）、治療組別（trt）、性別（sex）、生存時間（futime）、事件狀態（death）、造血干細胞移植時間（txtime）、響應時間（crtime）和復發時間（rltime），其中sex、txtime、crtime和rltime可能是影響結局的效應修飾因子或預后變量，見表1（由于數據過多，本研究僅展示前十行數據）。相應的，AgD資料筆者擬定男性比例為52%，造血干細胞平均移植時間為189.72，HRAC值為0.62［95%CI（0.54，0.71）］。具體代碼如下：
表1
myeloid數據示例前十行數據展示
AB.IPD<-myeloid #調用myeloid數據
AC.AgD<-data.frame(male.prop_AC=0.52,txtime.mean_AC=189.72,HR_AC=0.62,lHR_AC=0.54,uHR_AC=0.71) #擬定虛擬數據
如果讀者使用其他來源的數據，可將數據記錄在Excel文件中，并以CSV格式存儲，然后通過“read.csv”命令讀取相關數據。
3.2 數據預處理
為簡化匹配調整過程，只將sex和txtime兩個變量當作需調整的效應修飾因子或預后變量，因此，本研究最終用到的IPD資料只包含trt、sex、futime、death和txtime變量，在R語言中可通過select函數實現變量的選擇。同時使用drop_na函數將缺失數據剔除，因為標準MAIC無法處理含缺失值的數據。最后納入分析樣本量是364。具體代碼如下：
AB.IPD<-AB.IPD %>%
dplyr::select(trt,sex,futime,death,txtime) %>% #選擇目標變量
drop_na() %>% #剔除缺失值
mutate(sex=ifelse(sex=="m",0,1)) #數據轉換
然后通過summarise函數查看當前試驗AB中數據匯總情況，并將其與試驗AC中AgD資料相比較，結果見表2，對比可發現兩組試驗人群基本特征存在一定差異。具體代碼如下：
表2
試驗AB和試驗AC中人群基本特征
AB.IPD %>%
summarise(sex_AB=round(prop.table(table(sex)),2)[1], #計算男性比例，并保留2位小數
txtime_AB=round(mean(txtime),2)) #計算造血干細胞移植時間平均值，并保留2位小數
4 錨定MAIC方法實現
4.1 構建目標函數和梯度函數
在R語言中，可以通過function函數創建公式4和公式5。具體代碼如下：
objfn <- function(a1, X){sum(exp(X %*% a1))} #創建目標函數
gradfn <- function(a1, X){colSums(sweep(X, 1, exp(X %*% a1), "*"))} #創建梯度函數
為使得
，在試驗AB和試驗AC中均同時減去
，從而使得試驗AC中平均協變量影響為0，試驗AB中平均協變量影響為
。在R語言中可通過sweep函數實現上述操作。具體代碼如下：
AB.IPD_EM<-sweep(with(AB.IPD, cbind(sex, txtime)), 2,
with(AC.AgD, c(male.prop_AC, txtime.mean_AC)), '-')
然后通過optim函數求解目標函數最小化時的
，結果如圖1，具體代碼如下：
圖1
目標函數最小時求解結果
opt <- optim(par = c(0,0), #指定初始值
fn = objfn, #指定目標函數
gr = gradfn, #指定梯度函數
X = AB.IPD_EM, #指定所需的
時的額外數據
method = "BFGS") #選擇BFGS法作為優化方法
opt
該結果包含了目標函數取值最小時的
，也就是$par的值；$value的值是目標函數的最小值；$counts，是指收斂前目標函數和梯度函數的評估的次數；$convergence=0，表明已成功完成收斂。
4.2 計算權重
根據公式4可計算試驗AB中每個個體在試驗AC人群結構中的權重，具體代碼如下：
a1 <- opt$par
wt <- exp(AB.IPD_EM %*% a1) #計算權重
此外，還可以計算相對權重。如果相對權重大于1，意味著該個體在重新匹配的總體中比在原試驗AB總體中占更多的權重，反之，則占更少的權重。計算過程如公式6所示，具體代碼如下：
 |
wt.rs <- (wt / sum(wt)) * nrow(AB.IPD_EM) #計算相對權重
然后通過summary函數查看每個個體權重情況，同時可通過qplot函數將權重分布情況以直方圖形式呈現，如圖2。具體代碼如下：
圖2
權重分布直方圖
summary(wt)
summary(wt.rs)
plot1<-ggplot2::qplot (wt, #選擇數據集
geom="histogram", #指定使用直方圖
xlab = "weight", #設置X軸標簽
binwidth=0.25) #指定直方圖條形寬度
plot2<- ggplot2::qplot(wt.rs, geom="histogram",xlab = "Rescaled weight",binwidth=0.25)
ggarrange(plot1,plot2, #將圖組合在一張畫布上
ncol = 1, #指定組合后圖形列數
nrow = 2, #指定組合后圖形行數
heights = c(1,1)) #指定每一行的高度比例
4.3 計算近似有效樣本量
Signorovitch等[8]表明通過加權調整后的試驗AB總體形成的偽總體的有效樣本量（effective sample size，ESS）可近似表示為：
 |
計算結果約為349，與納入分析的364條數據相比損失了部分樣本。在R語言中，相應代碼如下：
ESS<-sum(wt)^2/sum(wt^2) #計算近似有效樣本量
ESS
4.4 查看匹配調整后結果
為方便查找和分析數據，可借助cbind函數，將AB.IPD、AB.IPD_EM、wt、wt.rs四個數據框按列合并成一個新的數據框combine_data_AB，具體代碼如下：
combine_data_AB<-cbind(AB.IPD,AB.IPD_EM,wt,wt.rs) #數據框按列合并
names(combine_data_AB)[6:7]<-c("sex_centred","txtime_centred") #對變量重新命名
然后通過summarise函數計算匹配調整后的試驗AB中男性的加權比例和造血干細胞移植加權平均時間。
combine_data_AB %>%
summarise(male.prop_AB=weighted.mean(sex,wt),
txtime.mean_AB=weighted.mean(txtime,wt))
將試驗AB調整前與調整后結果分別和試驗AC中基線情況相比較，比較結果如表3所示，在以試驗AC人群特征為基礎調整后的試驗AB結果與試驗AC中人群特征基本一致。
表3
調整前后試驗AB和試驗AC中人群基本特征
4.5 計算HR值
首先通過coxph函數進行Cox比例風險模型分析，并基于參數weights的設置，對匹配調整前后的試驗AB數據分別構建模型，模型結果可通過summary函數呈現。具體代碼如下：
#構建匹配調整后的Cox比例風險模型
cox_maic <- coxph(Surv(futime, death) ~ trt, #構建模型關系
data=combine_data_AB, #指定數據集
weights=wt) #加權調整
summary(cox_maic, conf.int=0.95) #給定置信水平下結果
#構建未經匹配調整的Cox比例風險模型
cox_unmaic<-coxph(Surv(futime, death) ~ trt,
data=combine_data_AB)
summary(cox_unmaic, conf.int=0.95)
為提高分析準確性，將HR值及對應95%上下限取對數，然后根據Bucher間接比較方法，計算得到干預措施B相對于干預措施C的HRBC值及對應95%上下限。因篇幅有限，本研究僅呈現匹配調整后的HR值和間接比較結果，如果想獲得未調整的結果，可將下列代碼中調整后參數替換為未調整參數。具體代碼如下：
logHR_AC<-log(AC.AgD$HR_AC) #對HRAC取對數
logHR_se_AC<-(log(AC.AgD$uHR_AC)-log(AC.AgD$lHR_AC))/(2*1.96) #計算log(HRAC)標準誤
logHR_AB_maic<-log(summary(cox_maic)$conf.int[1]) #對HRAB取對數
logHR_se_AB_maic<-(log(summary(cox_maic)$conf.int[4])-log(summary(cox_maic)$conf.int[3]))/(2*1.96) #計算log(HRAB)標準誤
HR_BC_maic<-exp(logHR_AC-logHR_AB_maic) #計算HRBC
HR_var_BC_maic<-logHR_se_AC^2+logHR_se_AB_maic^2 #計算HRBC方差
lHR_BC_maic<-exp(log(HR_BC_maic)-1.96*sqrt(HR_var_BC_maic)) #計算HRBC95%下限
uHR_BC_maic<-exp(log(HR_BC_maic)+1.96*sqrt(HR_var_BC_maic)) #計算HRBC95%上限
4.6 Bootstrapping法計算HR值
在實際情況中，通常無法考慮全部影響人群之間平衡的效應修飾因子和預后變量，往往只會對經統計方法檢驗顯著的影響變量進行調整；且不同觀察值可能被賦予不同權重，進而導致結果可能受到觀察值之間相關性的影響，直接使用絕對結果將產生偏倚。因此，需要考慮這些可能存在的偏倚，常用的方法有魯棒方差估計（robust variance estimation，RVE）和Bootstrapping法[9]。其中，Austin等[9]研究發現，使用Bootstrapping法可以通過重復抽樣較為正確地估計標準誤和置信區間。因此本研究使用Bootstrapping法來計算標準誤和置信區間，并通過boot函數在R語言中實現這一操作，具體步驟如下。
首先，通過function函數構建boot函數所需的統計量，具體代碼如下：
hr.fun <- function(data, indices) { #確定函數輸入參數
d <- data[indices,] #根據輸入的索引構建子數據集
fit <- coxph(Surv(futime, death) ~ trt ,weights=wt, data = d) #擬合Cox比例風險模型
ci <- confint(fit) #計算置信區間
c(exp(coef(fit)), exp(ci))} #將模型估計的HR值和置信區間以向量形式返回
然后通過boot函數完成Bootstrapping估計，具體代碼如下：
set.seed(123) #設定種子，確保結果可重現
HR_AB_bootstraps <- boot(data = combine_data_AB, #指定數據
statistic = hr.fun, #指定統計量
R=1000) #指定抽樣次數
然后從Bootstrapping結果中將HRAB值及95%上下限選擇出來并取對數和計算標準誤。具體代碼如下：
logHR_AB_boot <- log(median(HR_AB_bootstraps$t)) #對重復抽樣后中位HRAB取對數
boot_HR_ci_AB_HR <- boot.ci(boot.out = HR_AB_bootstraps, type="perc") #使用百分位數法計算置信區間
lHR_AB_boot<-boot_HR_ci_AB_HR$percent[4] #選擇百分位數法計算置信區間下限
uHR_AB_boot<-boot_HR_ci_AB_HR$percent[5] #選擇百分位數法計算置信區間上限
logHR_se_AB_boot<-(log(lHR_AB_boot)-log(uHR_AB_boot))/(2*1.96) #計算重復抽樣后中位HRAB標準誤
然后同樣根據Bucher法，計算得到更為精確的干預措施B相對于干預措施C的HRBC值及對應95%上下限（代碼類同本文“4.5 計算HR值”）。
最后，用data.frame函數，將未調整、匹配調整后以及Bootstrapping法精確估計計算間接比較結果匯總在一張表，結果如表4所示，具體代碼如下：
表4
不同方法下HRBC及95%置信區間
HR_BC<-c(HR_BC_unmaic,HR_BC_maic,HR_BC_boot)
lHR_BC<-c(lHR_BC_unmaic,lHR_BC_maic,lHR_BC_boot)
uHR_BC<-c(uHR_BC_unmaic,uHR_BC_maic,uHR_BC_boot)
method<-c("unmaic", "maic","boot")
outcome<-data.frame(method,HR_BC,lHR_BC,uHR_BC) #結果以數據框呈現
outcome
從表4結果中可以看出，與干預措施C相比，干預措施B降低了死亡風險，并且匹配調整改變了效應大小，符合我們對潛在效應修飾因子或預后變量（男性比例、造血干細胞平均移植時間）調整后預期的評估；同時，經過Bootstrapping法計算得到的HR值與經基礎加權得到的結果基本一致。
以上就是在R語言中實現基于錨定MAIC法的具體操作過程。如果缺少共同對照，則可以選擇使用非錨定MAIC，其匹配調整的過程與錨定的MAIC過程基本相似，只需將本研究中試驗AB和試驗AC的數據分別當做干預措施B（IPD資料）和干預措施C（AgD資料）的數據，從而獲得干預措施C人群結構下匹配調整后的干預措施B加權數據。由于缺乏共同對照，非錨定MAIC不能直接計算加權調整后干預措施B與干預措施C的相對效應，但可通過數字化工具提取干預措施C生存曲線數據，并重構干預措施C的IPD，然后與匹配調整后的干預措施B數據一同結合在Cox比例風險模型中，進而求得干預措施B與干預措施C的相對效應。有關重構IPD及生存分析的方法已有文獻詳細介紹了在R語言中實現的方法，可參考石豐豪等[10]的研究。
5 討論
本研究對MAIC原理和方法進行介紹，并結合survival內置的myeloid生存數據和虛擬AgD資料，展示了如何在R語言中實現錨定MAIC過程。結果顯示，匹配調整后改變了間接比較結果，這也證明了效應修飾因子或預后變量對結局指標的影響。
隨著依托單臂試驗附條件上市的創新藥物不斷涌現，以及交叉試驗的存在，如何有效評估不同臨床試驗之間的效應差異是一個需要重視的問題。相較于傳統的間接比較，在依托IPD資料的情況下使用MAIC，能夠充分利用已有數據，通過匹配調整有效地平衡好RCT間效應修飾因子的影響，降低比較結果的偏倚風險。鑒于MAIC的適應性，其應用也越來越廣泛。例如，在生存數據中，Chen等[11]在鱗狀非小細胞肺癌一線治療研究中，Maloney等[11]在復發或難治性大B細胞淋巴瘤研究中分別使用了基于錨定MAIC和非錨定MAIC法；在其他類型數據中（二分類或連續性數據），Mt-Isa等[13]采用錨定MAIC法比較不同肺炎球菌疫苗的免疫應答情況，Thom等[14]使用非錨定MAIC法比較克立硼羅軟膏與外用鈣調磷酸酶抑制劑治療輕中度特應性皮炎患者的療效。
MAIC的應用也有一定的局限性[2,7,15]：① 只能對比兩組情況，不適用于復雜證據體；② 默認將提供匯總資料的人群視作決策目標人群，其外推性有限，無法在指定的目標人群中進行效果預測；③ 調整多種基線因素的能力取決于擁有IPD的試驗中是否有足夠數量的患者，且匹配過程中可能存在數據丟失的問題，需要較大樣本量；④ 與傳統傾向得分匹配加權不同，MAIC中AgD資料的可用性阻礙了使用現有方法來檢查傾向評分模型的擬合和校準情況；⑤ 對生存數據進行分析時，更適合應用于等比例假設條件下。
此外，與MAIC法相類似的人群匹配調整方法還有模擬治療比較（simulated treatment comparison，STC）和多水平網狀Meta回歸（multi-level network meta-regression，ML-NMR）。STC與MAIC在適應范圍上無優劣之分，但STC法處理非線性模型的結果（例如Logit回歸或事件時間模型）時，可能會增加偏倚[16]，Veroniki等[17]研究也發現MAIC法應用的比例高于STC法；ML-NMR法雖然相較于MAIC法能夠實現包含多個研究的比較，但是其應用仍只局限于二分類結局數據，此外，能否將單臂研究整合到證據體中也需進一步探索驗證[2]。
綜上所述，在缺乏直接比較證據的前提下，MAIC是一種有效的替代手段，可以為相關研究和決策提供一定的證據支持。此外，Gregory等[18]已經開發出在R語言中進行MAIC的集成程序包MAIC，但該程序包中函數及處理過程多是集成式，在一定程度上簡化了語句，但在理解上可能會存在一定的難度，尤其是對于MAIC處理過程和相關原理的理解方面，讀者可在理解本研究的基礎上進一步探索學習。
聲明 本研究不存在任何利益沖突。
在特定人群中，頭對頭隨機對照試驗（randomized controlled trial，RCT）是無偏估計干預措施間療效差異的金標準。然而在醫療實踐中，由于創新治療技術快速上市，頭對頭試驗較為缺乏，導致決策所需直接比較證據不足。雖然間接比較和網狀Meta分析（network meta-analysis，NMA）可在一定程度上解決上述問題，但是方法本身未考慮人群基線特征差異，簡單地合并匯總數據（aggregated data，AgD），例如風險比（hazard ratio，HR）等效應量，可能產生偏倚較大的研究結果[1-2]。此外，間接比較和NMA應用的局限性明顯，當缺少共同對照時，難以構建合適“橋梁”完成對不同干預措施的比較[3]。
理想情況下，當能夠完全獲得擬對比干預措施RCT個體數據（individual patient data，IPD）時，可以通過傾向得分匹配調整（propensity score matching，PSM）法對人群間的基線特征進行調整。然而PSM法會受限于IPD的可獲得性，如多數情況下研究者僅能獲得某一研究的IPD和另一研究的AgD，無法開展PSM。為解決該問題，匹配調整間接比較（matching-adjusted indirect comparisons，MAIC）方法得以開發和使用[4]。
MAIC是一種基于對人群傾向得分加權匹配的間接比較方法，目前已廣泛應用于英國國家衛生與臨床優化研究所（National Institute for Health and Care Excellence，NICE）的衛生技術評估研究中[5]。MAIC能夠依據AgD的人群基本特征對IPD校準，然后完成比較分析。根據研究間有無共同對照，可將MAIC分為有共同對照組的錨定MAIC和無共同對照的非錨定MAIC[6]。由于錨定MAIC法有共同對照組作為橋梁，統計方法上認為錨定MAIC可能給出更穩健的結果，因此更推薦錨定MAIC法[7]。
因為錨定MAIC和非錨定MAIC法的原理和應用類似，為方便說明，本研究主要以錨定MAIC方法為例，基于Phillipo等[6]開發的理論和代碼，詳細闡述MAIC的原理，同時考慮到目前腫瘤藥物評價的廣泛開展，本研究將具體介紹如何基于R語言使用錨定MAIC法對生存數據進行分析。
1 方法簡介
錨定MAIC法是使得試驗AB人群（IPD資料，接受干預措施A或干預措施B）加權匹配調整后與試驗AC人群（AgD資料，接受干預措施A或干預措施C）可比。相應的，調整后試驗AB中干預措施t（干預措施t=A或B）的平均結果產出可用如下公式表達：
|
公式1中，代表試驗AB中接受干預措施t的人數；代表試驗AB中接受干預措施t的第i個患者被納入試驗AC人群結構中的權重。
1.1 明確效應修飾因子或預后變量
在加權匹配調整前，為降低匹配工作量并提高匹配精確度，首先需要根據結局變量特征建立相應的多變量回歸模型，然后依據顯著性水平對相關協變量進行判斷，以明確其是否為重要的潛在效應修飾因子或預后變量。
1.2 估計權重
權重的估計是按照試驗AC人群基本特征使用試驗AB人群數據估計傾向得分Logistic回歸模型，進而得到匹配調整權重。（公式2）
|
公式2中，是常數項，是干預措施t的回歸參數，是試驗AB中接受干預措施t的第i個患者的協變量。因此將公式2左右兩邊轉化為指數形式并代入公式1中，可得到公式3。
|
一般情況下，如果有足夠的IPD資料，可以使用類似于PSM中應用的標準方法估計回歸參數。在只有試驗AB的IPD資料，且無法充分獲取試驗AC中聯合協變量分布的情況下，可采用Signorovitch等[8]提出的矩估計法。在矩估計法中，為確保權重能平衡好匹配調整后的試驗AB總體與試驗AC總體之間的平均協變量值，這相當于要求實現最小化公式4，且此時試驗AC的平均協變量值。
|
將公式4作為目標函數，為求其最小值時對應的，可以通過構建對應的梯度（導數）函數（公式5），然后使用BFGS算法進行迭代獲得，進而通過公式4計算得到匹配調整權重。然后按照公式1對IPD校準并完成比較分析。
|
2 軟件安裝與程序加載
2.1 R語言安裝
截至筆者完稿時，R軟件最新版本為4.2.2，讀者可在其官網下載與計算機操作系統相匹配的版本。
2.2 程序包加載
本研究中需要安裝的程序包包括：dplyr、tidyr、boot、survival、survminer、epiDisplay、ggplot2和ggpubr。通過install.packages函數安裝成功后，可通過library函數加載調用。以dplyr程序包的安裝和加載為例，具體代碼如下：（注：代碼中“#”后面的文字是對代碼的解釋，不參與代碼運行）
install.packages("dplyr") #安裝程序包
library("dplyr") #加載程序包
3 數據加載與預處理
3.1 數據的讀取與錄入
本研究中IPD資料采用survival程序包中內置的myeloid數據，作為試驗AB數據（干預措施A vs. 干預措施B）；AgD資料是筆者擬定的虛擬數據，當作試驗AC數據（干預措施A vs. 干預措施C）。基于該數據集進行演示，目的是通過錨定MAIC方法，使試驗AB數據經調整后與試驗AC數據基線可比，從而獲得匹配調整后干預措施B相較于干預措施C的HR值。myeloid數據是一項基于急性髓細胞白血病試驗的模擬數據集，該數據集一共包含8個變量，646條數據。變量分別是受試者編號（id）、治療組別（trt）、性別（sex）、生存時間（futime）、事件狀態（death）、造血干細胞移植時間（txtime）、響應時間（crtime）和復發時間（rltime），其中sex、txtime、crtime和rltime可能是影響結局的效應修飾因子或預后變量，見表1（由于數據過多，本研究僅展示前十行數據）。相應的，AgD資料筆者擬定男性比例為52%，造血干細胞平均移植時間為189.72，HRAC值為0.62［95%CI（0.54，0.71）］。具體代碼如下：
表1
myeloid數據示例前十行數據展示
AB.IPD<-myeloid #調用myeloid數據
AC.AgD<-data.frame(male.prop_AC=0.52,txtime.mean_AC=189.72,HR_AC=0.62,lHR_AC=0.54,uHR_AC=0.71) #擬定虛擬數據
如果讀者使用其他來源的數據，可將數據記錄在Excel文件中，并以CSV格式存儲，然后通過“read.csv”命令讀取相關數據。
3.2 數據預處理
為簡化匹配調整過程，只將sex和txtime兩個變量當作需調整的效應修飾因子或預后變量，因此，本研究最終用到的IPD資料只包含trt、sex、futime、death和txtime變量，在R語言中可通過select函數實現變量的選擇。同時使用drop_na函數將缺失數據剔除，因為標準MAIC無法處理含缺失值的數據。最后納入分析樣本量是364。具體代碼如下：
AB.IPD<-AB.IPD %>%
dplyr::select(trt,sex,futime,death,txtime) %>% #選擇目標變量
drop_na() %>% #剔除缺失值
mutate(sex=ifelse(sex=="m",0,1)) #數據轉換
然后通過summarise函數查看當前試驗AB中數據匯總情況，并將其與試驗AC中AgD資料相比較，結果見表2，對比可發現兩組試驗人群基本特征存在一定差異。具體代碼如下：
表2
試驗AB和試驗AC中人群基本特征
AB.IPD %>%
summarise(sex_AB=round(prop.table(table(sex)),2)[1], #計算男性比例，并保留2位小數
txtime_AB=round(mean(txtime),2)) #計算造血干細胞移植時間平均值，并保留2位小數
4 錨定MAIC方法實現
4.1 構建目標函數和梯度函數
在R語言中，可以通過function函數創建公式4和公式5。具體代碼如下：
objfn <- function(a1, X){sum(exp(X %*% a1))} #創建目標函數
gradfn <- function(a1, X){colSums(sweep(X, 1, exp(X %*% a1), "*"))} #創建梯度函數
為使得，在試驗AB和試驗AC中均同時減去，從而使得試驗AC中平均協變量影響為0，試驗AB中平均協變量影響為。在R語言中可通過sweep函數實現上述操作。具體代碼如下：
AB.IPD_EM<-sweep(with(AB.IPD, cbind(sex, txtime)), 2,
with(AC.AgD, c(male.prop_AC, txtime.mean_AC)), '-')
然后通過optim函數求解目標函數最小化時的，結果如圖1，具體代碼如下：
圖1
目標函數最小時求解結果
opt <- optim(par = c(0,0), #指定初始值
fn = objfn, #指定目標函數
gr = gradfn, #指定梯度函數
X = AB.IPD_EM, #指定所需的時的額外數據
method = "BFGS") #選擇BFGS法作為優化方法
opt
該結果包含了目標函數取值最小時的，也就是$par的值；$value的值是目標函數的最小值；$counts，是指收斂前目標函數和梯度函數的評估的次數；$convergence=0，表明已成功完成收斂。
4.2 計算權重
根據公式4可計算試驗AB中每個個體在試驗AC人群結構中的權重，具體代碼如下：
a1 <- opt$par
wt <- exp(AB.IPD_EM %*% a1) #計算權重
此外，還可以計算相對權重。如果相對權重大于1，意味著該個體在重新匹配的總體中比在原試驗AB總體中占更多的權重，反之，則占更少的權重。計算過程如公式6所示，具體代碼如下：
|
wt.rs <- (wt / sum(wt)) * nrow(AB.IPD_EM) #計算相對權重
然后通過summary函數查看每個個體權重情況，同時可通過qplot函數將權重分布情況以直方圖形式呈現，如圖2。具體代碼如下：
圖2
權重分布直方圖
summary(wt)
summary(wt.rs)
plot1<-ggplot2::qplot (wt, #選擇數據集
geom="histogram", #指定使用直方圖
xlab = "weight", #設置X軸標簽
binwidth=0.25) #指定直方圖條形寬度
plot2<- ggplot2::qplot(wt.rs, geom="histogram",xlab = "Rescaled weight",binwidth=0.25)
ggarrange(plot1,plot2, #將圖組合在一張畫布上
ncol = 1, #指定組合后圖形列數
nrow = 2, #指定組合后圖形行數
heights = c(1,1)) #指定每一行的高度比例
4.3 計算近似有效樣本量
Signorovitch等[8]表明通過加權調整后的試驗AB總體形成的偽總體的有效樣本量（effective sample size，ESS）可近似表示為：
|
計算結果約為349，與納入分析的364條數據相比損失了部分樣本。在R語言中，相應代碼如下：
ESS<-sum(wt)^2/sum(wt^2) #計算近似有效樣本量
ESS
4.4 查看匹配調整后結果
為方便查找和分析數據，可借助cbind函數，將AB.IPD、AB.IPD_EM、wt、wt.rs四個數據框按列合并成一個新的數據框combine_data_AB，具體代碼如下：
combine_data_AB<-cbind(AB.IPD,AB.IPD_EM,wt,wt.rs) #數據框按列合并
names(combine_data_AB)[6:7]<-c("sex_centred","txtime_centred") #對變量重新命名
然后通過summarise函數計算匹配調整后的試驗AB中男性的加權比例和造血干細胞移植加權平均時間。
combine_data_AB %>%
summarise(male.prop_AB=weighted.mean(sex,wt),
txtime.mean_AB=weighted.mean(txtime,wt))
將試驗AB調整前與調整后結果分別和試驗AC中基線情況相比較，比較結果如表3所示，在以試驗AC人群特征為基礎調整后的試驗AB結果與試驗AC中人群特征基本一致。
表3
調整前后試驗AB和試驗AC中人群基本特征
4.5 計算HR值
首先通過coxph函數進行Cox比例風險模型分析，并基于參數weights的設置，對匹配調整前后的試驗AB數據分別構建模型，模型結果可通過summary函數呈現。具體代碼如下：
#構建匹配調整后的Cox比例風險模型
cox_maic <- coxph(Surv(futime, death) ~ trt, #構建模型關系
data=combine_data_AB, #指定數據集
weights=wt) #加權調整
summary(cox_maic, conf.int=0.95) #給定置信水平下結果
#構建未經匹配調整的Cox比例風險模型
cox_unmaic<-coxph(Surv(futime, death) ~ trt,
data=combine_data_AB)
summary(cox_unmaic, conf.int=0.95)
為提高分析準確性，將HR值及對應95%上下限取對數，然后根據Bucher間接比較方法，計算得到干預措施B相對于干預措施C的HRBC值及對應95%上下限。因篇幅有限，本研究僅呈現匹配調整后的HR值和間接比較結果，如果想獲得未調整的結果，可將下列代碼中調整后參數替換為未調整參數。具體代碼如下：
logHR_AC<-log(AC.AgD$HR_AC) #對HRAC取對數
logHR_se_AC<-(log(AC.AgD$uHR_AC)-log(AC.AgD$lHR_AC))/(2*1.96) #計算log(HRAC)標準誤
logHR_AB_maic<-log(summary(cox_maic)$conf.int[1]) #對HRAB取對數
logHR_se_AB_maic<-(log(summary(cox_maic)$conf.int[4])-log(summary(cox_maic)$conf.int[3]))/(2*1.96) #計算log(HRAB)標準誤
HR_BC_maic<-exp(logHR_AC-logHR_AB_maic) #計算HRBC
HR_var_BC_maic<-logHR_se_AC^2+logHR_se_AB_maic^2 #計算HRBC方差
lHR_BC_maic<-exp(log(HR_BC_maic)-1.96*sqrt(HR_var_BC_maic)) #計算HRBC95%下限
uHR_BC_maic<-exp(log(HR_BC_maic)+1.96*sqrt(HR_var_BC_maic)) #計算HRBC95%上限
4.6 Bootstrapping法計算HR值
在實際情況中，通常無法考慮全部影響人群之間平衡的效應修飾因子和預后變量，往往只會對經統計方法檢驗顯著的影響變量進行調整；且不同觀察值可能被賦予不同權重，進而導致結果可能受到觀察值之間相關性的影響，直接使用絕對結果將產生偏倚。因此，需要考慮這些可能存在的偏倚，常用的方法有魯棒方差估計（robust variance estimation，RVE）和Bootstrapping法[9]。其中，Austin等[9]研究發現，使用Bootstrapping法可以通過重復抽樣較為正確地估計標準誤和置信區間。因此本研究使用Bootstrapping法來計算標準誤和置信區間，并通過boot函數在R語言中實現這一操作，具體步驟如下。
首先，通過function函數構建boot函數所需的統計量，具體代碼如下：
hr.fun <- function(data, indices) { #確定函數輸入參數
d <- data[indices,] #根據輸入的索引構建子數據集
fit <- coxph(Surv(futime, death) ~ trt ,weights=wt, data = d) #擬合Cox比例風險模型
ci <- confint(fit) #計算置信區間
c(exp(coef(fit)), exp(ci))} #將模型估計的HR值和置信區間以向量形式返回
然后通過boot函數完成Bootstrapping估計，具體代碼如下：
set.seed(123) #設定種子，確保結果可重現
HR_AB_bootstraps <- boot(data = combine_data_AB, #指定數據
statistic = hr.fun, #指定統計量
R=1000) #指定抽樣次數
然后從Bootstrapping結果中將HRAB值及95%上下限選擇出來并取對數和計算標準誤。具體代碼如下：
logHR_AB_boot <- log(median(HR_AB_bootstraps$t)) #對重復抽樣后中位HRAB取對數
boot_HR_ci_AB_HR <- boot.ci(boot.out = HR_AB_bootstraps, type="perc") #使用百分位數法計算置信區間
lHR_AB_boot<-boot_HR_ci_AB_HR$percent[4] #選擇百分位數法計算置信區間下限
uHR_AB_boot<-boot_HR_ci_AB_HR$percent[5] #選擇百分位數法計算置信區間上限
logHR_se_AB_boot<-(log(lHR_AB_boot)-log(uHR_AB_boot))/(2*1.96) #計算重復抽樣后中位HRAB標準誤
然后同樣根據Bucher法，計算得到更為精確的干預措施B相對于干預措施C的HRBC值及對應95%上下限（代碼類同本文“4.5 計算HR值”）。
最后，用data.frame函數，將未調整、匹配調整后以及Bootstrapping法精確估計計算間接比較結果匯總在一張表，結果如表4所示，具體代碼如下：
表4
不同方法下HRBC及95%置信區間
HR_BC<-c(HR_BC_unmaic,HR_BC_maic,HR_BC_boot)
lHR_BC<-c(lHR_BC_unmaic,lHR_BC_maic,lHR_BC_boot)
uHR_BC<-c(uHR_BC_unmaic,uHR_BC_maic,uHR_BC_boot)
method<-c("unmaic", "maic","boot")
outcome<-data.frame(method,HR_BC,lHR_BC,uHR_BC) #結果以數據框呈現
outcome
從表4結果中可以看出，與干預措施C相比，干預措施B降低了死亡風險，并且匹配調整改變了效應大小，符合我們對潛在效應修飾因子或預后變量（男性比例、造血干細胞平均移植時間）調整后預期的評估；同時，經過Bootstrapping法計算得到的HR值與經基礎加權得到的結果基本一致。
以上就是在R語言中實現基于錨定MAIC法的具體操作過程。如果缺少共同對照，則可以選擇使用非錨定MAIC，其匹配調整的過程與錨定的MAIC過程基本相似，只需將本研究中試驗AB和試驗AC的數據分別當做干預措施B（IPD資料）和干預措施C（AgD資料）的數據，從而獲得干預措施C人群結構下匹配調整后的干預措施B加權數據。由于缺乏共同對照，非錨定MAIC不能直接計算加權調整后干預措施B與干預措施C的相對效應，但可通過數字化工具提取干預措施C生存曲線數據，并重構干預措施C的IPD，然后與匹配調整后的干預措施B數據一同結合在Cox比例風險模型中，進而求得干預措施B與干預措施C的相對效應。有關重構IPD及生存分析的方法已有文獻詳細介紹了在R語言中實現的方法，可參考石豐豪等[10]的研究。
5 討論
本研究對MAIC原理和方法進行介紹，并結合survival內置的myeloid生存數據和虛擬AgD資料，展示了如何在R語言中實現錨定MAIC過程。結果顯示，匹配調整后改變了間接比較結果，這也證明了效應修飾因子或預后變量對結局指標的影響。
隨著依托單臂試驗附條件上市的創新藥物不斷涌現，以及交叉試驗的存在，如何有效評估不同臨床試驗之間的效應差異是一個需要重視的問題。相較于傳統的間接比較，在依托IPD資料的情況下使用MAIC，能夠充分利用已有數據，通過匹配調整有效地平衡好RCT間效應修飾因子的影響，降低比較結果的偏倚風險。鑒于MAIC的適應性，其應用也越來越廣泛。例如，在生存數據中，Chen等[11]在鱗狀非小細胞肺癌一線治療研究中，Maloney等[11]在復發或難治性大B細胞淋巴瘤研究中分別使用了基于錨定MAIC和非錨定MAIC法；在其他類型數據中（二分類或連續性數據），Mt-Isa等[13]采用錨定MAIC法比較不同肺炎球菌疫苗的免疫應答情況，Thom等[14]使用非錨定MAIC法比較克立硼羅軟膏與外用鈣調磷酸酶抑制劑治療輕中度特應性皮炎患者的療效。
MAIC的應用也有一定的局限性[2,7,15]：① 只能對比兩組情況，不適用于復雜證據體；② 默認將提供匯總資料的人群視作決策目標人群，其外推性有限，無法在指定的目標人群中進行效果預測；③ 調整多種基線因素的能力取決于擁有IPD的試驗中是否有足夠數量的患者，且匹配過程中可能存在數據丟失的問題，需要較大樣本量；④ 與傳統傾向得分匹配加權不同，MAIC中AgD資料的可用性阻礙了使用現有方法來檢查傾向評分模型的擬合和校準情況；⑤ 對生存數據進行分析時，更適合應用于等比例假設條件下。
此外，與MAIC法相類似的人群匹配調整方法還有模擬治療比較（simulated treatment comparison，STC）和多水平網狀Meta回歸（multi-level network meta-regression，ML-NMR）。STC與MAIC在適應范圍上無優劣之分，但STC法處理非線性模型的結果（例如Logit回歸或事件時間模型）時，可能會增加偏倚[16]，Veroniki等[17]研究也發現MAIC法應用的比例高于STC法；ML-NMR法雖然相較于MAIC法能夠實現包含多個研究的比較，但是其應用仍只局限于二分類結局數據，此外，能否將單臂研究整合到證據體中也需進一步探索驗證[2]。
綜上所述，在缺乏直接比較證據的前提下，MAIC是一種有效的替代手段，可以為相關研究和決策提供一定的證據支持。此外，Gregory等[18]已經開發出在R語言中進行MAIC的集成程序包MAIC，但該程序包中函數及處理過程多是集成式，在一定程度上簡化了語句，但在理解上可能會存在一定的難度，尤其是對于MAIC處理過程和相關原理的理解方面，讀者可在理解本研究的基礎上進一步探索學習。
聲明 本研究不存在任何利益沖突。
