药物经济学|不使用VBA的现代Excel建模方法(一)
截至目前，绝大多数用于HTA的药物经济学模型都是使用Microsoft Excel搭建的，然而这种使用电子表格构建模型的方法在近几十年来没有实质性更新，许多经济学模型仍使用Briggs et al.在2006年出版的建模手册中的方法。近年来，得益于R语言强大的向量，矩阵等功能，基于其构建的药物经济学模型被验证具有相比于传统Excel模型执行速度更快且更容易验证等优势。自2020年起，Excel和Google Sheets等表格工具进行了多次更新，增加了对动态数组函数（Dynamic array function）功能的支持，这些函数不仅允许将R语言建模中的许多技术应用在电子表格中，同时也避免了使用复杂的VBA（Visual Basic for Applications）代码进行概率分析和其他复杂计算的需要，大幅简化了电子表格模型的构建并减少了开发、验证和执行所需的时间。本两篇推送将基于发表在PharmacoEconomics的文章《A Modern Approach for Constructing Decision Analytic Models in Microsoft Excel》（在Microsoft Excel中构建决策分析模型的现代方法）一文，讲解如何使用新版Excel的动态数组函数功能构建高效马尔可夫队列模型，包括如何使用单个公式实现从构建马尔可夫轨迹到计算总成本与产出完整的流程，以及如何不使用VBA代码进行抽样、蒙特卡洛模拟及绘制CEAC等。
本篇将以一个四状态（即16种转移可能性）的模型为例，使用动态数组函数构建确定性、时间齐次（Time-homogeneous）的马尔可夫队列模型。下一篇推送将进行概率性马尔可夫队列模型的讲解。首先为治疗方案A和治疗方案B创建转移矩阵，将治疗方案A和B的16个转移概率分别输入到（16*1）的列向量中，并为这个列向量赋予一个范围名称：v_p_A（v代表向量、p代表转移概率、A代表治疗组）。接下来使用转移概率创建4x4的转移矩阵以表示马尔可夫轨迹。在旧版Excel中需要多组公式的计算，但在现代版本的Excel中可以使用新的WRAPROWS函数，通过单个函数从列向量中为治疗方案A创建转移矩阵：输入这个函数后，公式会自动溢出至相邻的单元格以完成矩阵。对于治疗方案B，可以通过将函数内的“A”替换为“B”来完成修改。为每种治疗方案创建转移矩阵后，便可以通过在一个单元格内输入单一公式计算每种治疗方案的马尔可夫轨迹（即每个周期中归属的健康状态）。相较于以往方法，这种计算方式实现了显著效率提升，不仅允许周期数的即时更新，还会根据周期数动态延长或缩短马尔可夫轨迹（注：本节方法仅针对时间齐次，时间依赖模型转移概率的计算需要使用更复杂的函数）。在创建马尔可夫轨迹前，需要为模型设置一个初始状态分布的行向量。在本模型中假设所有患者在初始状态都处于第一个状态中。初始状态分布的1x4向量可以在最左侧的单元格中输入以下静态数组函数定义：为此行初始分布向量的1x4范围赋予一个范围名称：v_init。（若要创建列向量，则将函数中的逗号替换为分号）。接下来通过设定周期长度和时间范围计算模型的总周期数，并赋予这个数字cycles范围名：最后，为此前创建的转移矩阵设定范围名称，治疗方案A和治疗方案B的转移矩阵将分别命名为m_tm_A和m_tm_B（m代表矩阵、tm代表转移矩阵）。命名时应确保名称包括了整个4×4的单元格范围。为转移矩阵赋名后，便可使用新的SCAN和REDUCE函数通过单一公式得到每种治疗方案的马尔可夫轨迹。请注意，SCAN函数在Excel中不可用，因为该函数在Excel中灵活性较差，不支持将向量用作初始值。使用REDUCE函数的替代方法将在下一节给出。尽管在药物经济学中不常用，Google表格的SCAN函数具有更好的灵活性，通过在Google Sheets中输入以下函数可以简单地计算马尔可夫轨迹：= SCAN(v_init, SEQUENCE(cycles + 1), LAMBDA(dist, cycle, IF(cycle > 1, MMULT(dist, m_tm_A), dist)))首先，该函数读取初始状态分布（v_init），然后在SEQUENCE函数中为每个cycle（从1到cycles + 1）迭代执行LAMBDA函数中的计算，并将每次计算的结果存储在临时变量dist中。函数中使用cycles + 1在结果中创建了一个额外的周期，这也允许了后续应用半周期校正或辛普森法校正（Simpson's method）。在第一个周期（cycle = 1）中，LAMBDA函数返回初始状态分布，并在所有后续周期中使用MMULT函数将上一个周期的状态分布（dist）乘以治疗方案A的转移矩阵（m_tm_A），并返回一个四列（v_init行向量的宽度）和cycles + 1行的矩阵，最后将临时变量dist的每次过程计算结果返回在单独一行上。REDUCE函数与SCAN函数的原理相似，但相同结构下其只返回最终计算结果，不会给出所有过程计算结果。由于创建马尔可夫轨迹时需要保留所有过程结果（即每个循环的状态分布），因此SCAN函数更加适用。尽管如此，在Excel中的REDUCE函数相较于SCAN有一个主要优势：支持使用向量作为初始值。因此使用REDUCE函数计算马尔可夫轨迹时需要在SCAN函数式的基础上进行调整，使其在dist的每次计算中保留每个循环的状态分布，以便最终返回完整的马尔可夫轨迹。更新后的函数如下：= REDUCE(v_init, SEQUENCE(cycles), LAMBDA(dist, cycle, VSTACK(dist, MMULT(INDEX(dist, cycle,), m_tm_A))))此函数在每次迭代中执行矩阵乘法（使用MMULT），然后将每个周期的状态分布堆叠在前一个周期的状态分布下方（使用VSTACK）。通过这种方式，在第一个周期后，dist成为了一个矩阵（而非SCAN函数中的行向量），在LAMDA的每次迭代后会在dist的底部添加额外一行，并在后续迭代中读取前一个周期的状态分布（使用INDEX），执行矩阵乘法来计算当前周期的状态分布（使用MMULT），然后将此分布存储在dist的新行中（使用VSTACK）。当SEQUENCE函数中的每个周期都重复此操作后，最终的dist将输入至工作表中。在此函数中，dist和cycle是两个临时变量，更改变量名不会对函数造成影响。以上REDUCE函数在Google Sheets和Excel中的输出结果相同，优势是此方法也适用于Excel。与SCAN函数一样，此函数会创建一个cycle+1周期的马尔可夫轨迹，并允许应用周期内矫正。若不需要矫正，可以将SEQUENCE(cycles)替换为SEQUENCE(cycles - 1)以去掉VSTACK创建出的多余1周期。在继续下一步之前，需要对每个创建的马尔可夫轨迹分配一个合适的范围名。在本案例中将分别命名为m_mt_A和m_mt_B。在使用Excel时，建议为包含动态数组函数的任何单元格（包括马尔可夫轨迹中的左上角单元格）分配一个单独的名称。此名称应与分配给动态数组函数溢出的范围的名称相同，但附加下划线_。在此案例中，每个马尔可夫轨迹中的左上角单元格将分配名称m_mt_A_# 或m_mt_B_# 。这样做的原因是现代版本的Excel支持直接引用动态数组函数溢出的范围，方法是在包含函数的单个单元格后添加“#”。即使用m_mt_A_# 和m_mt_B_# 来引用马尔可夫轨迹。相较于使用m_mt_A或m_mt_B，包含#的引用会在周期数发生变化时动态更新马尔可夫轨迹的长度。为两组治疗方案创建马尔可夫轨迹后，便可以方便地使用单一公式创建按状态和周期细分的成本与产出表。首先为每种治疗方案创建（1*4）的行向量，用于存储每种治疗方案下四个状态中每个状态所产生的每周期总成本，包括直接成本和周期内存在一定时间的间接成本。接下来为每个行向量分配一个v_c_inp的范围名称，并在其后附加一个治疗组标识符（其中c代表成本，inp代表模型输入）。在本模型中，为治疗方案A和B分别命名为v_c_inp_A和v_c_inp_B。随后为健康产出创建一个相同的行向量。若模型中包括了多种健康产出（例如QALY和LY），则为每个产出分别创建行向量。在本模型中健康产出将包括QALY和LY，命名为v_q和v_ly。接下来创建一个列向量（“cycle+1”行*1列）来存储每个周期下的成本与健康产出的周期权重，该向量用于实现折现和周期内矫正。若成本和健康产出使用了不同的折现率，则应分别使用单独的列向量。在此模型中，成本的周期权重向量被命名为v_cw_c，健康产出的向量被命名为v_cw_h。每个向量的第一个单元格分别被命名为v_cw_c_和v_cw_h_，从而允许使用v_cw_c_# 和v_cw_h_# 来引用这些向量。现在便可以通过在成本和产出表中左上角单元格输入一个简短函数，方便的计算每个治疗方案的按状态和周期细分的折现成本和产出，并可应用任何周期内矫正。该函数将马尔可夫轨迹，相应的成本或健康产出向量和周期权重向量相乘，从而创建一个按周期（行）和状态（列）细分的成本和产出矩阵。本模型中成本矩阵的函数为：= m_mt_A_# ∗ v_c_inp_A ∗ v_cw_c_#本函数中使用的# 引用可以确保即使周期数量发生变化，公式仍然有效，并且生成的表格的长度会动态调整。随后可以使用类似的公式创建QALY矩阵，函数式为：= m_mt_A_# ∗ v_q ∗ v_cw_h_#计算LY时仅需将函数式中的 _q 替换为 _ly。新版Excel的BYROW函数可以实现在每个表格中用单一函数计算每个周期，跨所有状态的汇总成本或产出。首先，将治疗方案A的分解成本、QALY和LY的表格分别命名为m_c_A、m_q_A和m_ly_A，并为治疗方案B分配类似的名称（将_A替换为_B）。为计算治疗方案A每个周期、跨所有状态的汇总成本，将以下函数式输入到第一个周期对应的行中的一个新列的单元格中：= BYROW(m_c_A, LAMBDA(cycle, SUM(cycle)))这个函数将溢出并完成整个列向量，避免了在每行中分别输入SUM函数。此函数也可用于QALY和LY，并通过将函数式内的m_c_A替换为所需的范围名称来完成治疗方案B。（若将SUM替换为AVERAGE则可计算平均值，将BYROW替换为BYCOL则可跨数据行进行计算）。通过对分解成本或产出的相应表格求和，可以计算出所有状态和周期的总成本或产出。例如，治疗方案A的总成本可以使用SUM函数计算：尽管不常使用，但计算每个周期的汇总成本和结果（包括当前周期和所有先前周期）可以使决策者能够考虑治疗方案在净健康效益方面达到盈亏平衡所需的时间，从而有效地对模型的时间范围进行敏感性分析。将治疗方案A在每个周期下的汇总成本列向量（使用BYROW函数计算）命名为“v_c_agg_A”。可以通过以下SCAN函数计算每个周期汇总成本的列向量：= SCAN(0, v_c_agg_A, LAMBDA(Φ, cycle, Φ + cycle))对于v_c_agg_A中的每一行，此函数将当前周期的汇总成本（存储在cycle中）与所有先前周期的累积成本（存储在“Φ”中）相加（Excel和Google Sheets均已支持在公式中使用希腊字符和其他符号，“Φ”用于表示“累积”）。SCAN函数中的第一个参数是0因为计算累积成本时所需的起始值为零。与此前用于马尔可夫轨迹的 SCAN 函数不同，此处的SCAN函数在Excel和Google Sheets中都有效，因为初始值不是数组。使用BYROW函数创建汇总列向量“v_c_agg_A”、“v_q_agg_A”和“v_ly_agg_A”后，则可以使用以下LET函数创建包含治疗方案A在每个周期的成本、QALY和LY的三列汇总表格：= LET(c_Φ, SCAN(0, v_c_agg_A, LAMBDA(Φ, cycle, Φ + cycle)), q_Φ, SCAN(0, v_q_agg_A, LAMBDA(Φ, cycle, Φ + cycle)), ly_Φ, SCAN(0, v_ly_agg_A, LAMBDA(Φ, cycle, Φ + cycle)), HSTACK(c_Φ, q_Φ, ly_Φ))表格的最后一行（最后一个cycle）对应于使用SUM函数计算的总成本和结果，因此无需同时给出两个表格，但此表比使用SUM函数得出的总成本和产出表更全面。在本节，我们将基于以上所有函数式，为每种治疗方案创建一个整合函数式，使得输入一个公式便能实现计算马尔可夫轨迹、导出分解成本和健康产出表、并输出一个单行、多列的总成本和总产出表格。这种整合函数式可以输出每个周期的累积成本和产出，其中最后一行对应总和。这种方法避免了需要一个庞大的表格放置马尔可夫轨迹，分解成本和产出等表格，也允许在不使用VBA宏的情况下进行蒙特卡洛模拟，从而减少开发时间并提高概率模型的计算效率（后文）。为构建一个单一的函数式来计算任何治疗方案的总成本和产出，需要该函数执行以下步骤：- 根据转移概率矩阵计算马尔可夫轨迹，以矩阵形式存储在内存中（m_mt）；
- 使用储存的马尔可夫轨迹计算每个状态和周期的分解、折现成本与产出，并根据需要应用周期内校正。这些结果以矩阵形式存储在内存中（m_c、m_q和m_ly）；
- 对这些分解的成本和产出求和，以计算总折现成本和产出。这些以单个值的形式存储在内存中（c、q和ly）；
- 将总折现成本、QALY和LY以具有三列的行向量输出，仅有这个表格会保存到工作表中。
= LET(m_mt, REDUCE(v_init, SEQUENCE(cycles), LAMBDA(dist, cycle, VSTACK(dist, MMULT(INDEX(dist, cycle,), m_tm_A)))), m_c, m_mt * v_c_inp_A_# * v_cw_c_#, m_q, m_mt * v_q_inp_# * v_cw_h_#, m_ly, m_mt * v_ly_inp_# * v_cw_h_#, c, SUM(m_c), q, SUM(m_q), ly, SUM(m_ly), HSTACK(c, q, ly))为了进行向量之间的乘法，需要确保即使周期数量发生变化的情况下，向量也会保持正确的长度。因此在函数中使用“#”引用循环权重向量。m_c、m_q和m_ly的计算中引用马尔可夫轨迹（m_mt）时不需要“#”引用，因为m_mt是在同一个LET函数中计算的。此公式也可以通过整合SCAN函数，以输出每个周期的累积成本和产出。整合后的公式为：= LET(m_mt, REDUCE(v_init, SEQUENCE(cycles), LAMBDA(dist, cycle, VSTACK(dist, MMULT(INDEX(dist, cycle,), m_tm_A)))), m_c, m_mt * v_c_inp_A_# * v_cw_c_#, m_q, m_mt * v_q_inp_# * v_cw_h_#, m_ly, m_mt * v_ly_inp_# * v_cw_h_#, v_c_agg, BYROW(m_c, SUM), v_q_agg, BYROW(m_q, SUM), v_ly_agg, BYROW(m_ly, SUM), c_Φ, SCAN(0, v_c_agg, LAMBDA(Φ, cycle, Φ + cycle)), q_Φ, SCAN(0, v_q_agg, LAMBDA(Φ, cycle, Φ + cycle)), ly_Φ, SCAN(0, v_ly_agg, LAMBDA(Φ, cycle, Φ + cycle)), HSTACK(c_Φ, q_Φ, ly_Φ))若后续希望对时间范围进行敏感性分析，则应优先选择此公式。现在可以使用每种治疗方案的总成本和产出计算ICER。首先，为治疗方案A的总成本、QALY和LY赋名：c_A、q_A和ly_A，并为治疗方案B分配类似的名称。增量成本（c_inc）可以通过c_B减c_A计算得出，QALY和LY也以类似的方式计算，随后可计算ICER：生命年（LY）的增量成本可以通过将q_inc替换为ly_inc计算得出。最后，可以使用单个LET公式创建一个汇总结果表，其中包括每种治疗的总成本和结果、增量成本和结果以及ICER：= LET(c_inc, c_B - c_A, q_inc, q_B - q_A, ly_inc, ly_B - ly_A, icer_q, c_inc / q_inc, icer_ly, c_inc / ly_inc, tx_A, HSTACK(c_A, "N/A", q_A, "N/A", ly_A, "N/A", "N/A", "N/A"), tx_B, HSTACK(c_A, c_inc, q_A, q_inc, ly_A, ly_inc, icer_q, icer_ly), VSTACK(tx_A, tx_B))
下篇推送将基于本篇的确定性模型构建概率性马尔可夫队列模型，包括不使用VBA进行随机抽样、蒙特卡洛模拟，及绘制CEAC等功能的实现。
