使用heemod包进行基于马尔可夫模型的药物经济学确定性敏感性分析

原创 TwT的学习笔记 TwT的R语言学习笔记 2024年12月20日 04:10

前面两期介绍了如何使用heemod包构建马尔可夫模型并进行药物经济学成本-效益基线分析：

使用heemod包构建马尔可夫模型并进行药物经济学分析
TwT的学习笔记，公众号：TwT的R语言学习笔记使用heemod包构建马尔可夫模型并进行药物经济学分析（一）

使用heemod包构建马尔可夫模型并进行药物经济学分析
TwT的学习笔记，公众号：TwT的R语言学习笔记使用heemod包构建马尔可夫模型并进行药物经济学分析（二）

但在实际生活中，输入模型的参数（如成本、效用值等）并非一个确定的数值，它们往往会在一个区间内波动；此外，各个模型参数对于结果的影响力往往也不尽相同的。

敏感性分析（Sensitivity analysis）则可以在一定程度上反映模型的稳定性和各参数的影响力大小。敏感性分析根据分析过程中参数的取值情况又可以分为确定性敏感性分析（Deterministic sensitivity analysis）以及概率敏感性分析（Probabilistic sensitivity analysis）。在确定性敏感性分析中，各参数的取值为一个确定的值，而在概率敏感性分析中，各参数的取值并非一个确定的值，而是服从某种概率分布。

今天介绍如何使用heemod包进行基于马尔可夫模型的药物经济学确定性敏感性分析。我们将继续以之前一期中构建的马尔可夫模型为例。

1）首先，我们回顾一下模型参数的组成以及基线结果：

mod_par #查看模型参数，mod_par通过define_parameters()定义，详见第一期

这里显示我们有32个参数，其中有14个基础参数（rAE_A, rAE_B, cAE_A, cAE_B, cA, cB, cC, cDeath, uAE_A, uAE_B, uPFS, uPD, uDeath，Discount_rate），其他参数是根据这14个参数根据一定运算得到的。由于死亡状态的效用值（uDeath）应当始终为0，我们不将其纳入我们的敏感性分析。后续我们会使用剩余的13个参数进行确定性敏感性分析。

res_mod #查看基线结果，res_mod通过run_model()计算得来，详见第二期

2）使用define_dsa函数定义确定性敏感性分析中各参数的上下限：

