综述：多状态过程分析的概述与最新进展

**摘要**

多状态模型为研究疾病过程提供了一个强大的框架，可以用来构建基于强度的以及更具描述性的边际回归模型。它们也是构建疾病过程与动态标记过程联合模型、以及包含随机审查和间歇性观察时间的联合模型的自然基础。本文回顾了如何将多状态模型应用于生命历史数据的方法。还讨论了关于伪值的研究以及如何将随机效应纳入模型以模拟对过程历史的依赖性和过程间的异质性。文中列出了可用于此类分析的软件。

**1 引言**

**1.1 背景**

长期队列研究为研究慢性疾病过程的发作和进展提供了极好的数据来源。当关注某个特定事件发生的时间时，通常会使用生存分析方法来获得洞察。然而，在许多情况下，个体可能会经历多种类型的事件（例如，在心血管研究中，心肌梗塞、出血、非致命性中风和死亡），并且可能希望研究这些事件的共现及其发生时间之间的关系。多状态模型通过分析不同健康状态之间的转换率，为研究这些过程提供了一个灵活的框架。我们的目标是提供多状态分析的一般框架和方法论的综述，讨论替代模型的构建、估计和解释，并展示它们的应用。本文延续了我们第一篇关于基于强度的事件发生时间分析模型的指导性论文[1]。它主要面向熟悉生存分析概念和符号的读者，我们将这些概念和符号推广以适应更复杂的过程。我们强调了多状态似然与生存分析中常规使用的部分似然之间的联系，这有助于使用现代生存分析软件来分析多状态过程。我们希望这一概述能够促进多状态模型在人类健康研究中的广泛和明智的应用，并推进STRATOS计划[2]的总体目标。国际STRATOS（加强观察性研究的分析思维）计划旨在弥合统计方法学的最新进展与应用观察性研究中常用方法之间的差距[2]。这一差距主要是由于缺乏针对观察性数据分析的详细指导，而随机对照试验（RCT）则有广泛采用的CONSORT指南[3]。观察性研究通常涉及比RCT更多样化的目标、设计和数据结构，带来了需要多种统计方法的复杂分析挑战。尽管方法论创新不断涌现，但应用健康研究仍经常依赖于一套有限的传统方法。为了提高观察性研究结论的有效性，STRATOS为具有不同统计水平的研究人员开发了基于证据的指导。最初，STRATOS包括七个主题组（TGs），每个组都针对一个特定的统计分析领域[2]。2015年增加了生存分析组（TG8），提供了关于事件发生时间数据的指导，其首篇贡献发表于2021年，重点是基于强度的失效时间过程建模[1]。本文将这项工作扩展到了多状态设置。涉及不同阶段的疾病可以使用多状态模型自然地进行表征。这些状态可能代表某种进行性疾病的不同的阶段（如肝炎[4]）、间歇性疾病（如慢性支气管炎[5]）中的症状存在与否、癌症临床试验中对治疗的反应的不同阶段[6, 7]，或者COVID-19感染的进程，因为个体在医院中会经历中度、重度和危重疾病状态，然后出院或死亡[8, 9]。在其他情况下，多状态模型已被用于描述昆虫在其生命周期中连续发育阶段的通过[10]、婚姻的形成和解除[11]，或劳动力市场中个体就业状态的变化[12]。仔细建模过程动态可以提供有关疾病自然进程、疾病发作和进展的风险因素以及干预效果的宝贵科学见解。多状态模型还为指定时间依赖的协变量过程和感兴趣的疾病过程的联合模型提供了一个有用的框架；这样的模型可以增进对两个或多个过程之间动态关系的理解。多状态模型的构建始于指定一组代表过程不同状态的状态。这些状态通常对应于不同的健康状态或疾病过程的不同阶段。例如，在退行性疾病中，通常有分级系统来量化损伤程度——修改后的Steinbrocker评分系统基于放射图像来评估风湿病学中的联合损伤程度[13]，而丙型肝炎感染的肝损伤程度则根据纤维化、伴有门脉扩张的纤维化、桥接纤维化和肝硬化等类别进行分级[14]。状态也可以通过离散化连续标记来定义——这是一种在过程受到间歇性观察时常用的策略；例如，Satten和Longini[15]通过建模CD4细胞计数的变化以及艾滋病的诊断来研究HIV中的免疫功能。在其他情况下，状态是根据事件的发生来简单定义的。例如，在癌症试验中，通常区分存活且无复发的状态、复发后的存活状态和死亡状态；死亡状态还可以进一步区分是否在复发前或复发后发生。状态集通常具有有限数量的元素，我们用整数来标记，但也可以对具有可数状态的过程进行建模（例如，复发事件过程[16]）。通常有一部分状态可以从给定状态直接进入，这部分状态是根据上下文确定的。图1包含了一些使用多状态模型可以分析的丰富过程的示例状态空间图。状态空间图中的箭头表示可以直接进行的转换。当过程可以终止时，个体无法退出的吸收状态集用表示。

**1.2 多状态过程的例子**

以下公开可用的数据集在文中被用来说明多状态分析的各种方法。代码可以在支持信息文件中找到。

**1.2.1 结肠癌中的复发和死亡**

在20世纪80年代进行的一项临床试验中，研究了左旋咪唑和氟尿嘧啶作为切除结肠癌辅助疗法的使用[19]，共有929名C期疾病患者被随机分配到三个组之一：观察组、单独使用左旋咪唑组或左旋咪唑与氟尿嘧啶联合使用组。患者招募始于1984年3月，并于1987年10月完成。癌症复发时间和生存时间被认为是关键结果。我们考虑图1d中的疾病-死亡过程，其中状态1表示癌症复发时的状态，状态2表示无复发的死亡，状态3表示复发后的死亡。状态0对应于治疗分配时的初始状态，我们将其设为时间起点（）。该数据集可在R包survival中找到[20]。由于左旋咪唑没有效果，我们将观察组和左旋咪唑组合并作为对照组，并将治疗编码为5FU+Lev (1) 对比 Control (0)。共有468人经历了复发，452人被观察到死亡；大多数死亡（414例）发生在复发之后。我们在第2.2节中重新讨论了这个例子，以说明累积转换强度和随时间变化的状态占用概率的估计。在第2.4节中，我们展示了基于强度的回归模型的发现。

**1.2.2 银屑病关节炎中的关节损伤**

我们考虑的第二个研究涉及银屑病关节炎患者的关节损伤研究，这是一种以关节炎症和损伤以及皮肤受累为特征的自身免疫性疾病。多伦多银屑病关节炎诊所从1977年开始招募患者，并随时间跟踪他们以研究临床环境中的疾病过程[21]。临床和放射学评估定期进行，但不同患者之间的检查时间以及同一患者随时间的变化很大。这里我们考虑了从msm R包[23]中导出的数据集[22]中305名银屑病关节炎患者的关节损伤进展。共有305名患者，平均随访时间为5.5年（最小值=0.1年，最大值=19.2年），我们关注至少有2次就诊记录的患者——平均就诊次数为2.6次（最小值=2次，最大值=7次）。定义了四种损伤状态：轻度、中度、重度和非常严重的损伤。具体来说，状态0对应于0个受损关节，状态1对应于1-4个受损关节，状态2对应于5-9个受损关节，状态3对应于10个或更多受损关节。这是一个进行性过程（即转换是不可逆的），转换发生在连续时间内，但关于受损关节数量的信息仅在每两年进行一次的定期放射学检查时收集。实际检查时间在患者之间差异很大，导致随机访问时间。图4a展示了五名患者的招募时间和随后的放射学检查时间。在第4节中，我们考虑了图1b，其中包含四个状态和两个协变量——基于积液数量和沉降率——并模拟它们对转换到更严重损伤状态的影响。积液是指关节周围组织因液体积聚而发生的肿胀，这表明疾病处于活跃阶段；积液的数量是衡量关节或患者疾病活动程度的指标。红细胞沉降率是一种血液标志物，反映了全身炎症的水平。

1.2.3 鹿特丹肿瘤库数据
鹿特丹肿瘤库的数据包括1546名患有淋巴结阳性乳腺癌的患者，这些患者在1978年至1993年间接受了肿瘤切除手术，这些数据可以在survival R包[20]中找到。我们将肿瘤切除手术的日期作为图1d的起始时间；复发日期、无复发死亡日期和复发后死亡日期分别是进入状态1、2和的状态的时间。预后基线变量包括手术时的年龄、绝经状态、肿瘤大小、肿瘤分级、初次活检中的雌激素和孕激素受体水平、激素治疗以及化疗情况。在1546名患者中，有924人出现了疾病复发（63%），106人在没有复发迹象的情况下死亡（7%），771人在复发后死亡（占出现癌症复发的患者的79%）。这个数据集用于展示第5.1节中讨论的基于脆弱性的方法。

1.3 论文概述
在第2.1节中，我们定义了强度函数，这些函数是多状态过程完整模型的基础构建块。介绍了不同类别的强度函数，它们根据控制风险的时间尺度进行区分。然后我们讨论了马尔可夫过程的转移概率矩阵的计算，并考虑了强度的泛函，这些通常是推断的目标。第2.2节介绍了针对受右删失影响的单样本的非参数估计方法。第2.3节讨论了基于乘法强度的模型的构建，以研究固定或变化协变量对转移率的影响。第2.4节基于时间变化协变量、删失时间和感兴趣过程的联合模型推导了似然函数。第3节重点介绍了使用所谓伪值估计状态占用概率的回归模型。第4节讨论了在连续时间内运行的过程被间歇观察的情况，例如第1.2.2节中的例子。第5节探讨了在多状态疾病过程建模中使用随机效应的方法。第6节提供了可用于多状态分析的统计包和函数的综述。第7节总结了多状态分析的优点和局限性，第8节提出了结论性意见和未来研究的方向。为了便于参考，支持信息文件的第S5节提供了论文中使用的符号表。

