卫生技术评估Markov Model的敏感性分析

Author

沈迟

Published

October 23, 2025

1. 演示数据

数据来源：屈水令,潘晓平,王爱玲,等.Markov模型的R实现及在卫生经济学评价中的应用[J].中国卫生统计,2018,35(05):670-672+676.DOI:CNKI:SUN:ZGWT.0.2018-05-008.

2 Markov Model主分析

2.1 分析平台信息

print(version)

               _                           
platform       x86_64-apple-darwin17.0     
arch           x86_64                      
os             darwin17.0                  
system         x86_64, darwin17.0          
status                                     
major          4                           
minor          2.3                         
year           2023                        
month          03                          
day            15                          
svn rev        83980                       
language       R                           
version.string R version 4.2.3 (2023-03-15)
nickname       Shortstop Beagle

2.2 加载分析所需包

library(heemod)
library(ggplot2)
library(cowplot)

options(scipen = 100)

2.3 设置Markov Model参数

模型参数的设定有两种方式：一种是通过define_parameters函数统一指定模型各个环节需要的参数；另一种是在后续功能函数（如define_transition函数）中单独设置。对于参数较多时，建议对于重要参数选择前者，优势在于当某些参数需要变更时，可以极大减少工作量。对于非重要参数可以选择后者。两种方式可以同时使用。

2.3.1 健康状态及状态间转移概率

提前设置参数主要分三类：

不同健康状态的转移概率，每一种策略都需要设置。
成本和效用，同样是每一种策略都需要设置。
贴现率，这个也可以不设置。

param <- define_parameters(
  discount_rate = 0.3,

  p_H2H_trt = 0.89,
  p_H2A_trt = 0.06,
  p_H2D_trt = 0.05,

  p_A2H_trt = 0,
  p_A2A_trt = 0.93,
  p_A2D_trt = 0.07,

  p_D2H_trt = 0,
  p_D2A_trt = 0,
  p_D2D_trt = 1,
  
  #---#
  p_H2H_ctl = 0.55,
  p_H2A_ctl = 0.25,
  p_H2D_ctl = 0.20,

  p_A2H_ctl = 0,
  p_A2A_ctl = 0.40,
  p_A2D_ctl = 0.60,

  p_D2H_ctl = 0,
  p_D2A_ctl = 0,
  p_D2D_ctl = 1,
  
  #---#
  cost_H = 16000,
  cost_A = 36000,
  cost_D = 0,

  u_H = 0.88,
  u_A = 0.7,
  u_D = 0
)

2.3.2 不同策略的转移矩阵

mat_trans_trt <- define_transition(
  state_names = c("HIV", "AIDS", "Death"),
  p_H2H_trt, p_H2A_trt, p_H2D_trt,
  p_A2H_trt, p_A2A_trt, p_A2D_trt, 
  p_D2H_trt, p_D2A_trt, p_D2D_trt)

mat_trans_ctl <- define_transition(
  state_names = c("HIV", "AIDS", "Death"),
  p_H2H_ctl, p_H2A_ctl, p_H2D_ctl,
  p_A2H_ctl, p_A2A_ctl, p_A2D_ctl, 
  p_D2H_ctl, p_D2A_ctl, p_D2D_ctl)

mat_trans_trt

A transition matrix, 3 states.

      HIV       AIDS      Death    
HIV   p_H2H_trt p_H2A_trt p_H2D_trt
AIDS  p_A2H_trt p_A2A_trt p_A2D_trt
Death p_D2H_trt p_D2A_trt p_D2D_trt

2.3.3 不同策略的参数组合

在定义define_strategy之前，也可以使用define_state函数分别定义不同健康状态的成本和效果、效用等。这种方法在不同的策略的初始成本和效用效果等均一样时比较推荐，可以减少代码输入量。

strat_trt <- define_strategy(
    transition = mat_trans_trt,
    HIV = define_state(
        utility = u_H,
        cost = discount(cost_H, discount_rate)
    ),
    AIDS = define_state(
        utility = u_A,
        cost = discount(cost_A, discount_rate)
    ),
    Death = define_state(
        utility = u_D,
        cost = discount(cost_D, discount_rate)
    )
)

strat_ctl <- define_strategy(
    transition = mat_trans_ctl,
    HIV = define_state(
        utility = u_H,
        cost = discount(cost_H - 6000, discount_rate)
    ),
    AIDS = define_state(
        utility = u_A,
        cost = discount(cost_A - 16000, discount_rate)
    ),
    Death = define_state(
        utility = u_D,
        cost = discount(cost_D, discount_rate)
    )
)

2.3.4 运行模型

res_mod <- run_model(
  "Treat" = strat_trt, "Untreat" = strat_ctl, 
  parameters = param,
  init = c(10000, 0, 0), # 指初始HIV状态的样本为10000，其他AIDS和Death样本为0
  cycles = 10,
  cost = cost,
  effect = utility
)

