前言
之前的几篇文章介绍了在搭载临床试验的药物经济学评价中常见的一些统计学问题，如缺失数据、成本和效果的相关性、基线特征不平衡以及数据的偏态分布等。为了有效解决这些方法学问题，介绍了一系列统计方法，包括多重插补（Multiple Imputation, MI）、自举法（Bootstrapping）、似不相关回归（Seemingly Unrelated Regression, SUR）等，并提供了相应的R代码示例。
今天的文章基于 《Ben ÂJ, van Dongen JM, El Alili M, Esser JL, Broulíková HM, Bosmans JE. Conducting Trial-Based Economic Evaluations Using R: A Tutorial. Pharmacoeconomics. 2023;41(11):1403-1413 》一文来展示如何同时解决在在搭载临床试验的药物经济学评价中遇到的四个方法学挑战，主要包括以下三部分内容：
搭载临床试验的药物经济学在缺失数据、成本和效果的相关性、基线特征不平衡以及数据的偏态分布四个统计方法面临的挑战和该教程推荐的解决方法； 通过一个完整的R语言实例（从数据导入到ICER的计算和成本效果可接受曲线的绘制）来展示如何结合bootstrapping法的多重插补（MI）和调整的SUR模型来解决上述的四个统计问题； 讨论该教程未涉及的一些统计问题（多臂试验、贝叶斯方法、线性混合模型等等）
方法学挑战与解决方案
1. 缺失数据（Missing Data）
在搭载临床的药物经济评价中，缺失数据是一个常见的统计问题，尤其是和成本和效果相关的数据，一旦某些数据缺失，可能导致分析样本不完整。简单地排除缺失数据样本或使用均值插补、最后观察值前移等方法无法有效应对不确定性，可能导致偏倚。
解决方法：为了有效处理缺失数据，文章推荐使用多重插补（Multiple Imputation，MI）方法。通过多重插补，研究者能够根据其他观测到的数据预测缺失值，并生成多个完整的数据集，从而避免由于丢失数据而产生的偏差。MI特别适用于缺失数据是“随机缺失”（Missing at Random，MAR）的情况。通过这种方法，可以使得分析结果更加稳健和高效。
2. 成本与效果的相关性（Correlated Costs and Effects）
在搭载临床的药物经济评价中成本和效果通常存在一定的相关性，这种相关性可能会导致传统分析方法失效，例如，效果更好的治疗可能会伴随更高的成本，反之亦然。
解决方法：为了处理成本与效果之间的相关性，文章采用了似不相关回归（SUR）与自举法（Bootstrapping）相结合的方法。SUR模型可以同时分析成本和效果，保持两者之间的相关性，而Bootstrapping则可以有效估计参数的不确定性，确保结果的稳健性。
3. 基线不平衡（Baseline Imbalances）
虽然随机对照试验（RCT）在理论上保证了基线平衡，但实际中有限样本可能仍导致干预组与对照组之间的基线变量（如年龄、性别）存在显著差异。若这些变量与成本或效用相关，未调整基线失衡可能导致结果偏倚。
解决方法：为了解决基线不平衡问题，文章推荐使用回归调整（Regression Adjustment）方法。通过在分析中引入基线特征作为协变量，能够消除基线差异对成本和效果估计的影响，从而提高估计的精度。
4. 成本和效果的偏态分布（Skewness of Costs and/or Effects）
成本和效果数据通常具有偏态分布，尤其是成本数据，通常会有少数患者的成本异常高，这会影响估计的精确度。如果数据存在偏态分布，传统的统计分析方法可能不适用。
解决方法：为了处理数据的偏态问题，文章采用了非参数自助法（Non-Parametric Bootstrapping），该方法不要求数据服从特定的分布假设，可以通过重采样方法获得对成本效益比（ICER）的稳健估计。此外，该方法还能够保持成本和效果之间的相关性，提高分析的有效性。
R语言教程
数据和代码下载地址：https://github.com/angelajben/R-Tutorial
数据集情况
数据集包括 200 名参与者（106 名在对照组，94 名在干预组），涵盖以下内容：
基线变量：
id：患者标识符。 age：患者年龄（以年为单位）。 sex：性别（0 = 女性，1 = 男性）。 E：基线效用值（utility value）。 C：基线成本值（cost value）。 Tr：治疗分配（0 = 对照组，1 = 干预组）。 随访变量：
Em1 至 Em4：3、6、9、12 个月时的效用值。 Cm1 至 Cm4：3、6、9、12 个月时的成本值。
1. 加载相应程序包
使用library函数加载相应程序包，如果没有安装某个R包，可以使用install.packages("包名")进行下载
library(mice) # 用于数据插补
library(systemfit) # 用于似不相关回归（seemingly unrelated regression）
library(car) # 用于运行似不相关回归
library(boot) # 用于自助法（bootstrap）
library(ggplot2) # 用于绘制成本效过平面图和成本效果接受曲线（CEAC）
library(ggpointdensity) # 用于绘制成本效果平面图
library(readxl) # 用于导入 .xlsx 格式的数据
library(tidyverse) # 用于数据操作
2.数据导入
使用read_excel函数从指定路径读取Excel文件，并将其存储到名为dataset的变量中。
dataset <- read_excel("data/dataset.xlsx")
head(dataset)
# A tibble: 6 × 14
id age sex Tr C E Cm1 Em1 Cm2 Em2 Cm3
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 60 0 1 452. 0.562 1783. 0.709 128. 0.987 340.
2 2 44 0 0 226. 0.892 1766. 0.662 630. 0.904 1132.
3 3 46 0 1 6400. 0.356 1352. 0.807 720. 0.791 409.
4 4 49 1 0 420. 0.749 236. 0.624 2889. 0.323 223.
5 5 58 0 0 653. 0.836 624. 0.669 1071. 0.652 4784.
6 6 52 0 0 991. 0.581 306. 0.851 650. 0.652 289.
# ℹ 3 more variables: Em3 <dbl>, Cm4 <dbl>, Em4 <dbl>
3.缺失值的多重插补
使用mice用于多重插补，处理缺失数据。
# 1 按治疗组 (Tr) 拆分数据集
Tr0 <- subset(dataset, Tr == 0) # 提取对照组数据
Tr1 <- subset(dataset, Tr == 1) # 提取干预组数据
# 2. 创建预测矩阵，排除Tr作为预测变量
predMat <- make.predictorMatrix(dataset)
predMat[,'Tr'] <- 0 # 确保`Tr`不作为预测变量参与插补
# 3 按治疗组运行多重插补 (MI) 过程，并将结果合并
imp.Tr0 <- mice(Tr0, m=5, method="pmm", predictorMatrix = predMat, seed = 1234, printFlag = FALSE)
imp.Tr1 <- mice(Tr1, m=5, method="pmm", predictorMatrix = predMat, seed = 1234, printFlag = FALSE)
# 对两个治疗组分别执行多重插补（使用pmm方法），生成5个插补数据集。seed用于设置随机数种子，确保结果可重复。
# 4 合并并储存插补后的数据集
imp <- rbind(imp.Tr0, imp.Tr1) # 合并对照组和干预组的插补数据
impdat <- complete(imp, action = "long", include = FALSE) # 将插补后的数据合并为长格式
#5 提取用于 Rubin 规则的插补次数
M <- imp[["m"]] # 提取多重插补的次数（5次）
#6 计算插补后的成本
impdat$Tcosts <- (impdat$Cm1 + impdat$Cm2 + impdat$Cm3 + impdat$Cm4)
#7 计算插补后的QALY
impdat$QALY <- 1/2 * ((impdat$E+impdat$Em1)*1/4 +(impdat$Em1+impdat$Em2)*1/4+(impdat$Em2+impdat$Em3)*1/4+(impdat$Em3+impdat$Em4)*1/4)
# 8. 将插补后的数据分组存储
impdata <- split(impdat, f = impdat$.imp) # 按插补数据集分组
4. bootstrapping结合似不相关回归模型
systemfit用于拟合似不相关回归（SUR）模型,boot用于bootstrapping分析。
# 9. 定义SUR模型拟合函数
fsur <- function(x, i){
dataset <- x[i,] # 提取bootstrapping样本数据
r1 <- Tcosts ~ Tr + C + age + sex # 成本回归模型
r2 <- QALY ~ Tr + E + age + sex # 效用回归模型
fitsur <- systemfit(list(costreg = r1, effectreg = r2), "SUR", data = dataset) # 拟合SUR模型
betas <- fitsur$coefficients # 提取回归系数
return(c(betas[["costreg_Tr"]], betas[["effectreg_Tr"]])) # 返回增量成本和效用
}
# 10. 对每个插补数据集应用bootstrapping
bootce <- lapply(impdata, function(x) boot(data = x, statistic = fsur, R = 5000))# 对每个插补数据集应用bootstrapping，重复5000次抽样，估计增量成本和增量效用。
5. 提取bootstrapping结果的统计量
从bootstrapping结果中提取增量成本（cost_diff）和增量效用（effect_diff），并将它们合并成一个矩阵。
# 11. 提取插补数据集的统计量
imputed <- lapply(bootce, function(x) ((x[["t0"]]))) # 提取每个插补数据集的统计量
imputed <- lapply(imputed, setNames, c("cost_diff", "effect_diff")) # 给列命名
imputed <- as.matrix(reduce(imputed, bind_rows)) # 合并所有插补数据集的结果
# 12. 提取bootstrapping统计量
postboot <- lapply(bootce, function(x) as.data.frame((x[["t"]]))) # 提取引bootstrapping结果的统计量
postboot <- lapply(postboot, setNames, c("bootcost_diff", "booteffect_diff")) # 给列命名
6.使用鲁宾规则（Rubin's rules）合并成本-效果结果
计算合并后的增量成本、增量效用、ICER及计算置信区间（95% CI）。
# 13. 合并插补数据集的统计量
pooled <- apply(imputed, 2, mean) # 对所有插补数据集的结果取平均
cost_diff_pooled <- pooled[["cost_diff"]] # 合并后的增量成本
effect_diff_pooled <- pooled[["effect_diff"]] # 合并后的增量效用
ICER <- cost_diff_pooled / effect_diff_pooled # 计算增量成本效果比（ICER）
> ICER
[1] 19683.58
# 14. 计算协方差矩阵
cov <- lapply(postboot, function(x) cov(x)) # 对每个bootstrapping样本计算协方差矩阵
cov
$`1`
bootcost_diff booteffect_diff
bootcost_diff 254521.96878 -3.6475702785
booteffect_diff -3.64757 0.0001193538
$`2`
bootcost_diff booteffect_diff
bootcost_diff 223959.506235 -3.5798285202
booteffect_diff -3.579829 0.0001281691
$`3`
bootcost_diff booteffect_diff
bootcost_diff 248717.630435 -3.5619111780
booteffect_diff -3.561911 0.0001225613
$`4`
bootcost_diff booteffect_diff
bootcost_diff 225116.228121 -3.4153979134
booteffect_diff -3.415398 0.0001205374
$`5`
bootcost_diff booteffect_diff
bootcost_diff 249318.490809 -3.7443873655
booteffect_diff -3.744387 0.0001191412
# 15 插补内协方差
W <- 1 / M * (cov[["1"]] + cov[["2"]] + cov[["3"]] + cov[["4"]] + cov[["5"]]) # 计算每个插补数据集内的协方差，按插补数据集数量M求均值
> W
bootcost_diff booteffect_diff
bootcost_diff 240326.764875 -3.5898190511
booteffect_diff -3.589819 0.0001219526
# 16 插补间协方差
B <- matrix(0, ncol = 2, nrow = 2) # 创建一个2x2的零矩阵，用于存储插补数据集之间的协方差
for (i in 1:M) { # 对每个插补数据集进行循环
B <- B + (matrix(imputed[i,], nrow = 2) - pooled) %*% (matrix(imputed[i,], nrow = 1) - pooled) # 计算每个插补数据集之间的协方差
}
B <- 1 / (M - 1) * B # 计算插补数据集之间的协方差矩阵，按插补数M-1标准化
> B
[,1] [,2]
[1,] 1.422230e+04 -1.086904e-02
[2,] -1.086904e-02 4.551635e-06
# 17 汇总协方差
cov_pooled <- (1 + 1 / M) * B + W # 合并每个插补数据集的协方差，得到合并后的协方差矩阵
> cov_pooled
bootcost_diff booteffect_diff
bootcost_diff 257393.521400 -3.6028618979
booteffect_diff -3.602862 0.0001274145
# 18 使用 Rubin 规则估计成本效果的上下区间
Za = 1.95996 # 设定1.96的Z值，对应95%置信区间
LL_cost_pooled <- cost_diff_pooled - (Za * sqrt(cov_pooled[1, 1])) # 计算成本的95%置信区间下限
UL_cost_pooled <- cost_diff_pooled + (Za * sqrt(cov_pooled[1, 1])) # 计算成本的95%置信区间上限
LL_effect_pooled <- effect_diff_pooled - (Za * sqrt(cov_pooled[2, 2])) # 计算效果（QALY）的95%置信区间下限
UL_effect_pooled <- effect_diff_pooled + (Za * sqrt(cov_pooled[2, 2])) # 计算效果（QALY）的95%置信区间上限
> LL_cost_pooled
[1] 299.1533
> UL_cost_pooled
[1] 2287.884
> LL_effect_pooled
[1] 0.04359198
> UL_effect_pooled
[1] 0.08783927
# 19 计算效率损失
FMI = B / (B + W) # 计算每个插补数据集的效率损失（FMI），用来衡量插补的不确定性
LE = FMI / M # 计算整体的效率损失（LE），通过FMI除以插补数据集数量M
> LE
bootcost_diff booteffect_diff
bootcost_diff 0.0111745036 0.0006037201
booteffect_diff 0.0006037201 0.0071960221
7. 绘制增量成本效果图
使用ggplot2绘制增量成本效果平面图，
# 20. 绘制增量成本效果平面图
point <- data.frame(imputed) # 将插补数据转为数据框
boot <- reduce(postboot, bind_rows) # 合并所有bootstraping样本
CE <- ggplot(data = boot, aes(x = booteffect_diff, y = bootcost_diff)) +
geom_pointdensity(aes(booteffect_diff, bootcost_diff), size = 2, alpha = 0.75, show.legend = FALSE, adjust = 0.05) +
geom_point(data = data.frame(x = mean(point$effect_diff), y = mean(point$cost_diff)),
aes(x, y), color = "red", size = 2) + # 绘制红色的增量效用和成本点
labs(x = "Differences in QALY") + # x轴标签
labs(y = "Differences in costs") + # y轴标签
geom_vline(xintercept = 0) + # 绘制垂直参考线
geom_hline(yintercept = 0) + # 绘制水平参考线
scale_x_continuous(breaks = seq(0, 1, by = 0.02)) + # 设置x轴的刻度
scale_y_continuous(breaks = seq(0, 4500, by = 1000)) + # 设置y轴的刻度
theme_minimal() # 设置主题
对增量成本效果散点图的解读和更加丰富的绘图可参考以前的文章。
8. 绘制成本效果可接受曲线
使用ggplot2绘制成本效果可接受曲线（CEAC），显示不同支付意愿下，干预是否具有成本效果的概率。
# 21. 计算增量净效益（INB）
wtp <- seq(0, 80000, 1000) # 设置支付意愿（WTP）范围
INB <- (wtp * effect_diff_pooled) - cost_diff_pooled # 计算增量净效益（INB）
varINB <- wtp^2 * cov_pooled[2, 2] + cov_pooled[1, 1] - 2 * wtp * cov_pooled[1, 2] # 计算INB的方差
seINB <- sqrt(varINB) # 计算标准误
z <- INB / seINB # 计算z值
CEAC <- as.data.frame(wtp) # 转换为数据框
CEAC$prob <- pnorm(z, 0, 1) # 计算成本效益可接受的概率
# 22. 绘制CEAC曲线
CEAC <- ggplot(data = CEAC, aes(x = wtp, y = prob)) +
geom_line(colour = "black", size = 1) + # 绘制CEAC曲线
ylim(0, 1) + # 设置y轴范围
labs(x = "Willingness-to-pay: incremental costs per QALY gained") + # 设置x轴标签
labs(y = "Probability of cost-effectiveness") + # 设置y轴标签
geom_vline(xintercept = 0) + # 绘制垂直参考线
geom_hline(yintercept = 0) + # 绘制水平参考线
theme_minimal() # 设置主题
对成本效果可接受曲线的解读和更加丰富的绘图可参考以前的文章。
讨论
该教程通过逐步的R语言代码展示了如何实施可用的统计方法来同时解决在在搭载临床试验的药物经济学评价中通常遇到的四个方法学挑战。通过结合非参数自举法的多重插补（MI）和调整的SUR模型，以处理缺失数据、相关的成本和效果、基线不平衡以及偏态数据。上面的代码适用于大多数的药物经济学评价，但未涉及以下统计学问题：
线性混合模型（Linear Mixed Model, LMM）
在整群随机试验中（例如，医院或其他医疗单位作为随机化单元）线性混合模型（LMM）可能是一个有效的选择，该方法可以处理数据中的聚集效应（例如，同一医院中的患者可能具有相似的特征），并且能有效地分析成本和效果。
然而LMM方法有一个缺点，即它不能像SUR那样直接处理成本和效果之间的相关性。如果要使用LMM，必须分别分析成本和效果，可能会忽略它们之间的相关性。如果结合自助法，LMM同样可以有效处理相关性问题，该方法对于纵向数据尤其有效，因为它允许在多个时间点上估算成本和效果差异，而不需要事先估算总成本和QALYs。
多臂试验（Multi-arm Trials）
在某些情况下，试验设计可能涉及多个治疗组，而不仅仅是两组比较。例如，治疗组A、B和C之间的比较。在这种情况下，文章的R代码可以进行适应，生成虚拟变量来表示治疗组，并在回归分析中作为协变量进行调整。这样，可以在多个治疗组之间进行多重比较。
虽然文章提供的方法可以扩展到多臂试验，但在多组比较时需要进行多个成对比较，这可能会增加分析的复杂性，需要更加小心地处理数据和解释结果。
汇总估计（Pooling Estimates）
该教程首先使用多重插补（MI）处理缺失数据，然后对插补后的数据进行自举法（Bootstrapping）分析。自举法可以帮助在多重插补的基础上估算不确定性，最终通过Rubin's规则将多个插补数据集的结果进行池化，以得出最终的成本效益比（ICER）及其置信区间。
局限性：有研究指出MI和自举法的结合有一些不确定性，目前尚不明确哪种顺序（先进行MI再自助，还是先进行自助再MI）最优。不同的方法可能会影响估计结果，特别是当样本量较小时，可能会导致不同的结果。
贝叶斯方法（Bayesian Framework）
贝叶斯方法是一种更加灵活的统计方法，可以处理缺失数据、成本与效果的相关性、基线不平衡和偏态数据等问题。贝叶斯方法允许使用先验知识并结合数据来更新估计，尤其适用于复杂的经济评价模型。
局限性：然而贝叶斯方法需要对概率论和统计建模有较深的理解，并且模型通常需要根据特定问题进行定制，这使得它比频率派方法（如SUR和LMM）更加复杂。贝叶斯方法的实现也不像频率派方法那样可以直接在常用统计软件中进行，因此对于大多数研究者来说，学习和实施贝叶斯模型可能具有一定的难度。
处理偏态成本和/或效果数据（Handling of Skewed Costs and/or Effects）
如果成本和效果数据已知遵循某种特定的偏态分布，使用适当的模型（如广义线性模型，GLM）可能会提高效率。然而，正如文章所述，通常很难准确判断成本和效果数据的分布形式，尤其是在样本量较小时，数据可能不完全符合假设的分布。
局限性：该研究选择使用非参数自举法（Bootstrapping），避免了对分布的假设，直接对数据进行重采样以估计不确定性。尽管这种方法可以有效处理偏态数据，但在已知分布的情况下，可能存在一些效率损失。如果研究者更偏向于使用广义线性模型（GLM），也可以根据代码提供的模板进行修改。
参考文献
