多状态马尔可夫链、生存分析心脏同种异体移植血管病变（CAV）数据可视化|附数据代码
原文链接：https://tecdat.cn/?p=36216
原文出处：拓端数据部落公众号
临床研究和医疗经济学研究中客户经常关注于评估患者在疾病从一种状态发展到另一种状态时的生存预后。标准生存模型仅直接模拟两种状态：存活和死亡。多状态模型允许直接模拟疾病进程，在这些过程中，患者在随机的时间间隔内处于健康或疾病的各种状态，但除了死亡外，进入或离开状态的时间都是未知的。多状态模型在假设死亡时间是已知的同时，可以轻松地容纳具有区间审查的中间状态，但死亡时间可能是右审查的。
将具有区间审查状态的疾病进程动态建模的概念化为连续时间马尔可夫链是一种自然的方式。
以下图表展示了一个可能的疾病进程模型，其中在任何状态下都有一定的死亡可能性，但除此之外，患者将从健康状态发展到轻度疾病，再到重度疾病，最后可能是死亡。
## Loading required package: shape的确，马尔可夫假设意味着患者在各种状态下所花费的时间是呈指数分布的。然而，随机多状态过程的数学理论非常丰富，可以容纳随时间变化且依赖于状态的更现实的风险率模型，以及其他对马尔可夫假设的放宽。此外，R语言（以及其他语言）中有强大的软件使多状态随机生存模型变得实用。
数据集
所探索的数据集是心脏同种异体移植血管病变（CAV）数据集，该数据集包含了622名心脏移植受者的血管造影检查个人病史。这是一个包含2846行和多个协变量的丰富数据集，包括患者年龄和移植后的时间，这两者都可以用作时间尺度，以及四种状态之间的多次状态转换，并且没有缺失值。中间状态的观察是区间审查的，并已记录在不同的时间间隔内。死亡是“确切”的或右审查的。
以下代码创建了一个新变量，该变量保存了每个观察值的原始状态数据，并以tibble格式显示数据。
状态转换表显示了在不同观察时间点，每对状态被观察到的次数。从表中可以看出，有46次从状态2（轻度CAV）转移到状态1（无CAV）的情况，有4次从状态3（重度CAV）转移到健康状态，以及13次从重度CAV转移到轻度CAV的情况。
遵循van den Hout的方法，并假设这些逆向转换是错误分类的，然后修改状态变量以确保没有倒退的情况。
df1 <- df %>% group_by(PM) %>%
mutate(b_age = min(age),
state = cummax(ste)
)这种转换将使状态转换表符合上面的图示，但状态1代表无CAV而不是健康状态。
建立和运行模型
首先，我们为强度矩阵Q设置初始猜测值，该矩阵决定了连续时间马尔可夫链中状态之间的转移率。
# 强度矩阵Q:
qnames <- c("q12","q14","q23","q24","q34")对于这个模型，从状态1到状态2以及从状态1到状态4的转移取决于时间years、移植时患者的年龄b_age以及捐赠者的年龄dage。其他转移仅取决于dage。因此，我们可以看到可以处理随时间变化的协变量，并允许个体状态转移由不同的协变量驱动。
现在，我们为模型设置剩余参数。
obstype = 1 表示观测值是在任意时间点获取的。它们是面板数据中常见的进程“快照”。
center = FALSE 表示在最大似然估计过程中，协变量不会以其均值为中心。该参数的默认值为TRUE。
deathexact = TRUE 表示最终吸收状态是精确观测到的。这是生存数据的定义假设。这相当于为状态4（我们的吸收状态）设置obstupe = 3。
method = BFGS 指示optim()函数使用1970年由Broyden、Fletcher、Goldfarb和Shanno同时发表的优化方法。（查阅这里）这是一种准牛顿方法，它使用函数值和梯度来构建待优化表面的图像。
control = list(trace=0,REPORT=1) 表示将传递给optim()函数的更多参数。
REPORT 如果control$trace为正，则设置“BFGS”，“L-BFGS-B”和“SANN”方法的报告频率。默认为“BFGS”和“L-BFGS-B”每10次迭代一次，或“SANN”每100次温度变化一次。（注意：SANN是C. J. P. Belisle（1992）在《Rd上的一类模拟退火算法的收敛定理》中提出的模拟退火方法的变体，该文章发表于《应用概率杂志》。）
trace 也传递给optim()函数。trace必须是非负整数。如果为正，则生成关于优化进度的跟踪信息。更高的值可能会产生更多的跟踪信息。
首先，检查模型是否收敛。对于BFGS方法，optim()可能返回的收敛代码为：0表示收敛，1表示已达到最大迭代次数限制，51和52表示警告。
接下来，使用一个由Gentleman et al. (1994)提出的可视化测试来评估模型对数据拟合的好坏程度。
重塑观测患病率
pp <- prev_l %>% ggplot() +
geom_line(aes(time, number, color = type)) +
xlab("时间") +
ylab("") +
ggtitle("")
# 按状态进行分面绘图
pp + facet_wrap(~state)从图中可以看出，状态1到3的观测患病率和预测患病率吻合得相当好。
生存曲线和计算
现在，我们可以直接跳转到主要结果，并查看拟合的生存曲线。
p <- df_l %>% ggplot(aes(time, 1 - survival, group = state)) +
geom_line(aes(color = state)) + 
这些曲线表明，能够预防CAV或至少延缓从轻度CAV发展到重度CAV的治疗可能会延长生存时间。此外，模型的马尔可夫结构允许提取与疾病进展和在每个状态中花费的总时间相关的信息。
该表格表明，患者预期可以避免CAV的平均时间约为7年。在进展到轻度CAV后，患者可以再预期有大约5年的时间。
基于强度矩阵Q的更直接计算，给出了从每个“瞬态”生存状态到“吸收”状态（即死亡）的预期时间。
在协变量提供的信息范围内，还可以生成更个性化的预测。例如，以下是一个60岁无CAV的人，在移植后5年，从20岁的捐献者那里接受心脏，其预期死亡时间。
一个相关的量，平均逗留时间，是指每次访问每个状态时预期的平均持续时间。由于我们假设了一个渐进性疾病模型，其中每个患者只访问每个状态一次，因此估计值应接近在每个状态中花费的总时间。然而，Jackson指出，在渐进性模型中，疾病状态的逗留时间将大于疾病状态的预期停留时间，因为在一个状态中的平均逗留时间是基于进入该状态的条件，而患病状态的预期总时间是对一个个体（可能在患病前死亡）的预测。事实上，这正是我们在状态2和3中看到的情况。
风险（危害）比
model_1 还会产生危害比的估计值，这些估计值显示了每个状态转移强度的估计影响。
以下是危害比的表格：
model_1表格显示，时间协变量years对从状态1到状态2的疾病进展有影响，但对从状态1到状态4的过渡影响较小。
协变量b_age，即移植时患者的基线年龄，对在CAV发病前死亡的影响比对过渡到CAV的影响更大。
协变量dage对从状态1的过渡有较小的影响，但此后似乎没有影响。
危害比是通过计算存储在模型对象中的马尔可夫过程对数转移强度上估计的协变量效应的指数来计算的。
要了解这些是如何工作的，首先查看基线危害比。
model_1$Qmatrices$baseline这些基线危害比是根据模型的强度矩阵Q计算得出的，假设没有协变量。它们也可以通过将协变量设置为零，然后提取函数直接从模型中提取出来。
95%置信限是通过假设对数效应的正态性来计算的。
要获取该模型强度矩阵的一个更具代表性的值，可以将协变量设置为它们的预期均值。
在实际操作中，可以通过在模型评估或预测时设置协变量的平均或代表性值来生成强度矩阵的代表性样本。这样做可以提供在给定协变量分布下的模型行为的平均估计。这对于模拟研究或预测新数据点特别有用，尤其是当没有具体的协变量值可用时。
接下来，我们可能希望检查协变量对危害比的贡献。以dage为例，我们来看它对危害比的影响。
接着，我们关注dage对从状态3转移到状态4的强度矩阵的贡献，这在上面的表格中给出为-0.009439。取这个值的指数，得到了我们在上面打印model_1时得到的危害比表中dage从状态3到状态4的转移危害比。
其余协变量的危害比表格由以下给出：
并且，从上面的表格中我们可以看到，时间协变量years对从状态1转移到状态2和状态4的转移强度有显著的影响，而对从状态1转移到状态3的转移强度没有影响。这是因为，一旦患者进入状态2（CAV），他们就不能再返回到状态1或转移到状态3，只能停留在状态2或转移到状态4（死亡）。因此，从状态1到状态3的转移强度为0。
探索转移概率和转移强度
我们同样可以查看不同时间点的状态转移矩阵，并观察这些概率如何随时间变化。这里我们计算了1年时的转移矩阵。
以及5年时的转移矩阵。
此外，我们还可以检查协变量对转移概率的影响。以下是一个基线年龄为35岁的患者在5年时的转移概率。
此外，我们还计算了年龄为60岁时进行手术的患者在5年时的转移概率。
pmatrsm(model_1, t = 5, covates = list(years = 5, b_age = 60, dage = 20请注意，从CAV状态的转移不受影响。
总结来说：连续时间马尔可夫链为处理多状态生存模型提供了一个自然的框架。