summary(res_mod)

2 strategies run for 10 cycles.

Initial state counts:

HIV = 10000
AIDS = 0
Death = 0

Counting method: 'life-table'.

 

Counting method: 'beginning'.

 

Counting method: 'end'.

Values:

         utility      cost
Treat   64642.96 653758490
Untreat 21549.18 244969683

Efficiency frontier:

Untreat -> Treat

Differences:

      Cost Diff. Effect Diff.     ICER    Ref.
Treat   40878.88     4.309378 9486.028 Untreat

改变plot函数中的type参数可以获得不同的图形：

type = c(“counts”, “ce”, “values”),type = “counts” represents state memberships (corrected) by cycle, type = “ce” plots models on the cost-efficiency plane with the efficiency frontier, and type = “values” state values per cycle.
panels = c(“by_strategy”, “by_state”, “by_value”),
values = NULL,
strategy = NULL,
states = NULL,
free_y = FALSE,
bw = FALSE

plot(res_mod, 
     panel = "by_state",
     type = "counts") +
scale_color_brewer(
        name = "Strategy",
        palette = "Set1") + 
labs(title = "type = counts") +
theme_bw() +
    theme(text = element_text(family = "Times"))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

plot(res_mod, 
     panel = "by_state",
     type = "ce") +
scale_color_brewer(
        name = "Strategy",
        palette = "Set1") + 
labs(title = "type = ce") +
theme_bw() +
    theme(text = element_text(family = "Times"))

plot(res_mod, 
     type = "values") +
scale_color_brewer(
     name = "Strategy",
    palette = "Set1") + 
labs(title = "type = values") +
theme_bw() +
    theme(text = element_text(family = "Times"))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

3 敏感性分析

3.1 确定性敏感分析（Deterministic Sensitivity Analysis）

se <- 
    define_dsa(
        cost_H, 6000, 25000,
        cost_A, 12300, 66000,

        u_H, 0.58, 0.92,
        u_A, 0.47, 0.82,

        discount_rate, 0.1, 0.5)

res_dsa <- run_dsa(res_mod, dsa = se)

Running DSA on strategy 'Treat'...

Running DSA on strategy 'Untreat'...

summary(res_dsa)

A sensitivity analysis on 5 parameters.

Parameters:
  -cost_H
  -cost_A
  -u_H
  -u_A
  -discount_rate

Sensitivity analysis:

                             utility  cost       .n_indiv .par_value_eval
Untreat, cost_A = 12300      21549.18  113836361 10000    12300.00       
Treat, cost_A = 12300        64642.96  531843343 10000    12300.00       
Untreat, cost_A = 66000      21549.18  410961230 10000    66000.00       
Treat, cost_A = 66000        64642.96  808081459 10000    66000.00       
Untreat, cost_H = 25000      21549.18  365847469 10000    25000.00       
Treat, cost_H = 25000        64642.96  917329635 10000    25000.00       
Untreat, cost_H = 6000       21549.18  110661032 10000     6000.00       
Treat, cost_H = 6000         64642.96  360901661 10000     6000.00       
Untreat, discount_rate = 0.1 21549.18  304358899 10000        0.10       
Treat, discount_rate = 0.1   64642.96 1087786298 10000        0.10       
Untreat, discount_rate = 0.5 21549.18  212058626 10000        0.50       
Treat, discount_rate = 0.5   64642.96  481548458 10000        0.50       
Untreat, u_A = 0.47          19435.81  244969683 10000        0.47       
Treat, u_A = 0.47            60497.59  653758490 10000        0.47       
Untreat, u_A = 0.82          22651.81  244969683 10000        0.82       
Treat, u_A = 0.82            66805.77  653758490 10000        0.82       
Untreat, u_H = 0.58          16395.60  244969683 10000        0.58       
Treat, u_H = 0.58            46906.62  653758490 10000        0.58       
Untreat, u_H = 0.92          22236.33  244969683 10000        0.92       
Treat, u_H = 0.92            67007.81  653758490 10000        0.92       
                             Cost     Effect   ICER      Cost Diff.
Untreat, cost_A = 12300          0.00 0.000000 -         -         
Treat, cost_A = 12300        41800.70 4.309378  9699.937 41800.70  
Untreat, cost_A = 66000          0.00 0.000000 -         -         
Treat, cost_A = 66000        39712.02 4.309378  9215.256 39712.02  
Untreat, cost_H = 25000          0.00 0.000000 -         -         
Treat, cost_H = 25000        55148.22 4.309378 12797.256 55148.22  
Untreat, cost_H = 6000           0.00 0.000000 -         -         
Treat, cost_H = 6000         25024.06 4.309378  5806.885 25024.06  
Untreat, discount_rate = 0.1     0.00 0.000000 -         -         
Treat, discount_rate = 0.1   78342.74 4.309378 18179.593 78342.74  
Untreat, discount_rate = 0.5     0.00 0.000000 -         -         
Treat, discount_rate = 0.5   26948.98 4.309378  6253.567 26948.98  
Untreat, u_A = 0.47              0.00 0.000000 -         -         
Treat, u_A = 0.47            40878.88 4.106179  9955.456 40878.88  
Untreat, u_A = 0.82              0.00 0.000000 -         -         
Treat, u_A = 0.82            40878.88 4.415395  9258.261 40878.88  
Untreat, u_H = 0.58              0.00 0.000000 -         -         
Treat, u_H = 0.58            40878.88 3.051101 13398.074 40878.88  
Untreat, u_H = 0.92              0.00 0.000000 -         -         
Treat, u_H = 0.92            40878.88 4.477148  9130.562 40878.88  
                             Effect Diff. Ref.    .nmb      .dnmb    
