随着单臂临床试验日益增多，缺乏头对头临床研究，未调整的间接比较和网状Meta分析应用局限性日益突出。当个体数据数量有限时，匹配调整间接比较能够充分地利用单项研究个体数据并通过倾向得分匹配其他研究的汇总数据，有效地平衡不同试验中人群差异带来的潜在偏倚，并完成目标干预措施间的疗效比较。本研究一方面围绕匹配调整间接比较有关概念和原理进行介绍，另一方面，基于当前肿瘤药物评价的广泛应用，系统展示如何基于R语言使用锚定匹配调整间接比较法针对生存数据进行分析，以期为科学循证决策提供参考。
引用本文： 胡宏飞, 马越, 裘怡津, 李宇昕, 周挺. 匹配调整间接比较用于生存数据分析：在R语言中的实现. 中国循证医学杂志, 2024, 24(1): 105-111. doi: 10.7507/1672-2531.202304073 复制
在特定人群中，头对头随机对照试验（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等[12]在复发或难治性大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等[12]在复发或难治性大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处理过程和相关原理的理解方面，读者可在理解本研究的基础上进一步探索学习。
声明 本研究不存在任何利益冲突。
