简介hesim 是一套面向卫生经济评价与决策分析场景的 R 语言建模工具包适合卫生经济学研究者、药物经济学评价人员及具备一定 R 基础的数据分析师使用用于将统计分析与经济评估整合到统一框架中。资源包共 697 个文件压缩后约 9.72MB包含 60 个 R 源码文件、20 个 C 源文件、102 个 rd 帮助文档、258 个 html 页面及若干 rmd、pdf、rda 等覆盖源码、文档、示例与依赖配置便于直接阅读与二次开发。该包支持队列离散时间状态转换模型、N 状态分区生存模型以及个体级连续时间状态转换模型涵盖 Markov 与半 Markov 过程并借助 Rcpp 与 data.table 实现高效个体级仿真、概率敏感性分析与患者异质性建模。读者可据此完成成本效益分析、参数不确定性传播与决策不确定性表达并参考内置文档与示例快速搭建自己的模拟建模流程。目前已有 471 人学习下载。1. 卫生经济模拟建模到底在算什么从一笔账说起一款新干预方案上市前最常被问到的问题不是它有没有效而是它值不值。这个问题在卫生经济模拟建模和决策分析里被拆成了一套可计算的流程把患者按健康状态分层给每个状态配上成本和效用权重再让一群虚拟患者按转移概率在状态之间流动最后把多年甚至终身的成本和质量调整生命年QALY折算回当下算出一个增量成本效果比ICER。hesim 就是干这件事的 R 包——它把离散事件模拟、状态转移模型Markov / 半马尔可夫、PartSA 这几类主流结构统一到一套数据接口下让你不用为每种模型重写一遍模拟循环。它适合谁做卫生技术评估的从业者、药企市场准入团队、高校里做药物经济学课题的研究生以及需要给报销决策提供量化依据的分析岗。如果你已经会用 Excel 搭 Markov 模型但每次改参数都要手动重算、做概率敏感性分析PSA时只能硬跑几千次hesim 能把你从重复劳动里捞出来。下面按模型怎么建 → 参数怎么定 → 模拟怎么跑 → 坑在哪的顺序讲透。2. 用 hesim 搭一个 Markov 模型从状态定义到转移矩阵2.1 先想清楚三件事状态、周期、视角任何卫生经济模型开工前先把这三个东西钉死否则后面代码写得再漂亮也是白搭。状态是患者可能处于的健康情形比如无进展进展死亡状态必须互斥且穷尽。周期是模型推进的时间步长常见是 1 个月或 1 年周期长度直接决定转移概率怎么换算。视角决定算谁的钱——卫生系统视角只算直接医疗成本社会视角还要算生产力损失和照护成本视角一变结论可能翻转。hesim 的设计哲学是数据驱动你把状态、转移、成本、效用都整理成数据框它负责跑模拟。这跟很多 GUI 软件点鼠标建模型的路子不同好处是可复现、可版本控制、能塞进自动化流程。常见做法是先用data.table把参数表搭好再喂给 hesim 的建模函数。2.2 定义状态与转移把临床路径翻译成矩阵假设一个三状态的简化模型稳定、进展、死亡。死亡是吸收态进去就出不来。转移概率矩阵长这样library(data.table) # 状态顺序1稳定, 2进展, 3死亡 # 每行是一个起始状态每列是下一周期到达的状态 # 这里给的是每周期转移概率周期设为1年 p_stable_to_progress - 0.15 # 稳定-进展 p_stable_to_death - 0.02 # 稳定-死亡 p_progress_to_death - 0.30 # 进展-死亡 # 构造转移概率矩阵行和为1 tmat - matrix(c( 1 - p_stable_to_progress - p_stable_to_death, p_stable_to_progress, p_stable_to_death, 0, 1 - p_progress_to_death, p_progress_to_death, 0, 0, 1 ), nrow 3, byrow TRUE) colnames(tmat) - rownames(tmat) - c(稳定, 进展, 死亡) print(round(tmat, 4))这段代码的关键在于行和必须为 1否则模拟出来的患者会凭空消失或增殖。byrow TRUE保证你按起始状态读行这是最容易搞反的地方——很多人按列填结果模型跑出来人数对不上。参数说明p_stable_to_progress这类概率通常来自临床试验的 Kaplan-Meier 曲线需要先做参数分布拟合再转成周期概率不能直接把 5 年生存率当 1 年概率用。2.3 把转移概率变成 hesim 认的数据结构hesim 的核心是hesim_data对象它把患者、状态、转移、成本、效用打包在一起。下面构造一个最小可跑的例子library(hesim) # 患者数据这里只放1个队列实际可放多个亚组 patients - data.table( patient_id 1:1000, age rnorm(1000, 60, 8) # 年龄作为协变量影响后续概率 ) # 状态数据每个患者在每个状态下的成本与效用 states - data.table( state_id c(1, 2, 3), state_name c(稳定, 进展, 死亡) ) # 转移数据把上面的矩阵拉成长表 transitions - data.table( from c(1, 1, 2), to c(2, 3, 3), prob c(p_stable_to_progress, p_stable_to_death, p_progress_to_death) ) # 成本与效用按状态给死亡状态通常为0 costs - data.table( state_id c(1, 2, 3), cost c(5000, 20000, 0) # 每周期直接医疗成本 ) utilities - data.table( state_id c(1, 2, 3), utility c(0.85, 0.60, 0) # 每周期效用权重 ) hesim_dat - hesim_data( patients patients, states states, transitions transitions, costs costs, utilities utilities )逻辑说明hesim_data()只是把数据装进一个容器它不校验你的概率是否合理所以行和检查要自己做。参数说明cost和utility都是每周期的量如果周期是 1 年那 QALY 就是效用乘以年数如果周期是 1 个月记得把成本和效用都换算成月值否则结果会差 12 倍。这是血泪经验里排名第一的翻车点。3. 参数不确定性怎么塞进去PSA 与分布拟合3.1 为什么点估计不够决策者要的是区间只跑一组参数得到的 ICER 是个点估计决策者会问如果概率高一点、成本低一点结论还成立吗。概率敏感性分析PSA就是让所有不确定参数同时按各自分布抽样跑几千次看 ICER 的分布和成本效果可接受曲线CEAC。hesim 把这套流程做成了create_psa()加run_psa()两步但前提是你得先给每个参数指定分布。常见做法转移概率用 Beta 分布成本用 Gamma 或 Log-normal 分布效用用 Beta 分布。下面演示怎么把前面的固定概率换成分布抽样# 假设 p_stable_to_progress 来自 150 人中 22 人进展 # Beta 分布参数alpha 事件数, beta 非事件数 alpha_prog - 22 beta_prog - 150 - 22 # 抽样 1000 次 set.seed(42) psa_probs - rbeta(1000, alpha_prog, beta_prog) # 检查抽样分布是否合理 summary(psa_probs) quantile(psa_probs, c(0.025, 0.5, 0.975))逻辑说明rbeta()的两个参数直接决定分布形状alpha和beta都大于 1 时分布接近钟形小于 1 时会压向两端。参数说明set.seed()必须设否则每次跑出来的 PSA 结果都不一样审阅者无法复现。抽样次数一般 1000 起步正式分析建议 5000 到 10000但要注意运行时间——hesim 的模拟是向量化的10000 次 PSA 在普通笔记本上通常几分钟内能跑完如果超过半小时先检查是不是把患者数设得过大。3.2 用 hesim 跑 PSA从参数表到结果表把抽样后的参数整理成 hesim 认的格式再调用模拟函数# 构造 PSA 参数数据框每行是一次抽样 psa_params - data.table( sample_id 1:1000, p_stable_to_progress psa_probs, # 其他参数同理这里为简化只放一个 p_stable_to_death rbeta(1000, 8, 142), p_progress_to_death rbeta(1000, 45, 105) ) # 用 hesim 的模拟接口跑状态转移模型 # 这里用简化的自定义循环展示逻辑实际项目用 hesim::run_psa() run_one_sample - function(p1, p2, p3, n_patients 1000, n_cycles 20) { # 初始全部在稳定状态 state - rep(1, n_patients) total_cost - 0 total_qaly - 0 cost_map - c(5000, 20000, 0) util_map - c(0.85, 0.60, 0) for (cycle in 1:n_cycles) { # 按当前状态抽样下一状态 next_state - state idx_stable - which(state 1) idx_prog - which(state 2) # 稳定-进展或死亡 u - runif(length(idx_stable)) next_state[idx_stable] - ifelse(u p1, 2, ifelse(u p1 p2, 3, 1)) # 进展-死亡 u2 - runif(length(idx_prog)) next_state[idx_prog] - ifelse(u2 p3, 3, 2) # 累加本周期成本与效用按周期初状态计 total_cost - total_cost sum(cost_map[state]) total_qaly - total_qaly sum(util_map[state]) state - next_state } c(cost total_cost / n_patients, qaly total_qaly / n_patients) } # 对每个抽样跑一次 results - psa_params[, { r - run_one_sample(p_stable_to_progress, p_stable_to_death, p_progress_to_death) .(cost r[cost], qaly r[qaly]) }, by sample_id] head(results)逻辑说明这段代码用显式循环展示状态转移的抽样逻辑实际生产环境应直接用hesim::run_psa()它内部做了向量化和并行优化速度快一个数量级。参数说明n_cycles 20表示模拟 20 个周期如果周期是 1 年就是 20 年通常要覆盖患者终身或直到 99% 进入死亡状态。n_patients 1000是每个抽样下的虚拟患者数这个数影响的是模拟随机误差不是 PSA 的不确定性一般 1000 到 10000 够用。3.3 算 ICER 和 CEAC把结果翻译成决策语言拿到 PSA 结果后算增量成本和增量 QALY再画成本效果可接受曲线# 假设有两组干预组和对照组这里用模拟数据演示 # 实际中 results 里应有 arm 列区分两组 set.seed(1) results[, arm : rep(c(对照, 干预), each .N/2)] # 按抽样算增量 inc - results[, .( cost mean(cost[arm 干预]) - mean(cost[arm 对照]), qaly mean(qaly[arm 干预]) - mean(qaly[arm 对照]) ), by sample_id] # ICER icer - mean(inc$cost) / mean(inc$qaly) cat(增量成本效果比 ICER , round(icer, 2), \n) # CEAC不同支付意愿阈值下干预组具有成本效果的概率 wtp_seq - seq(0, 100000, by 5000) ceac - sapply(wtp_seq, function(w) { mean(inc$cost - w * inc$qaly 0) }) ceac_dt - data.table(wtp wtp_seq, prob ceac) print(ceac_dt[wtp %in% c(20000, 50000, 100000)])逻辑说明ICER 是增量成本除以增量 QALY但它在增量 QALY 接近 0 时会爆炸所以更稳健的做法是看净货币收益NMB或 CEAC。参数说明wtp_seq是支付意愿阈值不同地区的阈值不同常见参考值是人均 GDP 的 1 到 3 倍但具体取值要按当地指南来不能随便套。ceac算的是在这个阈值下干预组比对照组更划算的概率概率超过某个门槛比如 0.7才认为决策稳健。4. 避坑与排查hesim 建模里最容易翻车的五件事4.1 转移概率行和不为 1患者凭空消失现象模拟跑完后各状态人数加起来不等于初始人数或者死亡状态人数远低于预期。原因转移矩阵按列填了或者某个状态的转移概率漏了一项导致行和小于 1。解决每次构造矩阵后强制检查rowSums(tmat)不等于 1 就报错。更稳妥的做法是用hesim的转移数据长表格式让包内部帮你校验。4.2 周期长度和参数单位不匹配结果差 12 倍现象模型跑出来 QALY 高得离谱或低得离谱成本也对不上。原因转移概率是按年给的但周期设成了月或者成本和效用没跟着换算。解决在代码开头用注释写死本模型周期 1 年所有参数表都按这个单位整理。如果必须用月转移概率要用1 - (1 - p_annual)^(1/12)换算不能直接除以 12。4.3 PSA 抽样次数太少CEAC 曲线锯齿严重现象CEAC 曲线上下跳动看不出稳定趋势。原因抽样次数只有几百次蒙特卡洛误差太大。解决正式分析至少 5000 次跑之前先设set.seed()跑完看 ICER 的 95% 区间是否稳定。如果时间不够先减少患者数比如从 10000 降到 1000因为患者数影响的是模拟精度抽样次数影响的是不确定性刻画后者对决策更重要。4.4 效用值超过 1 或为负QALY 算出来是黑匣子现象QALY 结果异常或者效用权重出现 1.2、-0.1 这种值。原因效用数据来自不同来源量表不统一或者 Beta 分布抽样时参数设错导致越界。解决效用值必须在 [0, 1] 区间死亡状态通常设为 0。用 Beta 分布抽样时如果alpha和beta都很小抽样会压向 0 和 1这时要检查原始数据的均值是否合理。我一般会在抽样后加一道pmax(pmin(u, 1), 0)的截断但这只是后悔药根治方法是核对效用来源。4.5 忽略半周期校正成本和效用被高估现象模型结果比手工 Excel 算的高一点但找不到原因。原因状态转移发生在周期中间但成本和效用按整周期累加了。解决在累加成本和效用时乘以 0.5 做半周期校正或者在周期末统一调整。这个细节在指南里常被提到但实际建模时最容易忘。hesim 本身不强制你做校正得自己在成本效用表里处理。5. 把模型跑快、跑稳、跑可复现的三个进阶习惯第一个习惯是把参数表和分析代码分离。我见过太多项目把概率硬编码在模拟循环里改一个数要翻遍整个脚本。正确做法是建一个params.csv或params.xlsx所有概率、成本、效用都从表里读代码只负责逻辑。这样审阅者能直接看表核对参数来源你也能用data.table::fread()秒级切换场景。第二个习惯是用targets或drake做流程管理。卫生经济模型往往要跑多个亚组、多个场景、多次 PSA手动跑容易漏。把每个步骤定义成 target改一个参数只重跑受影响的部分既省时间又保证可复现。下面是一个最小示例library(targets) list( tar_target(params_file, params.csv, format file), tar_target(params, data.table::fread(params_file)), tar_target(psa_results, run_psa_model(params)), tar_target(ceac_plot, make_ceac(psa_results)) )逻辑说明tar_target()把每个计算步骤声明成依赖节点format file表示这个 target 是文件文件变了才重跑。参数说明run_psa_model()和make_ceac()是你自己写的函数targets 只负责调度不关心内部实现。这套东西上手要半天但一旦用顺后面改模型参数就是改一行 CSV 的事。第三个习惯是给每次分析打版本标签。模型结论会随参数更新而变化如果不同版本的结果混在一起根本说不清哪个是最终版。我的做法是在输出目录里用日期加简短描述命名比如20250115_basecase、20250120_psa_5000同时在脚本开头写一行注释记录本次运行的参数来源和关键假设。这样半年后回头看还能还原当时为什么得出那个 ICER。最后一个具体技巧用hesim::icea()直接出成本效果平面图不要自己手画。它内部处理了增量计算和置信椭圆比手动ggplot2少踩很多坑。跑完 PSA 后调一次把图存成 PDF审阅时直接贴进报告。我自己的习惯是每次模型定稿前先跑一遍 10000 次 PSA看 CEAC 在关键阈值附近是否稳定如果不稳定就加抽样次数或检查参数分布直到曲线平滑为止。这个动作花不了多少时间但能挡掉大部分结论不可靠的质疑。希望帮到你。本文还有配套的精品资源点击获取