PyHEOR：用Python做卫生经济学建模

,

为什么写 PyHEOR

做卫生经济学评价（Health Economics and Outcome Research, HEOR）的同学大多用 R（heemod、hesim、DARTH）或者 TreeAge。R 生态确实成熟，但我在实际项目中遇到了一些痛点：

  • R 的多个包 API 风格不统一，拼一个完整流程需要在不同包之间跳来跳去
  • TreeAge 是商业软件，价格不菲，而且不够灵活
  • Python 在数据处理、机器学习方面生态更好，但缺少一个完整的 HEOR 框架

所以我写了 PyHEOR——一个纯 Python 的卫生经济学建模框架，目标是用一套统一的 API 覆盖 HEOR 的全流程。

项目地址：https://github.com/lenardar/PyHEOR

能做什么

PyHEOR 目前支持以下功能：

建模引擎

  • Markov 队列模型（离散时间状态转移）
  • 分区生存模型（PSM）
  • 微观模拟（个体水平状态转移）
  • 离散事件模拟（DES，连续时间）

证据合成

  • IPD 生存曲线拟合（6 种分布，AIC/BIC 比较）
  • KM 曲线数字化重建（Guyot method，从文献 KM 图反推 IPD）
  • NMA 后验样本整合（导入 R 包产生的 MCMC 样本）

分析与决策

  • 基线分析、单因素敏感性分析（OWSA）、概率敏感性分析（PSA）
  • 多策略比较：效率前沿、ICER、NMB、CEAC、CEAF、EVPI
  • 预算影响分析（BIA）
  • 模型校准（Nelder-Mead / 随机搜索）

导出与可视化

  • 28 种专业图表
  • Excel 多 Sheet 导出 + Excel 公式验证模型（用于交叉验证 Python 结果）
  • Markdown 一键报告（generate_report()，自动运行全部分析并生成报告 + 配套图片）

设计思路

一个模型对象搞定所有分析

定义好模型后，run_base_case()、run_owsa()、run_psa() 一气呵成，不需要为不同分析类型重新组织代码。同一个模型对象，确定性分析、敏感性分析、概率分析全部搞定：

1
2
3
4
5
6
result = model.run_base_case()   # 确定性基线
owsa   = model.run_owsa()        # 单因素敏感性
psa    = model.run_psa(n_sim=1000)  # 概率敏感性

# 一键生成完整 Markdown 报告（含龙卷风图、CE 平面、CEAC 等）
ph.generate_report(model, "report.md")

在 R 里做同样的事情，往往需要手动把参数拆开重组，或者在不同的包之间传递数据。PyHEOR 把这些都封装在模型内部了。

ph.C 补数哨兵

写转移矩阵时最烦的就是对角线元素——要手动算 1 - sum(其余)，参数一多就容易出错。PyHEOR 用一个 ph.C 占位符自动补齐：

1
2
3
4
5
model.set_transitions("SOC", lambda p, t: [
    [ph.C,  p["p_HS"], p["p_HD"]],  # ph.C = 1 - p_HS - p_HD
    [0,     ph.C,      p["p_SD"]],  # ph.C = 1 - p_SD
    [0,     0,         1        ],
])

这个设计参考了 R heemod 的 C 常量，用过的同学应该很熟悉。

Lambda 定义一切

费用、效用、转移概率都用 lambda 函数定义。好处是天然支持时变逻辑和参数依赖，不需要额外的配置项：

1
2
3
4
5
6
7
8
9
# 前 24 个月用药费 5000，之后降为 2000
model.set_state_cost("Sick", "Trt", lambda p, t: 5000 if t < 24 else 2000)

# 年龄依赖的死亡概率（微观模拟中可以访问个体属性）
model.set_transitions("SOC", lambda p, t, attrs: [
    [ph.C,  p["p_HS"] * (1 + (attrs["age"] - 55) * 0.02), 0.01],
    [0,     ph.C,  p["p_SD"]],
    [0,     0,     1],
])

注意微观模拟的 lambda 有三个参数 (p, t, attrs)，第三个参数是患者属性。引擎会根据 lambda 的参数个数自动判断是队列模式还是个体模式。

统一的参数系统

一个 add_param() 调用，同时定义基线值、OWSA 范围和 PSA 分布。不同分析自动取用对应的值：

1
2
3
4
5
6
model.add_param("p_HS",
    base=0.15,                        # 确定性分析用
    low=0.10, high=0.20,              # OWSA 范围
    dist=ph.Beta(mean=0.15, sd=0.03), # PSA 抽样分布
    label="健康→生病概率",             # 图表显示名
)

