R 药物经济学评价：Markov模型构建及markov轨迹图绘制

前言

药物经济学评价中比较常用的模型包括决策树（Decision tree）模型、马尔科夫（Markov）模型、分区生存模型（Partitioned Survival Model,PSM）、微观仿真模拟（Microsimulation）模型、离散事件模拟（Discrete Event Simulation，DES）模型等。其中Markov模型是药物经济学评价中常用的一种建模方法，常用于长期慢性病经济学评估中。该模型是一种特殊的循环决策树模型，是一种将临床事件和相关干预实施的时间因素纳入模型的动态模型，对现实环境中患者健康状态连续变化的一种粗略模拟，是一种离散时点状态转移的模型。

在markov模型中，研究时限被划分为等长的循环周期（Cycles Length），模型中的患者被定义划分为有限个健康状态（Health States），模拟中的每一个患者在每一个循环周期中必须且只能处在其中一个状态。用初始概率（Initial Probabilities）定义模拟开始时一组患者在各种健康状态中的人数分布，并通过转移概率（Transition Probabilities)来定义每一个周期内患者从一种状态转移到另一种状态的可能性。今天的文章以<Modelling the cost effectiveness of lamivudine/zidovudine combination therapy in HIV infection >文献为例，简单介绍Markov模型的构建及markov轨迹图在R语言中的绘制。

模型结构图：

 

转移概率计算

程序包加载及模型状态名称和周期的设定

library(tidyverse)
library(ggsci)
library(gganimate)

states_names <- c("A", "B", "C", "Death") # 状态名称
cycles <- 20                              # 模型周期

通过状态转移人数计算从每个状态到其他状态的转移概率

A.total <- 1251+350+116+17
B.total <- 731+512+15
C.total <- 1312+437

p.A2A <- 1251 / A.total
p.A2B <- 350  / A.total
p.A2C <- 116  / A.total
p.A2D <- 17   / A.total
p.B2B <- 731  / B.total
p.B2C <- 512  / B.total
p.B2D <- 15   / B.total
p.C2C <- 1312 / C.total
p.C2D <- 437  / C.total

转移概率矩阵

通过转移概率（Transition Probabilities)来确定转移概率矩阵

mat_P = matrix(c(p.A2A, p.A2B, p.A2C, p.A2D,
                 0,     p.B2B, p.B2C, p.B2D,
                 0,       0,   p.C2C, p.C2D,
                 0,       0,    0,      1),
               nrow = 4, byrow = T,
               dimnames = list(states_names,
                               states_names))
mat_P

> mat_P
              A         B          C       Death
A     0.7214533 0.2018454 0.06689735 0.009803922
B     0.0000000 0.5810811 0.40699523 0.011923688
C     0.0000000 0.0000000 0.75014294 0.249857061
Death 0.0000000 0.0000000 0.00000000 1.000000000

 检查每一行概率加起来之和是否等于1.

rowSums(mat_P)

> rowSums(mat_P)
    A     B     C Death 
    1     1     1     1 

Markov  trace

使用初始概率（Initial Probabilities）定义模拟开始时一组患者在各种健康状态中的人数分布，然后通过矩阵乘法来计算患者队列人群

initial <- c(1,0,0,0)
M_trace = matrix(rep(0, 4*21),
          ncol=4,
          dimnames = list(0:20,
                          c("A","B","C","Death")))

M_trace[1,] <- initial

for (i in 2:21){
  M_trace[i,] <- M_trace[i-1,] %*% mat_P
}
M_trace

> M_trace
             A           B          C       Death
0  1.000000000 0.000000000 0.00000000 0.000000000
1  0.721453287 0.201845444 0.06689735 0.009803922
2  0.520494846 0.262910628 0.18059602 0.035998510
3  0.375512717 0.257831905 0.27729592 0.089359455
4  0.270914884 0.225616773 0.33806874 0.165399604
5  0.195452434 0.185784574 0.36354831 0.255214678
6  0.141009801 0.147407084 0.36140189 0.350181229
7  0.101731984 0.114117654 0.34053023 0.443620127
8  0.073394875 0.086845747 0.30869729 0.531062087
9  0.052950973 0.065278842 0.27182282 0.609947364
10 0.038201654 0.048620213 0.23401643 0.679161707
11 0.027560709 0.035963116 0.19788955 0.738586622
12 0.019883764 0.026460490 0.16492601 0.788729740
13 0.014345207 0.019389137 0.13581754 0.830448113
14 0.010349397 0.014162175 0.11073351 0.864754914
15 0.007466606 0.010318351 0.08952225 0.892692795
16 0.005386808 0.007502899 0.07185350 0.915256795
17 0.003886330 0.005447095 0.05731440 0.933352173
18 0.002803806 0.003949642 0.04547092 0.947775632
19 0.002022815 0.002860998 0.03590474 0.959211445
20 0.001459366 0.002070768 0.02823342 0.968236444  

Markov  trace  可视化

先把宽格式数据转换成长格式数据，然后使用ggplot（）函数进行markov轨迹图的绘制

M_trace <- as_tibble(M_trace) 
M_trace$Cycles <- 0:20 
t_trace <-
pivot_longer(
  M_trace,
  cols = A:Death,
  names_to = "State",
  values_to = "probability"
)

ggplot(data=t_trace,
       aes(x=Cycles, y=probability, color=State))+
  geom_line()+
  geom_point()+
  ylab("Proportion")+
  theme_bw()+
  scale_color_lancet()+
  theme(axis.ticks = element_blank(),
        panel.grid=element_blank())  

使用gganimate包来展示markov动态过程

ggplot(data=t_trace,
       aes(x=Cycles, y=probability, color=State))+
  geom_line()+
  ylab("Proportion")+
  theme_bw()+
  scale_color_lancet()+
  theme(axis.ticks = element_blank(),
        panel.grid=element_blank())+
  transition_reveal(along = Cycles)+
  labs(title = "Cycles={frame_along}")

通过geom_area()绘制另一种常见的患者分布面积图，每个状态的面积代表了在该状态下患者的生存时间

t_trace$State = factor(t_trace$State, levels=c("Death","C","B","A")) 
ggplot(data = t_trace, aes(x=Cycles, y=probability, fill = State))+
  geom_area()+
  ylab("Proportion")+
  scale_fill_jama()+
  theme_bw()+
  theme(axis.ticks = element_blank(),
        panel.grid=element_blank())

患者总生存及生命年

用1减去每周期的死亡率便得到每周期患者的总生存率.

M_trace$OS <- 1 -M_trace$Death 
ggplot(M_trace, aes(x = Cycles, y = OS)) +
  geom_line(color="steelblue") +
  ylab("Proportion alive") +
  theme_classic() 

患者生命年计算：

sum(M_trace$OS )

8.991207

参考文献

[1] 刘国恩. 中国药物经济学评价指南2020[M]. 北京：中国市场出版社，2020.

[2] Chancellor JV, Hill AM, Sabin CA, Simpson KN, Youle M. Modelling the cost effectiveness of lamivudine/zidovudine combination therapy in HIV infection. Pharmacoeconomics. 1997;12(1):54-66.