def_dsa <- define_dsa( #将参数输入到def_dsa这个对象中    ###接下来根据 参数名称，参数下限，参数上限，的顺序输入    rAE_A, 0.005, 0.02,    rAE_B, 0.04, 0.16,    cAE_A, 400, 600,    cAE_B, 50, 400,    cA, 600, 1100,    cB, ifelse(state_time <= 5, 2000, 0), ifelse(state_time <= 5, 4000, 0),    cC, 30, 70,    cDeath, ifelse(state_time == 1, 14000, 0), ifelse(state_time == 1, 16000, 0),    uAE_A, -0.01, -0.04,    uAE_B, -0.005, -0.02,    uPFS, 0.8, 0.9,    uPD, 0.65, 0.8,    Discount_rate, 0, 0.08)

各参数的上下限需要根据实际情况确定，最好的方法是通过已发表的文献或官方渠道获得，这些数据的来源和上下限的定义方法都需要在论文中提及。如果无法从外部渠道获得某些参数的上下限，通常的定义方法是：成本及其他参数±20%-25%，效用值±10%，但需要注意，无论成本还是效用值都不应出现负数（不良反应的效用值则相反，不应为正）。具体定义方法请参考官方指南（很重要）。

3）使用run_dsa函数进行确定性敏感性分析：

res_dsa <- run_dsa(res_mod, dsa = def_dsa) #将结果保存到res_dsa对象中

4）使用plot函数绘制龙卷风图：

plot(res_dsa,     type = "difference",     result = "icer" #以icer为结果指标     )

在上图中，横轴为A药：B药的增量成本-效果比（Incremental cost-effect ratio, ICER），纵轴为各参数，每行两侧标出的两个值分别是该参数的下限和上限，由于一些参数用到了ifelse函数，这里并没有直接显示数值而是显示了我们的定义方法，显得很不美观。

接下来我们尝试使用ggplot2制作龙卷风图：

1）下载安装并载入ggplot2：

install.packages("ggplot2")library(ggplot2)

2）查看确定性敏感分析的结果：

res_dsa$dsa #结果保存在res_dsa下的dsa对象中

我们可以看到，数据每4行为1组，从左向右每一列分别为某参数当前取值下的成本（cost）、质量调整生命年（qaly）、样本量（.n_indiv）、治疗策略的名称（.strategy_names）、参数名称（.par_names）、参数值（.par_value）、参数的计算值（.par_value_eval）、成本（.cost）和质量调整生命年（.qaly）。

3）可见，其中并没有ICER值，为了作图，我们需要计算每个参数同一取值下的ICER值。根据ICER的计算公式：

可以进行如下操作：

library(dplyr)
icer_results <- res_dsa$dsa %>%  group_by(.par_names, .par_value, .par_value_eval) %>%  summarize(    ICER = (sum(.cost[.strategy_names == "A"]) - sum(.cost[.strategy_names == "B"])) /            (sum(.effect[.strategy_names == "A"]) - sum(.effect[.strategy_names == "B"]))  ) %>%  ungroup()

查看ICER结果：

icer_results

现在，每个参数的上限和下限分别对应一个ICER值，其中.par_names为该参数的名称，.par_value是我们定义的参数的上下限，.par_value_eval为.par_value的计算值，ICER为该参数上限或下限对应的ICER。

4）作图前，我们还需要根据参数对ICER的影响力大小排序：

icer_results <- icer_results %>%  group_by(.par_names) %>%  mutate(icer_diff = max(ICER) - min(ICER)) %>%  ungroup() %>%  mutate(order_rank = dense_rank(desc(icer_diff))) %>%  arrange(icer_diff, .par_names) %>% #对ICER影响力最大的参数排在第一  mutate(.par_names = factor(.par_names, levels = unique(.par_names)))

查看排序结果：

icer_results

现在oder_rank一列代表了各参数的影响力排序。

5）由于ggplot并没有提供制作龙卷风图（Tornado plot）的功能，我们需要自行通过添加geom_rect以达到绘图的目的。geom_rect需要4个参数（xmin, xmax, ymin, ymax）来定义矩形的四个角，所以接下来我们需要根据已有数据来计算这四个参数。

在龙卷风图中，各矩形（如上图中cA对应的蓝色矩形）是由基线ICER横向扩展至该矩形所对应的参数取上限/下限时（如cA取1100时）得到的ICER值，而矩形的纵向高度没有实际意义，可自由调整。因此，首先我们需要得到基线的ICER：

ICER_base <-    (res_dsa$model$run_model[      res_dsa$model$run_model$.strategy_names == "A", #A药的Cost       ".cost"] - res_dsa$model$run_model[        res_dsa$model$run_model$.strategy_names == "B", #B药的Cost        ".cost"]) /    (res_dsa$model$run_model[      res_dsa$model$run_model$.strategy_names == "A", #A药的QALY      ".effect"] - res_dsa$model$run_model[        res_dsa$model$run_model$.strategy_names == "B", #A药的QALY        ".effect"])

查看结果：

ICER_base

与之前取得的结果一致，确认无误。

接下来，我们对每个ICER值计算xmin、xmax、ymin和ymax四个参数，并标记该ICER对应的是该参数的上限还是下限（从上表可见，每个参数对应两个ICER值）：

icer_plot_data <- icer_results %>%    mutate(xmin = pmin(ICER, ICER_base), #xmin取ICER和ICER_base中较小的那个      xmax = pmax(ICER, ICER_base), #xmax取ICER和ICER_base中较大的那个      ymin = as.numeric(.par_names) - 0.8 / 2, #ymin取该行的坐标减去一半的单位宽度      ymax = as.numeric(.par_names) + 0.8 / 2) %>% #ymax取该行的坐标减去一半的单位宽度    group_by(.par_names) %>%    mutate(type = ifelse(.par_value_eval == max(.par_value_eval), "Upper Limit", "Lower Limit")) %>%    ungroup()

查看一下：

icer_plot_data

6）最后，开始作图之前，我们需要把参数名称换成通俗易懂的名称（如将cA换成cost of A）：

transl <- c("cA" = "Cost of A",            "cB" = "Cost of B",            "cC" = "Cost of C",            "rAE_B" = "AE rate with A",            "rAE_B" = "AE rate with B",            "cAE_A" = "Cost on AEs with A",            "cAE_B" = "Cost on AEs with B",            "cDeath" = "1-time cost on death",            "uPFS" = "PFS utility",            "uPD" = "PD utility",            "Discount_rate" = "Discount rate")            icer_plot_data$.par_names <- recode(icer_plot_data$.par_names, !!!transl)

7）开始作图：

ggplot(data = icer_plot_data,       aes(y = .par_names)) +  geom_rect(aes(    ymax = ymax,    ymin = ymin,    xmax = xmax,    xmin = xmin,    fill = type  )) +  theme_bw() +  theme(    axis.title.y = element_blank(),    legend.position = c(0.85, 0.15),    legend.title = element_blank(),    text = element_text(size = 10)  ) +  geom_vline(xintercept = ICER_base) +  scale_fill_manual(values = c("#008080",  "#FF7F50"))

方便起见，我们将以上处理过程打包装进一个DSA_PLOT函数当中：

DSA_PLOT <- function(comparison = c(), #输入两个策略的名称                     res_dsa = NULL, #输入dsa模型                     transl = NULL,  #定义更改参数名的规则                     colors = c("#008080",  "#FF7F50"), #定义龙卷风图的颜色                     width = 0.8, #定义龙卷风图中每个矩形的宽                     legend.position = c(0.85, 0.15)) { #定义图例的位置
  icer_results <- res_dsa$dsa %>%    group_by(.par_names, .par_value, .par_value_eval) %>%    summarize(      ICER = (sum(.cost[.strategy_names == comparison[1]]) - sum(.cost[.strategy_names == comparison[2]])) /         (sum(.effect[.strategy_names == comparison[1]]) - sum(.effect[.strategy_names == comparison[2]]))    ) %>%    ungroup()
  icer_results <- icer_results %>%    group_by(.par_names) %>%    mutate(icer_diff = max(ICER) - min(ICER)) %>%    ungroup() %>%    mutate(order_rank = dense_rank(desc(icer_diff))) %>%    arrange(icer_diff, .par_names) %>%    mutate(.par_names = factor(.par_names, levels = unique(.par_names)))
  ICER_base <-    (res_dsa$model$run_model[      res_dsa$model$run_model$.strategy_names == "A", #A药的Cost       ".cost"] - res_dsa$model$run_model[        res_dsa$model$run_model$.strategy_names == "B", #B药的Cost        ".cost"]) /    (res_dsa$model$run_model[      res_dsa$model$run_model$.strategy_names == "A", #A药的QALY      ".effect"] - res_dsa$model$run_model[        res_dsa$model$run_model$.strategy_names == "B", #A药的QALY        ".effect"])
  icer_plot_data <- icer_results %>%    mutate(xmin = pmin(ICER, ICER_base), #xmin取ICER和ICER_base中较小的那个           xmax = pmax(ICER, ICER_base), #xmax取ICER和ICER_base中较大的那个           ymin = as.numeric(.par_names) - width / 2, #ymin取该行的坐标减去一半的单位宽度           ymax = as.numeric(.par_names) + width / 2) %>% #ymax取该行的坐标减去一半的单位宽度    group_by(.par_names) %>%    mutate(type = ifelse(.par_value_eval == max(.par_value_eval), "Upper Limit", "Lower Limit")) %>%    ungroup()
  icer_plot_data$.par_names <- recode(icer_plot_data$.par_names, !!!transl)
  plot <- ggplot(data = icer_plot_data,         aes(y = .par_names)) +    geom_rect(aes(      ymax = ymax,      ymin = ymin,      xmax = xmax,      xmin = xmin,      fill = type    )) +    theme_bw() +    theme(      axis.title.y = element_blank(),      legend.position = legend.position,      legend.title = element_blank(),      text = element_text(size = 10)    ) +    geom_vline(xintercept = ICER_base) +    scale_fill_manual(values = colors)
  return(plot)}

测试一下这个函数：

transl <- c("cA" = "Cost of A",            "cB" = "Cost of B",            "cC" = "Cost of C",            "rAE_B" = "AE rate with A",            "rAE_B" = "AE rate with B",            "cAE_A" = "Cost on AEs with A",            "cAE_B" = "Cost on AEs with B",            "cDeath" = "1-time cost on death",            "uPFS" = "PFS utility",            "uPD" = "PD utility",            "Discount_rate" = "Discount rate")            plot <- DSA_PLOT(comparison = c("A", "B"),         res_dsa = res_dsa, transl = transl)

plot

没有问题。最后记得保存图形：

ggsave(plot = plot,  filename = "/Path/to/your/csv.csv") #将/Path/to/your/csv.csv替换为你的文件路径

以上只是提供了一个简单的例子，实际操作往往要复杂得多，随着模型模型结构的复杂化以及参数的增加，运行确定性敏感性分析的时间会大幅度地增加。为了解决这个问题，heemod包提供了多线程运算的功能，可以大幅度缩减确定性敏感性分析的运行时间。

以下为使用多线程的方法：

use_cluster(num_cores = 4) #这时填入线程数量，可以根据CPU核心的数量自行调整，但4个线程以上对速度的提升不大res_dsa <- run_dsa(res_mod, dsa = def_dsa)close_cluster()

我们可以比较一下单线程和多线程的运行速度

###单线程start <- Sys.time()res_dsa <- run_dsa(res_mod, dsa = def_dsa)end <- Sys.time()end-start

###多线程start <- Sys.time()use_cluster(num_cores = 4) #这时填入线程数量，可以根据CPU核心的数量自行调整，但4个线程以上对速度的提升不大res_dsa <- run_dsa(res_mod, dsa = def_dsa)close_cluster()end <- Sys.time()end-start

明显快了许多。

微信扫一扫
关注该公众号

继续滑动看下一个