内置分布：Beta、Gamma、Normal、LogNormal、Uniform、Triangular、Dirichlet、Fixed。参数化方式用 mean/sd，不需要手动算 alpha/beta。

贴现率也纳入了参数系统——传入 Param 对象即可自动参与 OWSA 和 PSA，无需额外 add_param()：

1
2
3
4
5
model = ph.MarkovModel(
    ...,
    dr_cost=ph.Param(0.03, low=0.0, high=0.08, label="费用贴现率"),
    dr_qaly=ph.Param(0.03, low=0.0, high=0.05, label="效用贴现率"),
)

灵活的费用体系

除了基础的状态费用，还支持多种特殊费用场景：

1
2
3
4
5
6
7
8
9
10
11
# 首周期一次性费用（如入组筛查费）
model.set_state_cost("PFS", "Trt", lambda p, t: 50000, first_cycle_only=True)

# 限定应用周期（如 24 个月停药）
model.set_state_cost("PFS", "Trt", lambda p, t: p["c_drug"], apply_cycles=(0, 24))

# 转移费用：从 PFS 进展到 Progressed 时的手术费
model.set_transition_cost("surgery", "PFS", "Progressed", 50000)

# 转移费用计划表：手术 + 后续两周期随访
model.set_transition_cost("surgery", "PFS", "Progressed", [50000, 10000, 10000])

转移费用是基于每周期转移流量自动计算的，不受半周期校正影响。计划表通过卷积处理多批次转入患者的费用叠加。

Excel 公式验证

做 HEOR 审稿人经常要求提供 Excel 验证模型。PyHEOR 可以导出一个用 Excel 公式独立计算的完整模型文件：

1
ph.export_excel_model(result, "verification.xlsx")

导出的 Excel 里，转移矩阵是输入区（黄色底色），Trace、费用、QALY 全部用 SUMPRODUCT 等公式计算，Summary sheet 显示 Excel 结果 vs Python 结果的差异（应该是 ~0）。这样审稿人可以在 Excel 里直接验证模型逻辑。

示例：Markov 队列模型

用 Markov 模型做一个三状态（健康→生病→死亡）的成本效果分析：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import pyheor as ph

model = ph.MarkovModel(
    states=["Healthy", "Sick", "Dead"],
    strategies=["Standard", "New Treatment"],
    n_cycles=40,
    cycle_length=1,
    dr_cost=ph.Param(0.03, low=0.0, high=0.08, label="费用贴现率"),
    dr_qaly=ph.Param(0.03, low=0.0, high=0.05, label="效用贴现率"),
    half_cycle_correction=True,
)

# 添加参数
model.add_param("p_HS", base=0.15, dist=ph.Beta(mean=0.15, sd=0.03))
model.add_param("p_SD", base=0.10, dist=ph.Beta(mean=0.10, sd=0.02))

# 转移矩阵
model.set_transitions("Standard", lambda p, t: [
    [ph.C,  p["p_HS"], 0.02],
    [0,     ph.C,      p["p_SD"]],
    [0,     0,         1],
])

model.set_transitions("New Treatment", lambda p, t: [
    [ph.C,  p["p_HS"] * 0.7, 0.02],
    [0,     ph.C,             p["p_SD"] * 0.8],
    [0,     0,                1],
])

# 费用和效用
model.set_state_cost("medical", {"Healthy": 500, "Sick": 3000, "Dead": 0})
model.set_state_cost("drug", {
    "Standard": {"Healthy": 0, "Sick": 0, "Dead": 0},
    "New Treatment": {"Healthy": 2000, "Sick": 2000, "Dead": 0},
})
model.set_utility({"Healthy": 0.95, "Sick": 0.60, "Dead": 0.0})

# 运行
result = model.run_base_case()
print(result.summary())
print(result.icer())

# OWSA 龙卷风图（贴现率通过 Param 自动参与）
owsa = model.run_owsa()
owsa.plot_tornado()

# PSA
psa = model.run_psa(n_sim=1000)
psa.plot_scatter(wtp=50000)
psa.plot_ceac(wtp_range=(0, 150000))

# 一键报告
ph.generate_report(model, "report.md")

从定义模型到出结果，代码量很少，逻辑也比较清晰。

示例：分区生存模型 (PSM)