Untreat, cost_A = 12300      -            -            0.00 -        
Treat, cost_A = 12300        4.309378     Untreat  87480.65  87480.65
Untreat, cost_A = 66000      -            -            0.00 -        
Treat, cost_A = 66000        4.309378     Untreat  89569.32  89569.32
Untreat, cost_H = 25000      -            -            0.00 -        
Treat, cost_H = 25000        4.309378     Untreat  74133.13  74133.13
Untreat, cost_H = 6000       -            -            0.00 -        
Treat, cost_H = 6000         4.309378     Untreat 104257.28 104257.28
Untreat, discount_rate = 0.1 -            -            0.00 -        
Treat, discount_rate = 0.1   4.309378     Untreat  50938.60  50938.60
Untreat, discount_rate = 0.5 -            -            0.00 -        
Treat, discount_rate = 0.5   4.309378     Untreat 102332.36 102332.36
Untreat, u_A = 0.47          -            -            0.00 -        
Treat, u_A = 0.47            4.106179     Untreat  82306.48  82306.48
Untreat, u_A = 0.82          -            -            0.00 -        
Treat, u_A = 0.82            4.415395     Untreat  91582.98  91582.98
Untreat, u_H = 0.58          -            -            0.00 -        
Treat, u_H = 0.58            3.051101     Untreat  50654.16  50654.16
Untreat, u_H = 0.92          -            -            0.00 -        
Treat, u_H = 0.92            4.477148     Untreat  93435.57  93435.57

龙卷风图（tornado diagram）

通过修改result参数，可以在图中获得不同的结果。

result = c(“icer”, “cost”, “effect”)

plot(res_dsa,
     result = "icer",
     strategy = "Treat",
     type = "difference") +
    scale_color_brewer(
        name = "Strategy",
        palette = "Set1"
  ) + 
    theme_bw() +
    theme(text = element_text(family = "Times"))

3.2 概率敏感性分析（Probabilistic Sensitivity Analysis）

rsp <- define_psa(
  discount_rate ~ lognormal(mean = 0.1, sdlog = 0.15),
  
  cost_H ~ gamma(mean = 16000, sd = sqrt(2000)),
  cost_A ~ gamma(mean = 36000, sd = sqrt(3000)),
  
  u_H ~ normal(mean = 0.88, sd = 0.1),
  u_A ~ normal(mean = 0.7, sd = 0.15)
)


pm <- run_psa(
  model = res_mod,
  psa = rsp,
  N = 100
)

Resampling strategy 'Treat'...
Resampling strategy 'Treat'...

Resampling strategy 'Untreat'...
Resampling strategy 'Untreat'...

summary(
  pm, 
  threshold = c(1000, 5000, 6000, 1e4))

2 strategies run for 10 cycles.

Initial state counts:

HIV = 10000
AIDS = 0
Death = 0

Counting method: 'life-table'.

 

Counting method: 'beginning'.

 

Counting method: 'end'.

Values:

         utility       cost
Treat   64912.22 1098486383
Untreat 21568.27  305453369

Net monetary benefit difference:

            1000     5000     6000    10000
Untreat 74968.91 57631.33 53296.93 35959.35
Treat       0.00     0.00     0.00     0.00

Efficiency frontier:

Untreat -> Treat

Differences:

      Cost Diff. Effect Diff.     ICER    Ref.
Treat    79303.3     4.334395 18296.28 Untreat

3.2.1 成本效果平面散点图

绘制增量效果为横坐标，增量成本为纵坐标绘制散点图.

plot(pm, 
    type = "ce") + 
    scale_color_brewer(
        name = "Strategy",
        palette = "Set1"
  ) + 
    theme_bw() +
    theme(text = element_text(family = "Times")) +
    geom_abline(intercept = 0, slope = 25000, color = "blue", linetype = "dashed")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

3.2.2 成本效果可接受曲线

plot(pm, 
     type = "ac", 
     max_wtp = 25000, 
     log_scale = FALSE) + 
    scale_color_brewer(
        name = "Strategy",
        palette = "Set1"
  ) + 
    theme_bw() +
    theme(text = element_text(family = "Times"))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.