简介hesim 是一套面向卫生经济评价与决策分析场景的 R 语言建模工具包适合卫生经济学研究者、药物经济学评价人员及具备一定 R 基础的数据分析者使用。它提供将统计分析与经济评估相集成的通用框架支持队列离散时间状态转换模型、N 状态分区生存模型以及个体级连续时间状态转换模型覆盖 Markov 与半 Markov 过程并借助 Rcpp 与 data.table 实现高效仿真、概率敏感性分析与患者异质性建模。资源包共 697 个文件约 9.72MB包含 60 个 R 源码、20 个 C 实现、102 个 rd 帮助文档、258 个 html 页面及若干 rmd 教程与 pdf 说明便于查阅函数用法与模型示例。目前已有 471 人学习下载。读者可据此掌握从参数分布或自举采样、PSA 不确定性传播到成本效益分析与决策不确定性呈现的完整流程并参考源码理解模型构建与计算加速思路。1. 从一张成本效果散点图说起hesim 到底解决什么问题如果你做过卫生经济学评价大概率经历过这样的场景模型跑完ICER 算出来是每 QALY 三万块决策者问「这个结论稳不稳」你只能硬着头皮说「我们做了单因素敏感性分析」。问题在于单因素分析一次只动一个参数而真实世界里成本、效用、转移概率是同时波动的。hesim 这个 R 包就是冲着这个痛点来的——它把微观模拟、离散事件模拟和分区生存模型统一到一套框架里让你能跑概率敏感性分析、算成本效果可接受曲线还能做个体层面的异质性建模。这个包适合两类人一类是做卫生技术评估的从业者手里有 Markova 或 TreeAge 的模型想迁移到开源环境另一类是流行病学或公共卫生背景的研究者需要把干预措施的成本效果算清楚但不想被商业软件锁死。它的核心价值在于「概率化」——不是给你一个点估计而是给你一个分布。你拿到的不只是「这个方案更划算」而是「在 95% 的模拟里这个方案有 78% 的概率更划算」。hesim 的底层逻辑是「先定义模型结构再生成模拟数据最后做决策分析」。它不要求你一开始就写对全部参数而是允许你用define_model()搭骨架用define_parameters()填血肉最后用run_model()跑出结果。这种分层设计的好处是你可以在参数层反复调整而不动模型结构也可以在模型层切换微观模拟和队列模拟而不重写参数。接下来我会按「装环境 → 搭模型 → 跑模拟 → 排错 → 进阶」的顺序把每个环节的代码和参数含义拆开讲。2. 环境准备与核心对象把 hesim 装进工作流2.1 安装与依赖检查hesim 在 CRAN 上有稳定版本但如果你要用最新的微观模拟功能建议从 GitHub 装开发版。我一般会先检查 R 版本低于 4.0 的话data.table和R6的某些特性会报错。# 检查 R 版本低于 4.0 建议先升级 getRversion() # 从 CRAN 安装稳定版 install.packages(hesim) # 如果需要开发版用 devtools 从 GitHub 装 # install.packages(devtools) # devtools::install_github(hesim-dev/hesim) # 加载核心包 library(hesim) library(data.table)装完之后别急着跑模型先确认几个关键依赖data.table用于高效处理模拟数据R6用于面向对象建模flexsurv用于参数生存模型拟合。如果flexsurv没装后面做生存分析会直接卡住。# 检查关键依赖是否可用 pkgs - c(data.table, R6, flexsurv, ggplot2) sapply(pkgs, requireNamespace, quietly TRUE)返回全是TRUE才算环境干净。有一个FALSE就去补装别跳过。2.2 理解 hesim 的三个核心对象hesim 的建模流程围绕三个对象展开hesim_data管数据params管参数model管结构。很多人第一次用觉得绕是因为把这三个混在一起写了。我习惯先画一张表把职责分清对象类型创建函数职责常见输入格式数据对象hesim_data()存储患者特征、策略、时间点data.table / data.frame参数对象params_surv()/params_lm()存储回归系数或分布参数列表或矩阵模型对象define_model()定义转移结构、干预策略公式或矩阵hesim_data()是最先要建的。它要求你提供三个要素strategies干预策略名称、patients患者层面的协变量、states健康状态。这三个要素的维度必须对齐否则后面run_model()会报「维度不匹配」——这是新手最常见的翻车点。# 定义三个干预策略 strategies - data.table( strategy_id 1:3, strategy_name c(标准治疗, 新药A, 新药B) ) # 定义患者特征这里模拟 1000 个患者 set.seed(42) patients - data.table( patient_id 1:1000, age rnorm(1000, mean 65, sd 10), sex rbinom(1000, 1, 0.5), risk_score runif(1000, 0, 1) ) # 定义健康状态 states - data.table( state_id 1:3, state_name c(稳定, 进展, 死亡) ) # 组装成 hesim_data 对象 hesim_dat - hesim_data( strategies strategies, patients patients, states states ) # 打印结构确认 print(hesim_dat)这段代码的逻辑是先把三个维度的数据分别建好再用hesim_data()打包。set.seed(42)是为了让患者特征可复现你做概率敏感性分析时如果每次患者样本都变结果就没法对比。age用正态分布、sex用二项分布、risk_score用均匀分布这是卫生经济评价里最常见的协变量生成方式。参数说明strategies的strategy_id必须是整数且从 1 开始patients的patient_id同理states的state_id也要连续。这三个 ID 是 hesim 内部索引的锚点断了就会报错。3. 模型搭建从转移概率到微观模拟3.1 用define_model()定义分区生存结构分区生存模型的核心是转移矩阵。hesim 允许你用两种方式定义一种是直接给转移概率矩阵另一种是用公式指定每个转移的回归关系。我一般先用矩阵方式跑通流程再换成公式方式做异质性分析。# 定义三个状态之间的转移矩阵 # 行是起点状态列是终点状态 # 这里假设不能从死亡回到其他状态 transition_matrix - matrix( c(0, 0.3, 0.1, # 从稳定出发留在稳定0去进展0.3去死亡0.1 0, 0, 0.2, # 从进展出发去稳定0留在进展0去死亡0.2 0, 0, 1), # 从死亡出发只能留在死亡 nrow 3, byrow TRUE ) # 用 define_model 定义模型结构 mod - define_model( transition_matrix transition_matrix, n_cycles 20, # 模拟 20 个周期 cycle_length 1, # 每个周期 1 年 starting_state 1 # 所有患者从稳定状态开始 ) # 查看模型摘要 summary(mod)transition_matrix的每一行代表一个起点状态每一列代表终点状态。第 1 行第 2 列的 0.3 表示从稳定到进展的概率是 0.3。注意第 3 行第 3 列的 1表示死亡状态是吸收态进去就出不来。n_cycles 20意味着模拟 20 年cycle_length 1表示每年一个周期。如果你的研究时限是 10 年但每半年一个周期就把n_cycles设成 20、cycle_length设成 0.5。提示转移矩阵的每一行概率之和必须等于 1否则run_model()会报「概率不合法」。我见过有人把 0.3 和 0.1 填进去忘了补留在原状态的概率结果跑出来患者凭空消失。3.2 用define_parameters()绑定成本和效用模型结构搭好后下一步是把成本和效用参数绑上去。hesim 的参数对象支持三种形式固定值、分布、回归方程。做概率敏感性分析时你至少要把关键参数设成分布。# 定义成本参数每个状态每年的成本 cost_params - params_lm( formula cost ~ 1, data data.frame( state_id c(1, 2, 3), cost c(5000, 15000, 0) # 稳定5000进展15000死亡0 ) ) # 定义效用参数每个状态的健康效用值 utility_params - params_lm( formula utility ~ 1, data data.frame( state_id c(1, 2, 3), utility c(0.8, 0.5, 0) # 稳定0.8进展0.5死亡0 ) ) # 把参数绑定到模型 mod - set_parameters( mod, cost cost_params, utility utility_params )params_lm()是线性模型参数容器这里虽然只填了固定值但它的结构支持后续替换成回归系数。formula cost ~ 1表示截距模型data里的state_id对应健康状态。成本参数里死亡状态设为 0 是惯例因为死者不产生医疗费用。效用参数里死亡状态也设为 0这是 QALY 计算的标准做法。如果你要做概率敏感性分析把固定值换成分布# 用正态分布替换固定值 cost_params_psa - params_lm( formula cost ~ 1, data data.frame( state_id c(1, 2, 3), mean c(5000, 15000, 0), sd c(500, 1500, 0) # 标准差设为均值的 10% ), uncertainty normal # 指定分布类型 )uncertainty normal告诉 hesim 在每次模拟中从正态分布抽样。sd设成均值的 10% 是经验做法如果你有文献来源的置信区间用(上限 - 下限) / 3.92反推标准差更严谨。3.3 跑微观模拟与结果提取参数绑好后就可以跑模型了。hesim 支持两种模拟方式run_model()跑队列模拟run_microsim()跑个体模拟。做异质性分析必须用后者。# 跑微观模拟1000 个患者100 次概率抽样 results - run_microsim( mod, n_patients 1000, n_sims 100, seed 123 ) # 提取成本效果结果 ce_results - summarize_ce( results, type icer, comparator 标准治疗 ) # 查看前几行 head(ce_results)n_patients 1000控制个体模拟的样本量n_sims 100控制概率敏感性分析的抽样次数。这两个参数直接决定运行时间——1000 患者乘 100 次模拟等于 10 万次个体轨迹普通笔记本大概跑 2 到 5 分钟。seed 123保证结果可复现。summarize_ce()的comparator参数指定基线策略ICER 是相对于它算的。跑完之后你会拿到一个包含strategy_id、cost、qaly、icer的 data.table。如果icer列出现NaN说明该策略在所有模拟中都不比基线更有效这是正常的不是 bug。4. 避坑与排查那些让我重跑一整天的问题4.1 转移概率矩阵行和不为 1现象run_model()报错Error: transition probabilities must sum to 1。原因手动填矩阵时漏了留在原状态的概率。比如从稳定出发你填了去进展 0.3、去死亡 0.1但忘了留在稳定应该是 0.6。解决用rowSums()检查每一行或者用normalize()函数自动归一化。# 检查行和 rowSums(transition_matrix) # 自动归一化慎用会掩盖逻辑错误 transition_matrix_norm - transition_matrix / rowSums(transition_matrix)我一般不建议自动归一化因为它会把你的逻辑错误悄悄修掉。手动检查每一行确认概率来源有依据。4.2 患者 ID 与策略 ID 不连续现象hesim_data()创建成功但run_microsim()报Error: patient_id must be consecutive integers。原因你从外部 CSV 读入患者数据时ID 可能跳号或者从 0 开始。解决用setorder()排序后重新生成连续 ID。# 重新生成连续 ID patients[, patient_id : .I] setorder(patients, patient_id).I是 data.table 的内置变量返回当前行的序号。这招对策略 ID 和状态 ID 同样适用。4.3 效用值超过 1 或为负现象QALY 算出来大于生命周期年限或者出现负值。原因效用参数没有做边界约束正态分布抽样时抽到了大于 1 或小于 0 的值。解决在params_lm()里加lower和upper约束或者用 Beta 分布替代正态分布。# 用 Beta 分布约束效用值在 0 到 1 之间 utility_params_beta - params_lm( formula utility ~ 1, data data.frame( state_id c(1, 2, 3), mean c(0.8, 0.5, 0), sd c(0.05, 0.08, 0) ), uncertainty beta )Beta 分布的参数化方式与正态不同hesim 内部会用矩估计法转换。如果你的效用均值接近 0 或 1Beta 分布比正态更合适。4.4 模拟次数太少导致 ICER 不稳定现象每次跑run_microsim()得到的 ICER 差异很大置信区间宽得没法用。原因n_sims设得太小概率空间没覆盖全。解决先跑 10 次看 ICER 的变异系数如果超过 20%加到 100 次以上。# 快速检查 ICER 稳定性 test_results - run_microsim(mod, n_patients 500, n_sims 10, seed 1) test_ce - summarize_ce(test_results, type icer, comparator 标准治疗) sd(test_ce$icer) / mean(test_ce$icer) # 变异系数变异系数低于 0.1 说明模拟次数够了高于 0.2 就继续加。别一上来就设 1000 次先跑小样本探路。4.5 生存模型外推超出数据范围现象用flexsurv拟合生存曲线后外推到 20 年时生存概率变成负数。原因参数生存模型在尾部可能不稳定尤其是 Gompertz 或 Weibull 分布。解决限制外推年限或者换用限制性立方样条。# 检查外推边界 fit - flexsurvreg(Surv(time, status) ~ 1, data survival_data, dist weibull) summary(fit, t c(1, 5, 10, 20)) # 看 20 年时的生存概率如果 20 年概率为负把模型时限缩到 15 年或者在define_model()里加max_cycles限制。5. 进阶技巧用概率敏感性分析做决策曲线5.1 算成本效果可接受曲线概率敏感性分析的最终产出不是 ICER 点估计而是成本效果可接受曲线CEAC。它回答的是「在 willingness-to-pay 阈值下每个策略成为最优的概率是多少」。# 定义支付意愿阈值范围 wtp_range - seq(0, 100000, by 5000) # 算 CEAC ceac - calc_ceac( results, wtp wtp_range, comparator 标准治疗 ) # 画图 library(ggplot2) ggplot(ceac, aes(x wtp, y prob, color strategy_name)) geom_line(size 1) labs(x 支付意愿阈值 (元/QALY), y 成为最优策略的概率) theme_minimal()calc_ceac()的内部逻辑是对每个 WTP 阈值算每个策略的净货币收益NMB WTP × QALY - Cost然后统计每个策略在 100 次模拟中 NMB 最高的频率。wtp_range从 0 到 10 万步长 5000覆盖了大多数中高收入地区的阈值范围。如果你的研究针对特定地区把上限调到当地人均 GDP 的 3 倍即可。5.2 用情景分析替代单因素敏感性分析单因素敏感性分析一次只动一个参数但决策者往往关心「如果新药价格降 20% 且效用提高 0.05 同时发生」会怎样。hesim 的情景分析可以同时改多个参数。# 定义情景新药价格降 20%效用提高 0.05 scenario_params - params_lm( formula cost ~ 1, data data.frame( state_id c(1, 2, 3), cost c(5000, 12000, 0) # 进展状态成本从 15000 降到 12000 ) ) # 把情景参数绑到模型 mod_scenario - set_parameters(mod, cost scenario_params) # 跑情景模拟 scenario_results - run_microsim(mod_scenario, n_patients 1000, n_sims 100, seed 123) scenario_ce - summarize_ce(scenario_results, type icer, comparator 标准治疗)对比ce_results和scenario_ce的 ICER 变化你就能量化「降价 提效」的联合影响。这比单因素分析更接近真实决策场景。5.3 一个我踩过的坑模拟种子与并行计算最后说一个血泪经验。hesim 支持并行计算加速但如果你在run_microsim()里设了seed又开了并行不同核心上的随机数序列会冲突导致结果不可复现。我当时的解决方案是要么关并行要么用future包显式管理随机数种子。# 方案一关并行老老实实跑 results - run_microsim(mod, n_patients 1000, n_sims 100, seed 123) # 方案二用 future 管理并行种子 library(future) plan(multisession, workers 4) results - run_microsim(mod, n_patients 1000, n_sims 100, seed 123, parallel TRUE) plan(sequential) # 跑完恢复默认从那以后我每次跑概率敏感性分析都强制走一遍「小样本探路 → 检查变异系数 → 确认种子 → 再放大样本」的流程再也没出现过结果对不上的情况。希望帮到你。本文还有配套的精品资源点击获取