肿瘤经济学评价最常用的就是 PSM。模型基于两条参数化生存曲线（OS 和 PFS）划分三个状态的占比，不需要写转移矩阵：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
import pyheor as ph

psm = ph.PSMModel(
    states=["PFS", "Progressed", "Dead"],
    survival_endpoints=["PFS", "OS"],
    strategies=["SOC", "New Drug"],
    n_cycles=120,
    cycle_length=1/12,  # 月周期
    dr_cost=0.03,
    dr_qaly=0.03,
)

# 基线生存曲线
baseline_pfs = ph.LogLogistic(shape=1.5, scale=18)
baseline_os = ph.Weibull(shape=1.2, scale=36)

# SOC：直接用基线
psm.set_survival("SOC", "PFS", baseline_pfs)
psm.set_survival("SOC", "OS", baseline_os)

# New Drug：在 SOC 基础上应用 HR 和加速因子
psm.set_survival("New Drug", "OS",
    lambda p: ph.ProportionalHazards(baseline_os, hr=0.7))
psm.set_survival("New Drug", "PFS",
    lambda p: ph.AcceleratedFailureTime(baseline_pfs, af=1.3))

# 费用
psm.set_state_cost("treatment", {
    "SOC": {"PFS": 1000, "Progressed": 2500, "Dead": 0},
    "New Drug": {"PFS": 6000, "Progressed": 2500, "Dead": 0},
})

# 效用
psm.set_utility({"PFS": 0.80, "Progressed": 0.55, "Dead": 0.0})

# 运行
result = psm.run_base_case()
print(result.summary())
print(result.icer())

# 可视化
result.plot_survival()    # 生存曲线
result.plot_state_area()  # 状态面积图

PSM 的核心是生存分布的选择。PyHEOR 内置了 10 种参数化分布（Exponential、Weibull、LogLogistic、LogNormal、Gompertz、Generalized Gamma 等），还支持 ProportionalHazards 和 AcceleratedFailureTime 包装器，可以很方便地在基线曲线上应用 HR 或 AF。

从文献 KM 图到建模

做 HEOR 经常遇到一个问题：文献只给了 KM 曲线图，没有原始数据。PyHEOR 集成了 Guyot method（Guyot et al. 2012），可以从数字化的 KM 坐标反推 IPD，然后直接拟合参数分布用于建模：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 1. 从 WebPlotDigitizer 等工具获取 KM 坐标
t_digitized = [0, 2, 4, 6, 8, 10, 12, 15, 18, 21, 24]
s_digitized = [1.0, 0.92, 0.83, 0.74, 0.66, 0.58, 0.50, 0.40, 0.32, 0.25, 0.20]

# 2. 从文献中读取 number-at-risk 表
t_risk = [0, 6, 12, 18, 24]
n_risk = [120, 88, 60, 38, 22]

# 3. 重建 IPD（内置坐标预处理：去噪、排序、单调化）
ipd_time, ipd_event = ph.guyot_reconstruct(
    t_digitized, s_digitized,
    t_risk, n_risk,
    tot_events=96,  # 可选：文献报告的总事件数
)

# 4. 拟合 6 种分布，AIC/BIC 自动比较
fitter = ph.SurvivalFitter(ipd_time, ipd_event, label="OS")
fitter.fit()
print(fitter.summary())        # AIC/BIC 比较表
best = fitter.best_model()     # 自动选最优

# 5. 直接用于 PSM 建模
psm.set_survival("SOC", "OS", best.distribution)

手动数字化的坐标难免有噪声（手抖、非单调点），guyot_reconstruct 内部会自动调用预处理（参考了 R 包 IPDfromKM 的策略），包括异常值检测、重复时间处理、强制单调非递增等。

这个「文献 KM 图 → IPD → 参数拟合 → 建模」的流程在实际项目中非常实用。

后续计划

PyHEOR 目前的功能基本覆盖了 HEOR 常见需求。后续可能会做的事情：

  • 更多的示例和教程
  • 与 R 生态的结果对比验证
  • 性能优化（特别是微观模拟和 DES）

另外，未来将考虑增强 PyHEOR 与 AI 的适配性，包括：

  • 结构化输出（to_dict / to_json），让分析结果直接可被 LLM 读取和理解
  • 自动解读（interpret(wtp)），一键生成标准化的结论文本
  • 自然语言建模接口，通过 JSON Schema 定义模型，LLM 无需编写 Python 代码即可完成建模

如果你也在做卫生经济学评价，欢迎试用和反馈：PyHEOR on GitHub