2 方法论
2.1 符号和基础
我们用t表示所研究过程的时间起点，并假设除非另有说明，否则从开始就观察这些过程。这里我们考虑了三种表示多状态过程的符号。我们为所研究的过程定义了这种符号，并在后续章节中讨论了此类过程的数据符号；例如，参见第2.2节。对于每个状态最多只能进入一次的进行性过程，有时将t定义为进入状态s的时间是方便的。如果状态是反复出现的，那么可以用t_i表示第i次进入状态s的时间。这种表示的一个问题是，某些状态的进入时间不存在（即，当状态没有被进入时），在这种情况下，相关的分布是不正确的。一种紧凑的替代表示方法是使用t_i^*表示时间t_i时占据的状态s，t^*表示随机多状态过程，其中t_i^*是过程的历史记录——记录了在t_i期间转换的次数、类型和时间。除非另有说明，否则我们在后续内容中假设t_i^*是已知的。计数过程符号提供了第三种强大且方便的多状态过程表示方法——这里我们用t_i^*记录在t_i期间发生的直接转换次数；如果从状态s进入状态t_j，则t_i^*增加1。计数过程符号在现代生命历史过程研究中得到广泛应用，因为它可以轻松适应非常一般的过程，并且与用于证明估计量大样本性质的鞅表示和理论相一致[24]。强度函数是多状态过程的基本构建块，其中强度函数i(t)由以下公式给出：
(1)
其中Δt表示t之前的无限小时间量。因此，它表示在时间t时，给定状态s被占据的情况下，单位时间内的瞬时转移风险，以及记录了在t_i期间所有转换的次数和类型的history。表示(1)非常通用，建模需要明确指定history如何影响转换风险；虽然这是一个强大的框架，但正确指定这种历史依赖性的性质可能是一个艰巨的挑战，但有一些特殊模型适用于广泛的应用。对于马尔可夫过程，强度函数i(t)不依赖于过程历史，除了状态s在t_i时被占据这一事实，所以我们写作i(t) = 1。在早期对马尔可夫过程的使用探索中，Fix和Neyman[25]使用时间齐次的转移强度来模拟癌症复发、死亡和失访。时间齐次模型在早期多状态建模工作中被常规使用，这些模型以及相关的弱参数模型[26]在过程被间歇观察时仍然有用。我们在第2.2节中指出，使用右删失数据可以轻松获得自然非参数估计，并强调在回归设置中的半参数方法。对于半马尔可夫过程，强度依赖于自进入状态s以来的时间，所以我们写作i(t) = λ(t - t_i^*)，其中λ是状态s在t_i期间被进入的总次数，t_i是状态s最近被进入的时间。多状态过程可能涉及一些具有马尔可夫形式的强度，其他一些具有半马尔可夫形式，还有一些可能涉及混合时间尺度。例如，在慢性阻塞性肺疾病中，个体的症状可能会反复发作，随着疾病持续时间的延长，发作的频率可能会增加；发作缓解的时间也可能随着疾病持续时间的延长而增加。在这种情况下，强度可以涉及指定一个基本的时间尺度，并通过回归纳入对时间的其他方面的依赖性。例如，可以采用以下形式的模型：
(2)
其中f(t)是时间的函数，γ是表征对时间依赖性的参数；如果f(t)是时间的函数，则得到一个半马尔可夫模型。连续时间马尔可夫模型是一个值得特别关注的典型模型。如前所述，对于这样的过程，强度函数i(t)是基础构建块，我们定义累积转移强度为I(t) = ∑_{i=1}^N Σ_{j=1}^M f(t - t_i) γ^(t - t_j) δ^(t - t_j)，其中δ^(t - t_j)表示在t_j之前发生转换的概率。表示(1)非常通用，建模需要明确指定history如何影响转换风险；虽然这是一个强大的框架，但正确指定这种历史依赖性的性质可能是一个艰巨的挑战，但有一些特殊模型适用于广泛的应用。对于马尔可夫过程，强度函数i(t)不依赖于过程历史，除了状态s在t_i时被占据这一事实，所以我们写作i(t) = λ(t - t_i)。在早期对马尔可夫过程的使用探索中，Fix和Neyman[25]使用时间齐次的转移强度来模拟癌症复发、死亡和失访。时间齐次模型在早期多状态建模工作中被常规使用，这些模型以及相关的弱参数模型[26]在过程被间歇观察时仍然有用。我们在第2.2节中指出，使用右删失数据可以轻松获得自然非参数估计，并强调在回归设置中的半参数方法。对于半马尔可夫过程，强度依赖于自进入状态s以来的时间，所以我们写作i(t) = λ(t - t_i^*) / (λ(t_i^*) + δ^(t - t_i)，其中λ(t_i^*)是状态s在t_i期间被进入的总次数，δ^(t - t_i)是状态s最近被进入的时间。多状态过程可能涉及一些具有马尔可夫形式的强度，其他一些具有半马尔可夫形式，还有一些可能涉及混合时间尺度。例如，在慢性阻塞性肺疾病中，个体的症状可能会反复发作，随着疾病持续时间的延长，发作的频率可能会增加；发作缓解的时间也可能随着疾病持续时间的延长而增加。在这种情况下，强度可以涉及指定一个基本的时间尺度，并通过回归纳入对时间的其他方面的依赖性。例如，可以采用以下形式的模型：
(3)
其中f(t)是时间的函数，γ是表征对时间依赖性的参数；如果f(t)是时间的函数，则得到一个半马尔可夫模型。连续时间马尔可夫模型是一个值得特别关注的典型模型。如前所述，对于这样的过程，我们定义累积转移强度为I(t) = ∑_{i=1}^N Σ_{j=1}^M f(t - t_i) γ^(t - t_j) δ^(t - t_j)，其中δ^(t - t_j)表示在t_j之前发生转换的概率。考虑到区间[t, t_i]的划分，由t_i和t_j定义，定义I(t) = Σ_{j=1}^M δ^(t - t_j) δ^(t - t_i)。以单位矩阵I为单位，注意I(t)可以被视为在状态s被占据的情况下，表示状态s占用分布的转移概率矩阵的近似；具体来说，I(t, s)的(i)项近似于在状态s被占据的情况下，状态s在t被占据的概率；对角线项I(t, s)近似于在状态s被占据的情况下不发生转换的概率。然后根据划分，乘积
(4)
近似于转移概率矩阵P(t, s)。随着t_i值的增大，这种近似效果更好，如果我们取极限，I(t, s) → P(t, s)。
(5)
其中最后一项只是用来表示乘积积分的符号[24]。获得了具有(i, j)项的转移概率矩阵P(t, s)后，我们现在可以使用这个矩阵来描述多状态过程的特征。当(i, j)在(4)中时，P(t, s)是状态s在时间t被占据的概率。这使我们能够定义一系列边际特征，包括：
i. 过程处于一组状态s中的概率由P(t, s)给出。如果P(t, s) > 1，则对应于过程在时间t之前终止的概率；如果P(t, s) = 1，则对应于过程由于吸收进入特定子集的吸收状态而终止的概率。
ii. 在状态s中预期的总停留时间由T(s)给出，其中停留时间是指在移动到另一个状态之前在特定状态中花费的时间；我们使用“总停留时间”这个术语来表示在给定状态下允许多次停留的时间。可以通过指定积分的有限上限来定义区间[t, t_i]上的受限平均停留时间，表示为T(s, t_i)。
iii. 对于进行性过程，状态s的累积发生率函数由P(t, s)给出，其中s包括状态s以及在状态s停留后可以进入的任何状态。这是在时间t之前进入状态s的概率。这里的“边际”一词指的是在不考虑中间事件或在时间起点之后实现的时间变化协变量的情况下定义的多状态过程的泛函。当数据受到左截断时，直接估计边际特征可能更为复杂。正如我们在第3.1节中讨论的，分析应该基于过程历史来支持独立延迟进入的假设。在这样的条件下估计转移强度后，可以通过涉及模型假设的计算来估计边际特征。接下来我们讨论非参数估计，这在许多实际应用中作为描述性分析是相关的，如第2.2节所示。

2.2 被删失数据、非参数估计和描述性泛函
尽管时间齐次的马尔可夫过程假设转移强度是恒定的，但可以通过允许强度随时间变化来放宽这一假设，同时保持马尔可夫性质。特别是，参数模型可以很容易地扩展以适应转移风险的时间趋势；分段常数强度提供了一个简单实用的公式，可以估计在预先指定的时间间隔内的转移率，特别是在过程被间歇观察时特别有用；参见第4节。这里我们考虑一个单样本问题，其中过程从共同的时间起点开始被连续观察，但受到右删失的影响，对于这种设置，非参数估计相对直接。在下文中，我们考虑一个由最多观察时间t的独立过程组成的单样本，这是一个固定的行政删失时间。我们将其视为计划的固定和共同的随访持续时间，但原则上它可以在个体之间变化。为了适应失访，让t_c表示一个随机的右删失时间，以便个体在区间[t_c, t]内被连续观察。我们在这里假设t_c与过程无关，但将在第2.4节中对此进行更多评论。在一般的多状态设置中定义风险过程和事件计数过程比在单一失效时间设置中更为复杂。函数I(t)表示个体i的过程在时间t是否被观察（即，是否未被删失）。如果N(t)是过程i在时间t内的总转换次数，N(t, s)表示个体i在时间t内的转换次数，I(t, s)是一个指示他们在时间t经历了转换的指标。让I(t, s)表示个体i在时间t处于状态s，为了区分潜在的计数过程和在这种删失方案下观察到的过程，让I(t, s')表示个体i在时间t可能从状态s'转移出去。然后，让I(t, s")表示在时间t为过程i记录的转换次数，N(t, s")表示在区间[t, t]内为过程i观察到的转换总数。如果t_c = t，则对于所有t和s，N(t, s)都为零，因为吸收状态不能被退出。一个自然的估计量I(t, s)由以下公式给出：
(3)
其中N(t, s)是在时间t在样本中观察到的转换总数，N(t, s')是在时间t处于风险中的个体总数。这与生存分析中累积危险度在时间点的非参数估计类似，其形式也是“样本在时间点的事件数量除以时间点的风险集大小”。请注意，我们需要这个估计在时间点是有定义的，按照惯例，当这不满足时我们取零。Nelson-Aalen估计量则是：

(4)

这是[24]区间内不同转移时间的离散和的Stieltjes积分表示。从(3)可以看出，除了观察到转移的时间点外，被积函数将为零。Aalen-Johansen估计量[27]是通过用(4)给出的估计值替换(2)右侧的未知量来获得的转移概率矩阵：

(5)

其中是通过用Nelson-Aalen估计量替换得到的估计累积转移强度矩阵。如果从时间点开始观察到过程，矩阵的顶行包含了状态占用概率的Aalen-Johansen估计量。这反过来使得可以估计函数(i)—(iii)以及其他许多函数。特别需要注意的是，如果我们有一个有两个状态的生存过程，状态为0和1，那么估计量(4)就是累积危险度函数的Nelson-Aalen估计，通过应用(5)我们可以得到生存概率的Kaplan-Meier估计，即矩阵的左上角元素。重要的是，虽然非参数估计量(5)是基于马尔可夫假设的，但这些估计对于非马尔可夫过程也是有效和稳健的，前提是审查完全独立于多状态过程[28, 29]。无穷小刀法在这种设置中提供了一种非常精确的稳健方差估计方法；这在R语言的生存库中得到了实现[7]。这意味着估计量(5)及其基于它的推断可以适用于广泛的过程[?]。

2.2.1 结肠癌研究回顾，I

为了说明这一点，我们考虑了第1.2.1节中介绍的结肠癌数据集。代码可以在支持信息文件的S1部分找到。数据框采用“计数过程”格式，适合分析广泛的多状态过程。在这种格式中，每个个体的随访期被划分为几个间隔，在这些间隔内，个体有可能从一个状态转移到任何其他可能的状态。方程(3)分母中的项要求跟踪个体何时处于不同状态，特别是当他们有从状态转移出去的风险时。图2a显示了复发的累积转移强度的Nelson-Aalen估计以及逐点95%置信区间。这些估计的斜率表明了复发强度（即在存活且无复发的个体中复发的风险）随时间的变化情况。对照组的估计显示前两年的强度大致保持不变，随后是一个强度较低的时期，导致斜率下降。5FU+Lev组的Nelson-Aalen估计在前两年的初始强度较低，之后也趋于平稳。图2b展示了进入死亡状态（和）的Nelson-Aalen估计以及逐点95%置信区间。非常小的估计值反映了无论哪个治疗组，直接从状态0转移到状态2的个体都相对较少。然而，在经历早期复发的个体中，接受5FU+Lev治疗的个体死亡强度有所增加；可以看到在最初的六个月内有明显的增加。对照组的累积强度没有明显增加，因为在考虑的时间范围内，的估计值非常小。然而，对于状态，方程(3)的分母在最初的几个月可能非常小，这可能是一个伪影——在这最初的六个月之后，的斜率大致相似。图2b表明5FU+Lev组的死亡风险略高于对照组。5FU+Lev组术后死亡率的总体降低主要是由于复发风险的降低。Nelson-Aalen估计之间的差异主要是由于在早期随访期间观察到的差异，当时处于风险中的受试者数量相对较少。

结肠癌研究：(a) 复发的累积转移强度的Nelson-Aalen（NA）估计以及逐点95%置信区间。(b) 进入死亡状态（和）的累积转移强度的NA估计以及逐点95%置信区间。(c) 复发的累积发生率的Aalen-Johansen（AJ）估计。(d) 无复发死亡和复发后死亡的累积发生率的AJ估计。如果 是进入状态1的时间，那么复发的累积发生率可以表示为到时间为止进入状态1的概率，由 给出。同样，无复发死亡的累积发生率定义为，而复发后的死亡概率由 给出。请注意，这些都是子分布函数，因为由于存在竞争风险，它们不会随着 而趋近于1。然而，这样的函数具有可解释为概率的吸引人的特点。图2c显示了每个治疗组的Aalen-Johansen估计。显然，5FU+Lev组的复发风险明显较低。图2d展示了与死亡相关事件的累积发生率的估计。显然，两个治疗组的无复发死亡风险都很低，5FU+Lev与复发后死亡风险降低有关。

2.3 基于强度的回归模型

如果感兴趣的是评估时变协变量与多状态过程之间的关系，可以通过基于强度的回归模型[16, 24, 30, 31]来研究。设 表示在 时间点的协变量， 表示协变量过程。如果 是包括协变量路径信息的扩展历史，那么可以通过在方程(1)中用 替换 来修改强度函数。基于强度的回归模型旨在描述 转移的瞬时风险如何依赖于协变量过程的特征。时变协变量可以代表外部因素，如季节、空气污染等。或者可能基于疾病过程的辅助特征；通常感兴趣的是检查标记过程如何与感兴趣的疾病过程相关。例如，在心血管试验中可能会记录胆固醇水平，并且可能感兴趣的是将这些与心脏事件、住院或死亡的风险联系起来。人们还可以定义时变协变量来总结过程历史的重要方面，通过回归来建模历史依赖性。最常见的形式是乘法模型[16, 24, 31, 32]，其中

系数向量 表示协变量对 转移强度的影响。具体来说， 是与第 个协变量增加一个单位相关的 转移的相对风险，而所有其他协变量值保持不变。在生存分析中，通常使用危险比率这个术语，但相对风险是一个更广泛的概念，更适合基于强度的回归分析。如果 ，这被称为调制马尔可夫模型，其中协变量过程调制了基线马尔可夫强度。这个简单模型可以用来研究事件之间的关联；例如，对于疾病死亡过程，如果 与 不同，那么中间事件会影响死亡风险。这可以通过拟合带有约束 的模型来研究。同样，如果 在(6)中，这对应于调制半马尔可夫模型[16]。当协变量是时间固定的，条件是协变量 时，这些模型分别简化为马尔可夫模型和半马尔可夫模型。对于给定协变量的马尔可夫过程，可以定义一个转移概率矩阵 ，其条目为 。如果协变量是离散的，可以通过分层来估计，或者通过拟合像(6)这样的回归模型并应用(2)中的乘积积分来估计。协变量效应也可以表示为对过程强度的加性效应[28]。对于非马尔可夫过程，计算转移概率更具挑战性[33, 34]。在给定协变量的条件下，多状态过程的随机性质完全由转移强度集[24, 30]指定。当协变量是时变的时，协变量和疾病过程的联合模型通常很有用。接下来我们将在第4节讨论在存在时变协变量和随机审查的情况下构建似然。

2.4 基于强度的模型的似然

当过程涉及时变协变量和随访丢失时，重要的是要认识到这些是数据生成中的随机过程。在这里，我们考虑了这一点[16]。设 表示到时间点是否发生随机审查（例如，随访丢失）的指示器，并将相应的计数过程表示为 。关于时变协变量的信息通常在多状态过程进入吸收状态或过程被审查时停止，我们在下文中假设这一点。因此，应该相应地调整风险过程。设 表示个体 在时间点所处的状态是非吸收状态，如果个体 可能在时间点转移，则 等于1。向量 记录了从状态 到 的累积转移次数，其中理解这些计数对于无法直接从状态 进入的状态将为零。最后， 是所有非吸收状态计数过程的向量。我们定义 表示在 上观察到的协变量向量的增量。这里， 确保多状态过程尚未达到吸收状态或被审查，因此可以观察到协变量。此外， 和 。然后，观察到的多状态、协变量和随机审查过程的历史由 表示。随机审查的强度通常定义为

(7)

其中。方程(7)中的项确保一旦多状态过程被审查或达到吸收状态，审查强度为零。为了构建完整的似然，我们考虑由点 定义的 的划分。然后我们考虑在子区间 上的贡献。为此，我们让 表示协变量向量的增量，并让 表示在 上的转移次数。最后， 是划分上审查和联合多状态及协变量过程的历史。对于区间 ，过程 的贡献如下：

请注意，只有当个体没有被审查并且多状态过程在时间点 不处于吸收状态时，才会对 进行似然贡献。其次，只有当个体在时间点 没有被审查时，才会与多状态和协变量过程相关联的贡献。第三，通过采用这种特定的分解，协变量过程的增量的随机模型不仅依赖于 ，还依赖于 和 。这适应了在多状态过程中达到某些（通常是吸收的）状态时协变量可能不再定义的情况。这种公式要求协变量在连续时间内是可用的。这可能会限制可以处理的问题范围，但在许多应用中，协变量在可观察的变化点之间是恒定的。例如，它们记录特定事件是否发生的情况。当仅在间歇性评估时间测量离散时变协变量时，可以形成联合多状态模型，而连续时变协变量可能导致使用联合建模技术[35]。在划分下，基于个体 在 上的数据的完整似然是以下三个项的乘积：

(8)

涉及协变量过程，

(9)

涉及多状态过程，

以及

涉及随机审查过程。如果通过建模审查或协变量过程无法获得关于主要感兴趣的参数（即索引多状态过程的参数）的信息，则称审查和协变量过程为无信息的。因此，除非有兴趣对协变量（通常称为“标记过程”）和多状态过程进行联合建模，并且在假设审查和协变量过程是不提供信息的情况下，否则通常将注意力限制在（9）上。这需要为可观察到的计数过程指定强度函数。我们将（9）中特定区间的贡献概率写为

这可以更明确地写为

（10）
为了进一步进行，需要定义可观察计数过程的强度

（11）
为了用感兴趣过程的强度来表达这一点，我们需要一个额外的假设，即随机审查在给定历史的情况下与多状态过程条件独立[24, 30, 36]。这通常简单地称为独立审查。在这个假设下，（11）中的概率是，我们可以将强度（11）写为。然后，通过用表达（11）并用极限，我们得到

（12）
其中

（13）
其中是观察到的过渡时间集合。这里呈现的似然贡献对应于个体。对于一组独立的过程，总体似然是这些项的乘积。如果元素不共享任何参数，则可以通过分别优化（13）来执行（12）的优化，而这可以使用标准的生存分析软件来完成，前提是它能够处理左截断和右截断的数据。所得估计量的大样本属性直接来自标准生存分析的属性。有关技术细节，请参见Andersen等人[24]，以及Aalen等人[30]、Cook和Lawless[16]以及Andersen和Ravn[37]的相关材料。

2.4.1 结肠癌研究再探，II
为了说明基于强度的回归建模，我们检查了结肠癌研究中的三个风险因素：治疗组（5FU+Lev与对照组，表示为）；侵袭程度，定义为二元变量，包括黏膜下层或肌肉（数据集中的值为1或2）与浆膜或连续结构（3或4）；以及涉及超过4个淋巴结的指标（）。实际涉及的淋巴结数量可能更好地反映疾病阶段，但在这个示例中我们进行了二分。然后我们应用以下形式的基于强度的回归模型

其中和是一个回归系数向量，表示协变量对强度的影响；见图1d。这里，马尔可夫假设意味着，在当前状态和协变量的条件下，转移强度不依赖于之前的过程历史。虽然这对建模很方便，但在许多应用中这个假设可能不现实。可以在支持信息的第S2节中看到的代码中看到，报告了基于模型的标准误差；当指定稳健的标准误差时，也会得到类似的结果。表1和图3的结果显示，当控制疾病侵袭程度和淋巴结参与程度时，5FU+Lev治疗显著降低了复发率。同样，当控制治疗因素时，疾病更广泛的个体和涉及超过4个淋巴结的个体有显著更高的复发率。对于无复发生死，没有证据表明任何风险因素有影响。对于复发后的死亡，当控制疾病程度和淋巴结参与程度时，人们可能会倾向于得出治疗可能有害的结论。然而，我们警告不要这样解释，因为需要更全面地处理可能的时间依赖性混杂因素。如果有更多的可能的时间依赖性协变量，可以进行更细致的分析。在这种情况下，主要目的是调查不同类型的个体在两个组中是否经历复发，如果确实如此，在评估治疗对复发后死亡率的影响时需要考虑这些差异。或者，可以像我们在第3.4节中讨论的那样，基于伪值进行分析。图3在图形查看器中打开

结肠癌研究：基于强度的Cox回归分析得出的相对风险及其相应的95%置信区间。表1. 结肠癌研究：基于强度的分析得出的Cox回归系数（Est）、标准误差（SE）和相对风险（RR）。转移

协变量

Est

SE

RR

进入复发状态，
5FU+Lev
-0.508
0.106
0.603

进入复发状态，扩展程度3或4
0.649
0.168
1.914

进入死亡状态，
5FU+Lev
0.235
0.113
1.265
0.037

进入死亡状态，扩展程度3或4
0.304
0.179
1.355
0.091

进入死亡状态，
5FU+Lev
0.031
0.333
1.035
0.917

进入死亡状态，扩展程度3或4
0.108
0.449
1.115
0.809

3 其他建模考虑
3.1 延迟进入和过程历史的不完整数据
前一节讨论了在理想化情景下的基于强度的建模，其中个体从他们过程的开始就被观察。虽然在启动队列中这是典型的，但在许多研究中，个体在过程已经进行了一段时间后才被纳入。招募通常是一个两步过程：首先，识别符合纳入条件的个体；其次，获得他们参与的同意。我们最初假设有关预选择历史的信息是可用的。让表示个体的招募时间，之后我们打算观察他们的过程，比如说。在某些设置中，例如UK Biobank [38, 39]，有关于某些转移的信息是可用的（例如，癌症诊断），而其他事件可能未被记录（例如，首次诊断高血压）。让和表示个体正在研究中（即，已经被招募并且尚未被审查或进入吸收状态），并且表示个体正在研究中并且在时间有从状态转移的风险。注意我们在中使用，因为他们必须在时间有风险，但对于和，我们使用，因为如果我们要看到任何这样的转移，他们必须在时间被观察。与之前的符号类似，让、和。对于时间固定的协变量，包括延迟进入时间的信息的扩展历史可以写为。这假设有关于之前的过程的信息是可用的，从而能够对进行强度函数的建模。在条件独立的延迟进入[40]和条件独立的随访丢失的情况下

然后可以像（12-13）中那样构建似然，但是用替换，并且是观察到的过渡时间集合。当审查是条件独立的但这种等式不成立时，有证据表明延迟进入是相关的。处理相关的延迟进入比处理相关的随访丢失更具挑战性，因为未选择个体的信息可能是未知的。获取代表性样本或人口数据可以帮助衡量任何偏差的程度，并提供减轻选择偏差影响的途径[39, 41]。在信息关于完全缺失或高度粗略化的情况下[42]，拟合严重依赖于过程历史的基于强度的模型变得具有挑战性。有必要努力获取这些历史，否则将需要更强的简化建模假设。尽管马尔可夫模型应该基于科学合理性和充分性的证据来证明其合理性，但在这种情况下它们特别有吸引力，因为强度在给定当前状态的情况下与历史无关。
3.2 时间依赖性协变量和联合建模
第2.4节中给出的似然分解证明了仅基于多状态过程的似然使用的合理性。然而，在感兴趣的是两个过程之间的关系时，协变量和多状态过程的联合建模具有科学价值。例如，在研究骨骼健康标志物与骨折风险之间的关系时，可以将连续的骨骼形成和吸收标志物纳入首次和后续骨折发生的强度中。反过来，骨折可以影响骨骼标志物，这种效应可以在两个过程的联合模型中进行检查[43]。在这种情况下，可以考虑基于（8）和（9）的似然。如果连续的骨骼标志物被离散化，可以构建一个联合的多状态模型，状态由标志物水平和骨折状态的组合定义，可能包括一个吸收状态（代表死亡）。在这个例子中，当协变量受到间歇性观察时，会出现额外的挑战。我们将在第4节讨论如何解决这个问题。
3.3 边缘参数和伪值的推断
强度函数是多状态过程的基本组成部分，如上所述，指定所有强度函数可以构建似然。这进一步意味着过程的所有概率方面都是确定的——至少当强度模型不包括引入“额外随机性”的时间依赖性协变量时（即，当基于方程（8）的似然是直接的时候）。因此，可以基于估计的强度函数来估计边缘特征，例如状态占用概率和状态中的预期停留时间，无论是通过“插入”方法（如果可以指定数学关系）还是通过模拟。然而，在回归设置中，插入方法并不提供直接描述时间固定协变量之间关联的参数。此外，如果主要的科学兴趣在于这种关联，那么通常需要建模所有强度函数，模型错误指定就成为一个问题。因此，直接指定一个边缘模型是有趣的，即不依赖于强度函数的模型。一般来说，不可能以这样的方式指定多状态过程的强度函数，以便简单的边缘模型（如（14）成立。因此，直接边缘模型应该被视为一个“工作模型”，有助于评估治疗对复发后死亡率的影响，但并不一定反映真实的数据生成机制。对于任何回归模型，应该通过适当的诊断仔细评估简化假设，如可加性和线性。在这里，我们讨论了两种直接的边缘建模方法：使用逆概率加权（IPCW）的直接二项回归和广义估计方程（GEE）[44, 45]以及伪值（PV）方法[46, 47]；另见Andersen和Ravn[37]的最新书籍。
3.4 其他建模考虑
3.1 迟延进入和过程历史的不完整数据
前一节讨论了在理想化情景下的基于强度的建模，其中个体从他们过程的开始就被观察。虽然这在启动队列中是典型的，但在许多研究中，个体在过程已经进行了一段时间后才被纳入。招募通常是一个两步过程：首先，识别有资格纳入的个体；其次，获得他们参与的同意。我们最初假设有关预选择历史的信息是可用的。让表示个体的招募时间，之后我们打算观察他们的过程，比如说。在某些设置中，例如UK Biobank [38, 39]，有关于某些转移的信息是可用的（例如，癌症诊断），而其他事件可能未被记录（例如，首次诊断高血压）。让和表示个体正在研究中（即，已经被招募并且尚未被审查或进入吸收状态），并且表示个体正在研究中并且在时间有从状态转移的风险。注意我们在中使用，因为他们必须在时间有风险，但对于和，我们使用，因为如果我们想看到任何这样的转移，他们必须在时间被观察。与之前的符号类似，让、和。对于时间固定的协变量，包括延迟进入时间的信息的扩展历史可以写为。这假设有关于之前的过程的信息是可用的，从而能够对进行强度函数的建模。在条件独立的延迟进入[40]和条件独立的随访丢失的情况下

然后可以像（12-13）中那样构建似然，但是用替换，并且是观察到的过渡时间集合。当审查是条件独立的但这种等式不成立时，有证据表明延迟进入是相关的。处理相关的延迟进入比处理相关的随访丢失更具挑战性，因为未选择个体的信息可能是未知的。获取代表性样本或人口数据可以帮助衡量任何偏差的程度，并提供减轻选择偏差影响的途径[39, 41]。在关于的过程信息完全缺失或高度粗略化的设置中[42]，拟合严重依赖于过程历史的基于强度的模型变得具有挑战性。努力获取这些历史是必要的，否则将需要更强的简化建模假设。尽管马尔可夫模型应该基于科学合理性和充分性的证据来证明其合理性，但在这种情况下它们特别有吸引力，因为强度在给定当前状态的情况下与历史无关。
3.2 时间依赖性协变量和联合建模
第2.4节中给出的似然分解证明了仅基于多状态过程的似然使用的合理性。然而，在感兴趣的是两个过程之间的关系时，协变量和多状态过程的联合建模具有科学价值。例如，在研究骨骼健康标志物与骨折风险之间的关系时，可以将连续的骨骼形成和吸收标志物纳入首次和后续骨折发生的强度中。反过来，骨折可以影响骨骼标志物，这种效应可以在两个过程的联合模型中进行检查[43]。在这种情况下，可以考虑基于（8）和（9）的似然。如果连续的骨骼标志物被离散化，可以构建一个联合的多状态模型，状态由标志物水平和骨折状态的组合定义，可能包括一个代表死亡的吸收状态。在这个例子中，当协变量受到间歇性观察时，会出现额外的挑战。我们将在第4节讨论如何解决这个问题。
3.3 边缘参数和伪值的推断
强度函数是多状态过程的基本组成部分，如上所述，指定所有强度函数可以构建似然。这进一步意味着过程的所有概率方面都是确定的——至少当强度模型不包括引入“额外随机性”的时间依赖性协变量时（即，当基于方程（8）的似然是直接的时候）。因此，可以基于估计的强度函数来估计边缘特征，例如状态占用概率和状态中的预期停留时间，无论是通过“插入”方法（如果可以指定数学关系）还是通过模拟。然而，在回归设置中，插入方法并不提供直接描述时间固定协变量之间关联的参数。此外，如果主要的科学兴趣在于这种关联，那么通常需要建模所有强度函数，模型错误指定就成为一个问题。因此，直接指定一个边缘模型是有趣的，即不依赖于强度函数的模型。一般来说，不可能以这样的方式指定多状态过程的强度函数，以便简单的边缘模型（如（14）成立。因此，直接边缘模型应该被视为一个“工作模型”，有助于评估边缘参数和之间的直接关联，但不一定反映真实的数据生成机制。对于任何回归模型，应该通过适当的诊断仔细评估简化假设，如可加性和线性。在这里，我们讨论了两种直接的边缘建模方法：使用逆概率加权的直接二项回归（IPCW）广义估计方程（GEE）[44, 45]和伪值（PV）方法[46, 47]；另见Andersen和Ravn[37]的最新书籍。我们通过研究固定时间点的来说明这些方法，但强调类似的方法也可以应用于多个时间点的联合推断，或者用于状态中的受限平均停留时间。此外，可以使用这些方法通过地标法[48-50]研究条件概率，如或。考虑一个回归模型

（14）
其中，是一个指定的链接函数，系数向量包括特定于时间点的截距。因此，系数将特定于状态和时间点，尽管为了方便表示，我们将系数向量表示为。典型的链接函数包括cloglog，对应于两状态模型中的比例风险模型（图1a），或logit函数。直接二项回归基于在时间观察到状态指标的受试者。这些是那些在时间之前必须发生或的受试者（即，或达到吸收状态的时间）必须在随机审查的时间之前发生。然后在这些受试者的状态指标中用作GEE的响应，其中

（15）
每个项都有一个权重，反映了未被审查的概率，通常包含对的部分导数。显然，这种方法需要一个随机审查的模型，在其最简单的形式中，得到的权重可以由Kaplan-Meier估计器给出。然而，如果协变量影响审查，则需要一个回归模型来估计权重。（15）中的项是独立的，通常使用三明治估计器来估计（15）解的方差，表示为，其中包含来自需要估计的贡献[45]。边际回归模型（14）也可以使用PV进行分析。采用这种方法，每个观测值的结果变量是通过边际状态占用概率的基估计器计算得出的（忽略协变量），用表示。即使不假设多状态过程是马尔可夫的，Aalen-Johansen估计器也是一致的[29, 51]。受试者的PV由下式给出：

（16）

其中是应用于从完整样本中移除受试者后得到的样本大小的（Aalen-Johansen）估计器。其直觉是量化了基估计器受到受试者数据影响的程度，在没有审查的特殊情况下（Aalen-Johansen估计器简化为时间在状态的过程的相对频率），则简单地等于[37]。注意，即使观察到了，也是通过对所有受试者使用（16）来计算的。然后PV被用作GEE中的响应变量，其中

（17）

（17）中的项不是独立的[52]，需要特殊技术来评估大样本属性，作为的解。这些取决于将观测到的多状态过程的数据映射到的功能的影响函数的属性[52, 53]。这些属性成立的一个必要条件是审查不依赖于协变量。如果协变量影响审查，则（16）中的Aalen-Johansen基估计器可以被IPCW估计器替代[54]。应当注意的是，当基估计器基于延迟进入的数据时，通常不满足影响函数的所需属性[55]。（17）中项之间缺乏独立性的另一个后果是，标准的GEE三明治估计器用于的方差应该被一个也涉及功能的二阶影响函数的校正估计器替代。然而，在实际应用中，校正项往往很小[52]。尽管在基估计器是Kaplan-Meier的特殊情况下已经进行了这样的比较[56]，但尚未基于（15）或（17）对估计器和进行系统比较。在大样本中，使用（16）计算PV可能非常耗时。通过R中的survival包实现的无限小刀法PV方法的近似提供了一个更有效的替代方案。此外，对于某些特定的多状态模型，如图1a、b和f所示，可以直接为所有时间点建模。这些模型每个都需要专门的估计方程，例如基于部分似然原理的方程[57-60]。

4 连续时间过程的间歇性观察

在许多情况下，状态之间的转换并不直接观察到，只记录在间歇性评估时间所处的状态。例如，在视网膜病变研究中，视力在诊所访问期间进行评估[61]；在糖尿病肝病学中，通过血液测试或活检评估肝功能[62]；在骨质疏松症研究中，定期放射学检查可以检测无症状的椎体骨折[63]。为了适应似然构建中的间歇性观察，我们考虑了多状态过程以及时间独立的协变量。评估过程由一个计数过程表示，它记录了直到时间点的评估次数。在时间点的评估结果为，否则为。由于访问只发生在仍在研究中的个体身上，因此评估过程在终止，所以我们观察到和。如果，则评估过程的强度为

其中。这种一般公式涉及对的依赖，但这个过程并没有被完全观察到。在这种情况下，必须指定疾病和访问过程的联合模型。这些模型通常在给定具有假设分布的潜在变量的条件下构建条件独立性。Lange等人[64]和Cook和Lawless[65]提出了考虑局部依赖性的联合模型，其中访问强度可能取决于在时间点所处的状态，并讨论了专注于仅涉及多状态过程强度的部分似然贡献所需的独立性条件；另见Grüger等人[66]。假设个体在时间点有次访问，并让表示个体在时间点的观察历史。如果访问过程是非信息性的，意味着访问和多状态模型之间没有共享参数，我们可以忽略访问过程，专注于以下形式的似然

（18）

如前所述，对于一般过程来说，用强度函数表示可能具有挑战性，但马尔可夫模型相对容易处理。例如，如果对于所有转换，我们可以构建转换强度矩阵，其中非对角线元素由给出，对角线元素为。然后，转换概率矩阵有元素。Kalbfleisch和Lawless[67]讨论了这种模型的Fisher评分算法，该算法以及其他优化方法在R的msm[23]包中实现。可以放宽时间均匀基强度的假设，允许它们在指定切点的数量和位置后是分段常数率。Titman[68, 69]考虑了具有阶段型分布的逗留时间的半马尔可夫模型[70]；另见Yang等人[71] Satten[72]考虑了具有共同乘法随机效应的渐进模型，该效应同时适应逗留时间的序列依赖性和个体间进展率的异质性。为了适应更一般的异质性形式，还开发了扩展，包括更高维的随机效应[73]和移动者-留居者组分[74]。随机效应（或脆弱性）模型对于间歇性观察的过程特别有用，其中怀疑存在历史依赖性，但缺乏详细信息，使得直接依赖性建模具有挑战性。转换信息通常通过双重观察方案获得。例如，在模拟认知和生存的痴呆症研究中，认知状态仅在评估期间观察到，而生存时间则连续记录。图1c中显示了一个简单的模型来说明这种设置。如果状态1的评估过程涉及随机访问时间，则进入状态1的时间是区间审查的。此外，在一些研究中，可能不确定是否发生了转换。例如，如果疾病在最后一次访问时没有被诊断出来，那么不清楚它是在那次访问和死亡时间之间发展的。因此，必须调整似然以考虑这种不确定性。在这种情况下，Leffondre等人[75]探讨了当主要兴趣是总体生存时，疾病-死亡模型的有用性。Joly等人[76]开发了用于拟合此类数据的基于强度的模型的方法，使用基于样条的方法来建模强度或Weibull形式[77]。当这些变量不是连续测量时，包括时间依赖的解释变量需要额外的假设。因此，这种情况通常使用时间依赖但区间恒定的协变量来解决，其中值变化的时间是已知的。Boruvka和Cook[78]检查了可识别性问题，并应用了筛选最大似然方法来估计转换强度和协变量的效应。更一般地，Commenges等人[79]探讨了在间歇性观察下使用样条估计多状态过程的方法。有关多状态模型中时间依赖协变量的进一步讨论，请参见[16, 31]。

4.1 痛风关节炎数据的重新审视

在这里，我们分析了第1.2.2节中的痛风关节炎联合损伤数据。该登记处的患者安排每年进行临床检查和每两年进行一次放射学检查，但不同患者之间的访问时间差异很大。图4a显示了五名患者的原始数据，水平线代表从第一次诊所访问到失去随访或死亡的时间段。垂直虚线标记表示放射学检查，突出了个体之间影像数据收集频率的显著差异。我们考虑了Gladman和Farewell[22]提出的四状态模型，如图1b所示。具体来说，状态0对应于没有受损关节，状态1对应于1-4个受损关节，状态2对应于5-9个受损关节，状态3是一个吸收状态，表示10个或更多受损关节。时间尺度从痛风关节炎诊断的年龄开始。假设访问过程是非信息性的，我们使用似然（18）在马尔可夫模型下估计分段常数转换强度，切点分别在疾病发作后的5年、10年和20年。图4b和表2基于支持信息文件S3中提供的代码。图4在图形查看器中打开

图4多伦多大学痛风关节炎队列的联合损伤数据：基于拟合的马尔可夫模型的状态占用估计，这些模型的基线强度在痛风关节炎发作后的5年、10年和20年处有切点。(a) 五名个体的UTPC访问时间示意图。(b) 状态占用的估计。表2. 基于拟合的多伦多大学痛风关节炎队列关节损伤数据的乘法强度基马尔可夫模型的估计，基线强度在痛风关节炎发作后的5年、10年和20年处有切点。

(a) 估计的回归系数（Est）和相对风险（RR）

转移
协变量
Est
SE
RR
95% CI

积液
0.742
0.400
2.100
(0.960, 4.597)
0.063

升高的ESR
0.239
0.278
1.271
(0.737, 2.188)
0.389

积液
0.536
0.297
1.710
(0.955, 3.062)
0.071

升高的ESR
0.774
0.281
2.169
(1.250, 3.759)
0.006

积液
0.306
0.311
1.358
(0.739, 2.497)
0.325

升高的ESR
-0.359
0.364
0.698
(0.342, 1.425)
0.323

(b) 估计的转换强度（Int），基线是将协变量设置为0

基线
Int
95% CI
Int
95% CI
Int
95% CI

0.092
(0.052, 0.161)
0.718
(0.387, 1.333)
0.419
(0.187, 0.939)
1.566
(0.763, 3.215)

0.097
(0.048, 0.198)
0.828
(0.394, 1.743)
0.796
(0.406, 1.559)
1.126
(0.502, 2.529)

0.244
(0.072, 0.821)
1.326
(0.434, 4.048)
1.174
(0.381, 3.618)
1.373
(0.431, 4.375)

图4b是通过使用msm()包最大化似然得到的（同时忽略协变量），其中每个是具有在疾病诊断后5年、10年和20年切点的分段常数强度的参数向量。使用这些估计值，我们计算并在疾病发作后的30年内绘制。随着患者病情恶化，过渡状态的状态占用概率上升然后下降，预计大约50%的患者在20年后会有10个或更多受损关节。表2总结了基于强度的回归模型的结果，包括基线协变量——这是广泛积液的指标和升高的红细胞沉降率（ESR）的指标，后者是炎症的标志。表2a显示了每个转换的回归系数，表2b提供了每个时间间隔的估计转换强度。有趣的是，基线时升高的ESR是更快转换的强预测因子，但它对其他两个转换没有显著影响。

5 具有脆弱性或copula方法的多状态模型

有时，即使考虑了解释个体间疾病进程变异的协变量，生命轨迹的变异性仍然比基于模型假设的预期要大。为了解决这个问题，可以考虑具有潜在随机效应的模型——在生存分析中，这些随机效应被称为脆弱性[80]。脆弱性通常被引入来模拟依赖性，但依赖性也可以使用copula函数来建模，copula函数是具有均匀边际分布的多元分布函数[81]。通过相应失效时间变量的概率积分变换，可以从copula函数构建更一般的联合分布[82]。已经有很多关注使用基于脆弱性和copula的模型来模拟个体或群体内的依赖性。在多状态生存模型的框架内，有两个关键场景是相关的：

i. 在受试者内部的依赖性：在这里，随机效应或copula模型解释了同一受试者内事件时间之间的未观察到的协变量。例如，在图1d中，随机效应可能解释了和之间的未观察到的依赖性。

ii. 在受试者之间的依赖性：在这种情况下，如家庭或研究中心这样的聚集数据涉及同一集群内个体之间的相关失效时间。随机效应或copula模型可以解决这种未观察到的依赖性。以下将讨论这两种情况，我们专注于脆弱性模型。在一般的多状态框架中应用脆弱性或copula模型可能很复杂，因为未观察到的异质性可能在不同转换之间变化。因此，为了简化问题，我们主要关注疾病-死亡模型（图1c或d）。如第5.1节和第5.2节所讨论的，受试者内部和受试者之间的依赖性是一个重要的领域，需要进一步的方法学研究。

5.1 受试者内部随机效应

例如，考虑一种慢性疾病的疾病-死亡模型，该模型涉及三种可能的转变：从健康到疾病（），从健康到无疾病死亡（），以及从疾病到疾病后的死亡（）。诊断后的生存时间被诊断时间所截断，假设观察到的协变量能够捕捉诊断时间和疾病后死亡时间之间的所有依赖性通常是不现实的。这促使我们引入一个未观察到的受试者特定随机效应，以解释受试者事件时间之间的剩余依赖性。以下讨论了两种不同的方法：一种方法关注在未观察到的随机效应条件下观察到的协变量的回归系数，另一种方法则不考虑这些随机效应。每种方法对观察到的协变量效应的解释不同，具体细节如下。这些方法将使用Cox模型（第5.1.1节）和加速失效时间模型（第5.1.2节）进行演示。软件实现将在第6.7节中概述。用于失效时间结果的随机效应模型通常被称为脆弱性模型，在这种模型中，一个代表与失效时间相关的未观察到的风险因素的随机受试者特定效应被称为脆弱性。Xu等人[83]提出了一个单参数伽马脆弱性模型，在该模型中，一个脆性变量贯穿疾病-死亡模型的所有三种转变（详见下面的方程（19-21））。这种建模策略在后续工作中被广泛采用[39, 84-88]，并且提供了一个简洁的框架；然而，假设所有转变都具有共同的未观察到的风险结构可能并不适用于每个应用。在涉及额外转变的更复杂的多状态模型中，假设一个单一的脆性变量可能过于简化了潜在的异质性。其他替代方法包括使用受试者特定的向量值脆弱性项[16, 89]，或者对共享的脆弱性成分进行转变特定的转换[89]。大多数现有的基于脆弱性的多状态模型方法都集中在疾病-死亡和渐进（向前）过程上，如图1b所示，并包括适应竞争风险和层次聚类结构的扩展[90]。脆弱性方法的一个关键特点是假设，在给定观察到的协变量和脆弱性变量的情况下，各个转变时间是独立的（或准独立的[39]）。这一假设允许对每个转变进行单独建模。此外，通常还假设未观察到的脆弱性变量和观察到的协变量是独立的。然而，为了确保模型的可识别性，需要某些假设[91]，并且用现有数据验证这些假设可能具有挑战性。在处理不可观察的随机效应时，区分两种方法是重要的：条件建模，它同时考虑观察到的协变量和未观察到的脆弱性变量；以及边际（总体平均）建模，它仅考虑观察到的协变量。在线性模型中，如果观察到的协变量与脆弱性变量独立，这两种方法会得到相同的估计结果。然而，在非线性模型中则不是这样，因此这种区分在实际应用中非常重要。选择这两种方法中的哪一种取决于分析的具体目标。在具有脆弱性的多状态模型中，时间尺度的选择至关重要。考虑两个疾病-死亡模型：在第一个模型中，状态是“健康”（状态0）、“疾病”（状态1）和“死亡”（状态2和）。在第二个模型中，状态是“手术”、“复发”和“死亡”。在第一种情况下，预期诊断时的年龄与从诊断到死亡的时间之间存在负相关，因为年龄较大的个体在诊断后的寿命通常较短。因此，鉴于我们打算使用共享的随机效应来捕捉这三种转变之间的相互作用，自然会使用相同的年龄尺度来表示所有三种转变。在第二种情况下，需要采用时钟重置方法（即半马尔可夫模型），其中时间在每次转变时重置，以防早期复发与剩余寿命呈正相关。在这里，每次转变的时间尺度取决于在前一个状态中度过的时间。

5.1.1 受试者内部随机效应：条件化与边际化的疾病-死亡Cox模型

由于图1d中的疾病-死亡模型描述了一个每个状态最多访问一次的渐进过程，我们采用了第1节中的简化符号，对于独立的观察值，用 和 分别表示非终止事件（例如，疾病诊断）和终止事件（例如，死亡）的时间。 的联合分布支持在 上。对于那些在非终止事件之前经历终止事件的个体，我们设置 。用 表示右删失时间，用 表示一组时间固定的协变量。虽然这里讨论的一些模型可以很容易地扩展以包含时间依赖的协变量，但其他模型则需要大量的额外修改。定义 。用 表示 是否对应于非终止事件，用 表示是否对应于终止事件。同时，用 表示终止事件或删失的时间，当观察到非终止事件时为1，否则为0。用 表示在非终止事件之后是否观察到了终止事件。观察到的数据是 。我们考虑一个标量潜在脆弱性变量 对于个体 ，其累积分布由未知参数 索引。Xu等人[83]提出了一个基于Cox函数的疾病-死亡模型，其中包含了一个伽马分布的共享脆弱性变量，该变量对每个强度函数具有乘法影响。这适应了非终止事件和终止事件时间之间的未观察到的依赖性。在时间独立的协变量 和 的情况下，控制三种转变的条件强度函数表示为

（19）
（20）
（21）

其中 和 ， 分别是特定于转变的基线危险函数和回归系数。鉴于受试者在时间 被诊断（即进行了 转变），诊断后的死亡时间 被 裁断。非终止事件的时间 不被纳入 的协变量向量中。相反，潜在事件时间 和 之间的依赖性来自两个关键因素：所谓的解释性危险比率 [83] 和潜在脆弱性分布参数 。解释性危险比率表征了脆弱性未捕捉到的非终止事件和终止事件时间之间的局部依赖性。如果 和 是独立的，解释性危险比率恒为1，脆弱性变量也恒为1（即脆弱性分布是简并的）。在均值为一且方差未知的伽马分布脆弱性的情况下，模型（19-21）的估计可以通过半参数最大似然估计器（MLE）进行，其中似然是通过在 的分布上对给定 的观察数据的条件似然进行平均得到的[83]。也可以使用半参数贝叶斯方法[84]。Jiang和Haneuse [87]以及Lee等人[92]对（19-21）进行了有趣的修改。基于条件危险（即给定 ）的生存预测需要知道未观察到的脆弱性变量，这限制了它们的实际应用性。通过积分掉脆弱性，我们得到了相对于 的所谓的边际化危险。这些边际化危险强烈依赖于假设的脆弱性分布及其参数。尽管这种方法允许在没有直接观察到的情况下进行预测，但它引入了对脆弱性错误指定的敏感性，并使协变量效应的解释复杂化。为了解决这些问题，Gorfine等人[39]提出了一种基于脆弱性的疾病-死亡模型的替代策略，其中使用Cox模型对相对于 的边际危险进行建模，并为指定的脆弱性分布推导出产生这种比例规格的条件危险。这些条件危险依赖于 和边际危险的参数，偏离了比例危险结构。这种方法允许在纳入脆弱性以解释未观察到的受试者特定协变量的同时估计边际模型参数。具体来说，给定 的疾病-死亡模型的条件危险表示为

（22）
（23）

相对于 的相应边际化危险函数定义为

（24）
（25）

为了保持边际化危险和条件危险函数之间的关系，必须根据假设的脆弱性分布将条件危险映射到它们的边际化对应物。实际上，非负函数 由脆弱性 的分布和相应的边际化危险 决定。例如，在期望值为1且方差为 的伽马脆弱性模型下，可以证明 ， ，其中 ， ， 和 。总之，上述基于脆弱性的模型（22）和（23）偏离了比例危险结构，它们用模型（24）和（25）的主要参数 和 来表示。基于脆弱性的模型特别适合开发估计程序，因为它们假设在给定观察到的协变量和脆弱性变量的情况下， 和 是准独立的。在标准脆弱性假设下，开发了一种伪似然方法来估计参数及其标准误差[39]。

5.1.2 受试者内部随机效应：加性 versus 乘性疾病-死亡AFT模型

与基于强度的多状态过程建模相比，加速失效时间（AFT）回归技术的适应是一个重要的替代方案。与基于乘性强度的模型不同，AFT模型直接在时间尺度上参数化协变量效应，提供了对协变量效应的独特且通常更易于解释的描述[93]。在竞争风险、复发事件和多变量失效时间设置中，AFT模型通常根据不同的转变对应的潜在失效时间来构建，承认由于竞争事件的发生，其中一些潜在事件时间可能不会被观察到；例如，参见[94-97]。这种潜在时间的构建已经得到了很好的确立。在本节中，我们回顾了疾病-死亡过程的AFT模型的最新发展，特别强调了纳入随机效应以解释未观察到的异质性的扩展。Lee等人[86]提出了以下加性尺度变化模型

其中 的第一个分量中的值为1，以允许截距。术语“加性”指的是在所有转变的对数线性规范中加入了一个共享的脆弱性项 。该模型是根据潜在的、特定于路径的失效时间来定义的，观察到的终止时间由实现的转变路径决定。该模型与典型的疾病-死亡进展一致：一个受试者首先经历一个事件，其类型决定了随后的转变。如果失败对应于非终止事件，那么随后的转变时间 根据从状态1到状态 的路径进行建模。随机误差 在 和 之间是独立的，并且对于 是特定于转变的，可能具有未指定的分布。 和 之间的关联由 的分布决定，表现为对数失效时间尺度中的加性成分。鉴于 假设遵循正态分布，既有参数化的也有半参数化的贝叶斯估计方法可用[86]。另一种方法是提出了基于乘性脆弱性的AFT模型[88]，在这种模型中，未观察到的脆弱性变量没有直接在对数失效时间线性模型中表达。相反，它影响随机误差 的分布。特别是，该模型由以下公式定义

通过以下共享脆弱性模型纳入 和 之间的依赖性。给定个体 的脆弱性变量 ，假设相应的条件基线危险函数为

其中每个 是 的未指定的基线危险函数。重要的是，该模型清晰地区分了观察到的协变量对事件时间的影响和通过脆弱性捕获的未观察到的异质性。在假设 服从均值为1且方差未知的伽马分布的情况下，可以使用半参数MLE（基于核平滑似然与EM算法的结合）[88]。Kats和Gorfine [88]的第2.2节深入探讨了上述加性和乘性方法之间的概念区别。阐明了加性方法的危险可能相对于 显示出非单调行为。相比之下，乘性模型的危险作为 的函数显示出单调增加。因此，乘性脆弱性模型提供了更简单的解释。

5.1.3 重新审视鹿特丹肿瘤银行数据

Kats和Gorfine [88]使用三种基于脆弱性的模型分析了鹿特丹肿瘤银行数据（第1.2.3节）：乘性AFT模型[88]、边际化Cox模型[39]和条件Cox模型[84]。添加剂-脆弱性AFT模型[86]在R包SemicompRisks中实现，用于模拟在复发后经历死亡的受试者的时间。然而，当将该模型应用于Rotterdam数据时遇到了收敛性问题（更多细节见[88]）。在当前的分析中，我们通过比较八个模型来扩展这一评估：四个基于Cox框架的模型和四个基于AFT框架的模型。基于Cox的模型包括边缘化Cox模型[39]以及以下三个不考虑脆弱性的模型：

Cox模型（无脆弱性I）：为每个转换分别拟合三个独立的Cox模型。对于竞争风险转换（复发和死亡），通过将竞争事件视为右删失来进行处理。对于复发时间的左截断，使用标准的风险集调整方法来处理。Cox模型（无脆弱性II）：该模型通过在转换模型中添加标准化的复发时间作为时间独立的协变量来扩展模型I。具体来说，让表示受试者的复发时间（进入状态1），并定义，其中和是的平均值和标准差，即在观察到复发的个体中计算得出。这种中心化和缩放将复发时间放在一个可比较的数值尺度上，有助于解释回归效应。包含这一变量遵循[98]的方法，旨在解决由复发引起的依赖性左截断问题。Cox模型（无脆弱性III）：作为模型II的扩展，该模型包括一个具有四个节点（在第20、40、60和80百分位数）的线性截断样条，以灵活地模拟标准化复发时间的非线性效应。具体来说，它包含了形式为的项，其中，，是回归系数。

在当前的分析中，我们通过比较八个模型来扩展这一评估：四个基于Cox框架的模型和四个基于AFT框架的模型。基于Cox的模型包括边缘化Cox模型[39]以及以下三个不考虑脆弱性的模型：

Cox模型（无脆弱性I）：为每个转换分别拟合三个独立的Cox模型。对于竞争风险转换（复发和死亡），通过将竞争事件视为右删失来进行处理。对于复发时间的左截断，使用标准的风险集调整方法来处理。Cox模型（无脆弱性II）：该模型通过在转换模型中添加标准化的复发时间作为时间独立的协变量来扩展模型I。具体来说，让表示受试者的复发时间（进入状态1），并定义，其中和是的平均值和标准差，即在观察到复发的个体中计算得出。这种中心化和缩放将复发时间放在一个可比较的数值尺度上，有助于解释回归效应。包含这一变量遵循[98]的方法，旨在解决由复发引起的依赖性左截断问题。Cox模型（无脆弱性III）：作为模型II的扩展，该模型包括一个具有四个节点（在第20、40、60和80百分位数）的线性截断样条，以灵活地模拟标准化复发时间的非线性效应。具体来说，它包含了形式为的项，其中，，是回归系数。

为了评估模型拟合度，我们使用了基于随机生存概率（RSPs）[88, 99]的视觉拟合优度（GOF）诊断方法。具体来说，对于疾病-死亡模型，评估了两个RSPs：（i）保持在状态0的概率，表示为；以及（ii）在经历非终端事件的个体中保持在状态1的概率，表示为。这些RSPs是根据观察到的协变量为每个观测值估计的；对于基于脆弱性的模型，通过积分脆弱性分布来获得边际概率。在正确指定的模型下，这些RSPs应该遵循一个均匀分布。因此，可以通过将估计的RSPs的直方图与标准均匀分布进行比较来视觉评估模型拟合度。在考虑的八个模型中，基于乘法脆弱性的AFT模型似乎提供了对数据的最佳整体拟合。总之，尽管Cox和AFT模型对大多数协变量提供了一致的方向性解释，但AFT框架提供了基于时间的协变量效应理解，这在临床背景下可能更直观。此外，考虑脆弱性的模型似乎更好地捕捉了该数据集中的复杂性和异质性，从而提高了拟合度和可解释性。

5.2 受试者间依赖性

在家庭或双胞胎研究等聚类数据的背景下，对疾病-死亡模型的探索有限。基于脆弱性的模型和估计方法已被开发用于处理聚类失效时间数据中的竞争风险（即图1f）[77, 100-102]。然而，将这些方法扩展到更复杂的多状态模型环境中是具有挑战性的，主要是由于不同转换状态下两个簇成员之间的依赖强度各不相同。Lee和Cook [92]首次使用潜在变量公式开发了一个疾病-死亡模型来处理第一个转换，或者采用copula模型来适应簇内从状态0转换到状态1的（可能是潜在的）时间内的依赖性。

6 软件可用性
本节概述了可用于分析多状态模型的软件和包（无论是R语言还是Python），并强调了每个包的独特贡献。

6.1 R语言包survival和mstate
survival包[103, 104]中的survfit函数[32]可以计算并绘制任何多状态情景下的Aalen-Johansen累积发生函数估计值。特别是，研究对象在研究过程中可以访问多个状态，可以在时间0之后开始（即延迟进入），并且可以从任何状态开始。Aalen-Johansen估计值的标准误差是使用无穷小刀法计算的。coxph函数为多状态模型中的每个转换提供Cox回归分析，无论是否有共享系数。mstate包[105, 106]包含了数据准备、描述性分析、危险函数估计、使用Aalen-Johansen估计值进行预测以及Cox回归建模的实用工具。由于其模块化方法，可以拟合不同的模型，例如加性危险模型，同时仍然允许基于Aalen-Johansen进行预测。此外，还提供了测试马尔可夫假设的函数[107]。

6.2 R语言包msm
msm包[23]提供了一套用于模拟和分析具有分段常数强度函数的连续时间马尔可夫过程的功能，以及包括未观察到的（隐藏的）状态的隐马尔可夫模型；后者在这里不进行讨论。该包的一个独特特点是，当过程受到右删失或间歇性观察时，它便于对具有分段常数基线强度的马尔可夫过程进行比例强度回归分析。因此，该包非常适合用于建模第4节中描述的数据类型，如第4.1节中的示例所示。该包提供了转换强度、转换概率矩阵和在每个状态中预期时间的估计值。它还提供了对未来状态占用的预测。参数通过最大似然法估计，从而能够使用标准的推断技术来计算置信区间和进行假设检验。

6.3 Python包PyMSM
PyMSM[108]是一个用于拟合竞争风险和多状态模型的包，提供灵活的模型指定、个体和群体级别的预测，以及全面的统计摘要和可视化。关键特性包括：(1) 多状态回归模型拟合：支持各种生存分析技术，如Cox回归、随机生存森林[109]或用户定义的机器学习模型。(2) 通过蒙特卡洛模拟进行预测：使用拟合的多状态模型，PyMSM通过蒙特卡洛模拟生成样本路径。给定协变量、初始状态和时间，它根据估计的模型顺序采样后续状态和每个状态中的持续时间，直到达到终止状态或超过预定义的最大转换次数。在每个观察值采样多个路径后，可以获取摘要统计信息，如状态占用概率和中位状态持续时间。(3) 预定义模型和数据模拟：允许加载或配置预定义的多状态模型，并通过随机路径生成模拟的生存数据，为研究提供了有价值的工具。

6.4 R语言包SmoothHazard
SmoothHazard包[110]旨在拟合疾病-死亡模型中的区间删失数据回归模型。它包括算法，可以同时拟合疾病-死亡模型的三个转换强度的回归模型，其中转换到状态1的时间（见图1c）可能是区间删失的，所有事件时间都可能是右删失的。这三个基线转换强度函数可以通过Weibull分布建模，或者，在半参数框架中通过M-splines建模。给定特定协变量，可以组合估计的转换强度来产生累积发生函数和预期寿命的估计值。

6.5 R语言包pseudo和eventglm
pseudo R包[111]包括计算各种感兴趣的边际参数（见第2.3节）的伪值的功能，例如累积发生函数或状态中的受限平均时间。R语言包eventglm也应用了伪值框架，并包括残差绘制、采样权重的使用和校正方差估计。

6.6 R语言包simMSM
R语言包simMSM[112]模拟多状态模型的事件历史。它能够生成具有潜在非线性基线危险函数的事件历史，以及非线性时间依赖或时间独立协变量的影响，同时考虑了对过去历史的依赖性。通过应用于累积全因危险率的逆向采样来实现事件历史的随机生成。

6.7 R语言包Targeting Illness-Death Models With Within-Subject Random Effects
R语言包SemiCompRisks[113]使用贝叶斯估计技术处理基于脆弱性的疾病-死亡模型，包括第5.1.1节和5.1.2节中描述的条件模型和加性模型。该包提供了结合伽马和正态脆弱性的Cox类型和加速失效时间模型。frailty-LTRC包（可在https://github.com/nirkeret/frailty-LTRC获取）使用第5.1.1节中概述的边际化Cox模型的伪似然方法，采用伽马脆弱性。另一方面，semicompAFT（可在https://github.com/leakats/semicompaft获取）在第5.1.2节讨论的乘法脆弱性设置下实现了半参数AFT模型，也使用伽马脆弱性。frailtypack[114]是一个特别关注处理脆弱性模型的包，适用于复杂场景中的时间至事件数据，包括具有复发事件或竞争风险的多状态模型。

7 多状态建模的优缺点
基于强度的多状态过程建模框架与信息随时间展开和揭示的方式非常吻合。特别是，它认识到过去（历史）影响未来，个体同时面临多种类型事件的风险，且过程可能因多种不同原因终止。通过纳入内部和外部时间变化的协变量，它们可以提供有关动态因素与疾病进展或死亡之间关联的有用见解。基于强度的框架用于检查一种事件的发生如何改变另一种事件的风险，这是局部依赖建模的基础[27]。在2.2.1节中，我们提到了如何在疾病-死亡模型的框架内研究疾病复发与死亡之间的局部依赖性。更一般地，可以形成联合模型来研究疾病的共存（即共病），或复发事件与终止事件之间的关系——Cook等人[115]使用这一框架研究了骨转移患者试验中复发骨骼事件与死亡之间的关系。这可以被视为疾病-死亡过程的扩展，在这种设置中可能发生可数数量的非致命事件。在这种情况下，可以使用Aalen-Johansen估计器[27]来联合估计转换概率矩阵和受限平均事件数量及生存概率。在非动态设置中，两个变量之间的关联通常被认为是对称的，即如果依赖于，则也依赖于。基于强度的模型可以适应不对称关系，其中一个事件的发生可能会改变另一种事件的风险，但第二种事件的风险在第一种事件发生时不会改变。当一个事件是死亡且依赖性必然是不对称的时候，这种做法具有自然吸引力。Aalen[116]指出，局部依赖建模可以提供Granger学派内的因果效应的见解。尽管全面的多状态模型在捕捉复杂事件历史方面具有优势，但它们在实践中并未被广泛采用。这种犹豫可能源于应用研究人员对这些方法的熟悉程度有限，以及担心这些模型需要更复杂的建模假设。这种复杂性可能会引起对稳健性的担忧，特别是当数据稀少或模型假设难以验证时。在癌症临床试验中，感兴趣的时间包括癌症进展时间、无进展死亡时间和癌症进展后的死亡时间。如2.2.1节所示，这些事件时间通常与三状态或四状态疾病-死亡过程一起建模；分别见图1c和d。然而，更常见的是基于无进展生存时间来评估治疗效果，这是一个复合终点，定义为在状态0中度过的时间。虽然避免了竞争风险问题，并允许进行简单的时间至事件分析，但治疗指标的回归系数并不产生具有明确解释的估计量[117]。最近的工作受到ICH-E9附录发布的启发，旨在定义在这些设置中清晰可解释的治疗效果估计量，多状态模型可以在获得不同状态预期时间的估计[118]、不同路径的概率或纳入共干预或治疗转换[119]方面发挥重要作用。然后可以通过分配不同状态的相对值[120]或更明确地指定效用[121]来进行摘要统计的综合。在任何应用中，都必须权衡制定解决潜在过程复杂性的全面模型的愿望与希望获得简单可解释估计量和稳健性的愿望。可能有一个明确的理解，即过程是复杂的，但如果数据有限，则复杂的模型可能难以拟合。这在第4.1节中有所讨论，我们讨论了银屑病关节炎的联合损伤的多状态分析，但个体仅受到间歇性观察。总之，多状态模型提供了几个优势：(i) 它们提供了一个表示复杂过程的全面框架；(ii) 它们能够建模事件之间的依赖性；(iii) 支持对多种类型结果的预测；(iv) 许多来自单变量生存分析的诊断工具可以扩展到这种设置。然而，这些模型也有局限性：(i) 指定正确的历史依赖形式可能很困难；(ii) 相关时间尺度有时不清楚；(iii) 基于强度的建模需要比简单模型更多的假设；(iv) 可靠的估计需要每个转换有足够的事件，这使得样本大小计算比单事件生存分析更为复杂。虽然本文的重点是多状态过程的统计建模和预测，但重要的是要注意，本文讨论的回归模型和方法并不产生具有因果解释的转换动态的估计量。为多状态过程定义因果估计量——特别是在存在时间依赖协变量、中间事件和状态间反馈的情况下——需要指定潜在结果和治疗分配机制的假设。这样的因果估计量通常最自然地用边际属性来表达，例如基于状态占用的状态概率或函数，而不是特定于转换的强度。在这种情况下，基于伪值的方法对于多状态数据的因果分析具有前景。最近的工作已经开始沿着这些线索开发多状态和事件历史的因果框架[119, 122-125]，但多状态模型中的因果推断的全面处理超出了本综述的范围，仍然是正在进行的研究领域。

8 讨论
在本文中，我们介绍了多种多状态过程建模方法，重点介绍了基于强度的模型和边际模型。第2节和第3节阐述了理解这些模型及其实际用途的基础概念。第4节在此基础上进一步讨论了连续时间过程中间歇性观察所带来的挑战，这在临床研究中是一个常见问题，因为状态之间的转换并不连续观察。我们在第5节中对基于脆弱性的模型的探索展示了它们捕捉未观察到的异质性的能力。本文讨论的方法为处理现实世界数据的复杂性提供了关键工具，使研究人员能够更准确地推断过程（如疾病进展）的动态。虽然基于脆弱性的模型为处理特定于受试者的未观察到的协变量提供了有用的框架，但它们也引入了额外的复杂性，特别是在关于依赖结构的假设方面。未来的研究应该专注于将这些方法扩展到更复杂的多状态过程。基于强度的回归模型无法提供对重要协变量缺失具有鲁棒性的模型参数估计。在因果关系的讨论中，以复发状态的占用为条件会导致“碰撞者偏差”（collider bias）；详见Cook和Lawless [16] 的第8.4节。对于时间固定的协变量，评估它们与多状态过程关系的另一种方法是通过分层，并计算每个层中不同状态的预期停留时间。验证多状态模型面临重大挑战，尤其是在使用真实世界数据评估其预测性能或拟合优度时。对现有模型检验方法进行全面的文献回顾，并识别出关键未解决的问题，对于推动这一领域的发展将非常有价值。另一个未在此讨论的重要领域是将机器学习方法应用于多状态生存数据。虽然将机器学习算法应用于多状态生存分析的方法仍在发展中，但它在提高预测准确性方面显示出巨大潜力，尤其是在医疗保健领域 [126-128]。一个关键挑战是在机器学习的预测性能与可解释性和鲁棒性之间找到平衡，这两者对临床决策都至关重要。现代基于生存的机器学习方法还提供了解释工具，包括变量重要性度量和特定状态的风险总结 [129]。另一个主要挑战是针对高级机器学习方法（如深度神经网络）的不确定性量化。最近的一项研究 [130] 提出了一种可以适应深度学习方法的生存分析不确定性量化方法。通过第3.3节讨论的伪值（pseudo-value）方法，可以直接将机器学习算法应用于多状态过程，因为这些算法通常是为横截面数据而非动态数据设计的。当可以指定相关的时间范围时，可以利用机器学习来预测状态占用概率。早期的相关工作集中在竞争风险和半竞争风险过程上 [131]。随着该领域的发展，将机器学习整合到多状态建模框架中可能会为分析复杂的生存数据开辟新的途径。

致谢：
M.G. 的工作得到了以色列科学基金会（项目编号767/21）和特拉维夫大学人工智能与数据科学中心（TAD）的资助。R.J.C. 获得了加拿大自然科学与工程研究委员会（RGPIN-2017-04207）和加拿大卫生研究院（FRN 13887）的发现基金支持。M.A. 是麦吉尔大学生物统计学的杰出詹姆斯·麦吉尔教授。M.P.P. 的工作得到了斯洛文尼亚研究与创新机构（项目编号P3-0154）的资助。作者感谢匿名副编辑和审稿人的宝贵意见和详细建议，这些意见极大地提升了我们手稿的内容质量和可读性。

资金支持：
本项工作得到了以色列科学基金会（项目编号767/21）、特拉维夫大学人工智能与数据科学中心（TAD）、加拿大自然科学与工程研究委员会（项目编号RGPIN-2017-04207）以及加拿大卫生研究院（项目编号FRN 13887）的资助。M.A. 是麦吉尔大学生物统计学的杰出詹姆斯·麦吉尔教授。M.P.P. 的工作得到了斯洛文尼亚研究与创新机构（项目编号P3-0154）的资助。

利益冲突：
作者声明没有利益冲突。

数据可用性声明：
支持本研究结果的数据可在CRAN包中找到：https://github.com/therneau/survival。这些数据来源于以下公共领域的资源：– survival R包，https://github.com/therneau/survival。