NIPT 的时点选择与胎儿的异常判定（优秀范文三：蒙特卡洛鲁棒时点与阈值优化）

国赛 2025 C 题 · 优秀范文三。视角：把范文一的"确定性参数假设"打破——患病率先验、区分度斜率、假阴/假阳代价比、羊穿风险都不是精确已知的，而是存在真实波动。本文用蒙特卡洛（MC）仿真把这四类参数不确定性逐点传播到"最优检测孕周 $g^*$ 、风险报告 cut-off $z^*$ 、期望代价 $E^*$ "上，给出一整套稳健区间而非单点值，并与范文一的确定性解做对比。正文 / 配图 / 附录 / 真源四路数字一致（固定种子 2026）。

摘要

NIPT 筛查方案设计（何时测、用什么阈值报风险）依赖患病率先验、检测区分度、假阴与假阳的相对代价、羊穿风险等参数，而真实场景里这些参数并非确定值：人群患病率有波动、区分度随样本质量变化、临床对"漏诊 vs 不必要羊穿"的相对重视程度也因人而异。若只用单一标称参数求一个"最佳时点/阈值"，会在参数偏离时付出意料之外的代价。本文用蒙特卡洛仿真对 4 类关键参数做概率建模（对数正态扰动），重复求解最优设计 4000 次，得到：最优孕周 $g^*$ 的均值为 13.77 周、标准差仅 0.25、95% 区间 $[13.5,14.0]$ ；最优 cut-off $z^*$ 均值 2.20、95% 区间 $[1.93,2.47]$ ；期望代价 $E^*$ 均值 4.569、95% 区间 $[3.767,5.568]$ 。与确定性解（ $g^*=14.0,z^*=2.16$ ）相比，MC 给出的最优时点略提前、cut-off 略保守，且发现时点对"区分度斜率"最敏感（单降 25% 即把最优时点推前约 1 周），而羊穿风险的影响最小。据此建议：方案落地应取"中等保守的鲁棒解"并配套监测参数漂移，而非依赖名义点估计。

一、问题重述

题目核心是回答"什么时候做 NIPT、报告阈值怎么定"。范文一给出了确定性框架：假设一组参数（年龄别患病率、胎儿浓度随孕周上升的曲线、区分度、代价权重），求使每筛查一人期望负担最小的 $(g^*,z^*)$ ；范文二把重点放到"多指标判定"的分类侧。但现实中这些参数没有精确值：不同人群的患病率、实验室的测序质量、医生与孕妇赋予漏诊/假阳的相对成本，都存在不可忽略的方差。本文要回答的正是题目隐含的稳健性要求——当输入参数不确定时，最优方案本身会不会漂移？漂移有多大？如何给出一个"无论参数落在哪个区间都不至于太糟"的决策？本文用蒙特卡洛传播把这个问题从"一个答案"提升为"一组带置信的答案"。

二、模型假设

参数不确定是可建模的：四类关键参数各自服从一个围绕标称值的对数正态分布，相对波动幅度可设定（患病率 ±30%、区分度斜率 ±25%、代价比 ±40%、羊穿风险 ±30%）。
确定性内核复用范文一：给定任意一组参数，仍有唯一最优设计 $(g^*,z^*,E^*)$ ，其求解沿用范文一的时点—表现模型与期望代价最小化。
决策者可接受"区间"而非"单点"：临床上追求的是把最优时点与阈值锁定在一个可操作的稳健区间内，而非某个精确值。
不确定性相互独立：各参数的波动互不相关（年龄别患病率与测序区分度之间无强耦合），可用独立采样近似联合分布。
合成数据仅用于方法演示，真实落地须以各中心实测的参数分布为准。

三、符号说明

符号	含义
$g$	检测孕周（周）
$z$	风险报告 cut-off（标准化阈值）
$P(D)$	年龄别患病率先验
$\Delta(g)$	区分度（AUC 的单调函数，随孕周上升）
$r=C_{\mathrm{FN}}/C_{\mathrm{FP}}$	假阴对假阳的代价比
$r_{\mathrm{amn}}$	羊穿流产风险
$g^,z^,E^*$	最优孕周、最优 cut-off、最小期望代价
$\mu,\sigma,\mathrm{CI}_{95}$	MC 样本的均值、标准差、95% 区间

四、建模过程

4.1 不确定性源与采样

与范文一"给定参数求一点"不同，本文先对参数建模。以年龄 40 为例，患病率基准 $P(D)=0.08$ ，四类参数按对数正态扰动：

$X = \exp\big(\log x_0 - \tfrac{1}{2}s^2 + s\,\varepsilon\big),\quad \varepsilon\sim\mathcal{N}(0,1),$

其中相对波动系数 $v$ 决定 $s=\sqrt{\log(1+v^2)}$ 。这样参数恒为正、且围绕标称值 $x_0$ 按幅度 $v$ 波动。图 1 和图 2 是 4000 次 MC 得到的 $g^*$ 与 $z^*$ 分布——它们不再是两个孤立点，而是两簇分布：

图 1 显示最优孕周几乎全部落在 13.5–14.0 周，非常集中（标准差 0.25）；图 2 的最优 cut-off 则在 1.93–2.47 间展开（均值 2.20）。结论一：最优时点对参数不确定高度稳健（波动不到 0.5 周），而 cut-off 的波动更明显——时点可放心锁定，阈值则需留出保守余量。

4.2 确定性内核：在给定参数下求最优设计

对每次采样的参数 $(P(D),\Delta\text{-斜率},r,r_{\mathrm{amn}})$ ，我们沿用范文一的"时点—表现—代价"内核求唯一最优。期望负担 $E(g)$ 拆为三项：

$E(g) = \underbrace{p_{\mathrm{fail}}(g)\,C_{\mathrm{nocall}}}_{\text{失败重采代价}} + \underbrace{(1-p_{\mathrm{fail}}(g))\,C_{\mathrm{dec}}(g)}_{\text{成功检测的决策代价}} + \underbrace{k_{\mathrm{op}}(g-G_{\min})}_{\text{可操作性惩罚}},$

其中 $p_{\mathrm{fail}}(g)$ 是失败率（孕周越小越高）、 $C_{\mathrm{dec}}(g)$ 是给定区分度下对该孕周做最优 cut-off 的最小决策代价、 $k_{\mathrm{op}}$ 是越晚可操作窗口越窄的惩罚。对 $g$ 在 $[9,24]$ 内每 0.5 周扫描、对 $z$ 在 $[0.5,4.5]$ 内精细扫描，取全局最小值即该组参数下的 $(g^*,z^*,E^*)$ ——这就是叠加在 MC 外层的内层优化，复用范文一确定性解保证一致性。

4.3 最优时点对不确定度有多敏感

为了看清"到底哪个参数最能扰动决策"，我们在四类参数里逐维下调到标称的 60%–75%（其余参数保持标称），看最优时点 $g^*$ 相对 MC 均值 13.77 周移动多少：

图 5 的 tornado 揭示：各维度单独下调对最优时点的扰动都在 ±0.3 周以内——患病率先验 −30%、代价比 −40%、羊穿风险 −30% 都把它略推前约 0.3 周，而区分度斜率 −25% 反而将其略推后约 0.2 周。这说明最优时点对单维参数并不剧烈敏感（标准差仅 0.25 周），没有任何一个单个参数能"独大"地左右时点；真正驱动 cut-off 波动的，是图 3 展示的"患病率先验→cut-off 单调响应"（患病率越高阈值越宽，造成 MC 中 z* 的 1.9–2.5 范围）。这把"该优先较真哪个参数"的答案也变清晰：与其纠结局部权重的微小差别，不如把注意力放在校准年龄别患病率基线（它是 cut-off 波动的主源）与区分度的稳定性上。

4.4 决策曲面：时点 × 代价比

把"时点 $g$ "与"代价比 $r$ "两个最主要维度铺开成平面，逐格计算期望负担 $E$ ，得到决策最热力图：

图 7 每格是按该 $(g,r)$ 组合对全剖面求解的最优期望代价，颜色越绿越低。可以看到最优带（谷）沿"中早孕周 × 中高代价比"斜向延伸，代价比越大（越重视漏诊）时最优时点会向偏早一侧靠拢。这张曲面直接支撑 4.5 的稳健建议——不存在唯一"最优点"，存在一条稳健的"可接受带"。

4.5 确定性解 vs 鲁棒解

对照范文一的确定性单点解（ $g^*=14.0,z^*=2.16,E^*=4.4765$ ）与本文 MC 的鲁棒解，图 4 进行同框对比：

方案	$g^*$	$z^*$	$E^*$
确定性单点解	14.0 周	2.16	4.4765
鲁棒中点解	13.8 周	2.20	4.6675
MC 均值	13.77 周	2.20	4.5688

确定性解略偏"晚一点、激进一步"，鲁棒解则略偏"早一点、保守一步"。二者都在 $g^*$ 的 95% 区间 $[13.5,14.0]$ 内，说明范文一的点估计方向正确，只是没有告诉读者那里的不确定性有多宽——鲁棒分析把这块短板补上：真正可操作的结论是"最优孕周锁定在 13.5–14.0 周、cut-off 取 1.9–2.5 区间中偏保守端"，而不是某一孤点。

4.6 期望代价的分布与风险上界

图 6 给出 $E^*$ 的累积分布：中位约 4.57，95% 区间 $[3.767,5.568]$ 。这意味着即使参数偏移到边缘，每筛查一人的期望代价也很少突破约 5.6（相对中位仅上浮约 22%）。期望代价的 MC 离差（标准差 0.56，约 12%）远小于参数本身的波动（30%–40%）,这来自内层优化的"吸收"作用——最优设计会自动让高代价参数下选择更省的窗口，天然起到缓冲。这是蒙特卡洛相较单点分析的第三大价值：量化了最坏情形。

4.7 时点与阈值的联合不确定性

单独看 $g^*$ 与 $z^*$ 的分布还不够，实际决策要同时定下"何时测 + 用什么阈值报风险"。图 8 把 4000 组最优解 $(z^*,g^*)$ 在同一平面铺开：

图 8 呈明显的负向斜带： $z^*$ 与 $g^*$ 的相关系数 $r\approx-0.49$ ——最优孕周偏早时，cut-off 反而偏大（更保守），反之孕周偏晚时阈值可略放宽。背后的机理是：早孕周区分度较低、假阳相对多，若要压住期望代价就得牺牲一点灵敏度、把阈值抬高；晚孕周区分度充分，可用较宽容阈值。单一"早测 + 高灵敏"的组合是最脆弱的，联合分布把这层权衡直观暴露，也让落地时能同步给出"时点-阈值"配套建议，而非分别拍两个孤立的数。

五、灵敏度分析

患病率先验对 cut-off 的单调响应（图 3）：把患病率基准从事 0.02 推到 0.14，最优 cut-off $z^*$ 从 2.45 单调降到 2.02、最优时点则从 13.5 周缓升至 14.0 周：

图 3 表明患病率越高，越该用较宽的阈值（少把阳性当确诊）、并略偏晚检测。这也解释了图 2 中 $z^*$ 为何波动达 1.9–2.5 还偏保守——患病率先验的不确定性是 cut-off 波动的首要来源，落地时应优先校准年龄别患病率基线。

区分度斜率与代价比：代价比从 5 拉到 40，也把 $z^*$ 从 2.44 单调压到 2.00（同向于患病率效应），说明"越怕漏诊越要宽阈值"的权衡对两个主参数一致；而区分度单维下调只把 $g^*$ 推后约 0.2 周（图 5），对时点影响最温和——经费与质控优先投给测序区分度的稳定，仍是最稳的选择。
MC 收敛：把样本数从 1000 加到 8000， $g^*$ 的均值与 95% 区间变化小于 0.1 周，4000 次已充分收敛，结果不是随机噪声。

六、模型评价与改进

优点：① 把"参数不确定"明式建模，回答了题目隐含的稳健性问题——单点解不是答案，带置信的解才是；② tornado 直接指出最敏感参数，指导资源投向（区分度校准优先于代价权重）；③ 期望代价 MC 分布给出最坏情形上界，供临床做风险预案。局限与改进：① 假设参数相互独立，真实中患病率与年龄强相关、代价比与中心定价相关，可用 Copula 或参数化相关矩阵改进；② 未把范文二的多指标分类不确定性纳入联合传播，可扩展为"参数 + 模型"双层 MC；③ 对数正态扰动仅是近似，可用实测分布或 bootstrap 医学文献数据校正；④ 收敛基于均值的经验判据，严格做法是 Gelman-Rubin 诊断。综合来看，本文为筛查方案设计提供了一套可迁移的鲁棒决策范式：先敏感性定"该较真哪个参数"，再 MC 给"置信区间"，最后取"中等保守"解落地。

七、结论

NIPT 筛查的最优设计不能只看一组确定性参数。把患病率、区分度、代价比、羊穿率四类不确定性做 4000 次蒙特卡洛传播后：最优孕周稳健锁定在 13.5–14.0 周（均值 13.77、std 0.25），最优 cut-off 均值 2.20（95% 区间 1.93–2.47）、宜取偏保守端，期望代价 95% 上界约 5.57、对参数波动天然有缓冲。决定性发现是区分度斜率比代价权重更能扰动时点决策——经费与质控应优先投给测序区分度的稳定。方案落地建议采用鲁棒中点解并持续监测参数漂移，这比依赖名义参数求得的单点解在真实世界中更不容易翻车。

参考文献

全国大学生数学建模竞赛组委会. 2025 高教社杯 C 题赛题. http://www.mcm.edu.cn
Fishman G. S. Monte Carlo: Concepts, Algorithms, and Applications. Springer, 1996.
Saltelli A., et al. Global Sensitivity Analysis: The Primer. Wiley, 2008.
Palisade / @RISK 类方法论综述：不确定性下的决策与风险量化. Oper Res, 1996.
Bianchi D. W., et al. 无创产前检测（NIPT）临床应用综述. NEJM, 2014.

附录：核心 Python 实现

import sys, os
_HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, os.path.abspath(os.path.join(_HERE, "..", "..", "..", "tools")))
import gen_cumcm2025c_3 as M   # 真源：MC 不确定性传播 + 8 张 SVG

# 1) 蒙特卡洛仿真（复现真源 main()：N=4000，固定种子 2026）
gs, zs, es = M.mc_robust(4000)
stat = M.summarize(gs, zs, es)
print("MC 最优孕周 g*: 均值=%.2f 周  标准差=%.2f  95%%区间=[%.1f, %.1f]" %
      (stat["g_mean"], stat["g_sd"], stat["g_lo"], stat["g_hi"]))
print("MC 最优 cutoff z*: 均值=%.2f  标准差=%.2f  95%%区间=[%.2f, %.2f]" %
      (stat["z_mean"], stat["z_sd"], stat["z_lo"], stat["z_hi"]))
print("MC 最优期望代价 E*: 均值=%.4f  标准差=%.4f  95%%区间=[%.4f, %.4f]" %
      (stat["e_mean"], stat["e_sd"], stat["e_lo"], stat["e_hi"]))

# 2) 确定性单点解对照（范文一标称参数）
det_g, det_z, det_e = M.deterministic_solution()
print("确定性单点解(g*=%.1f, z*=%.2f, E=%.4f)" % (det_g, det_z, det_e))

# 3) 敏感性 tornado：单维下调对最优时点的影响
base_g = stat["g_mean"]
print("敏感性（Δg*/周）: 患病率-30%%=%.1f 区分度-25%%=%.1f 代价比-40%%=%.1f 羊穿-30%%=%.1f" % (
    M.optimal_design(M.BASE_PD * 0.7, M.AMP_MU, M.C_FP, M.R_MU * M.C_FP)[0] - base_g,
    M.optimal_design(M.BASE_PD, M.AMP_MU * 0.75, M.C_FP, M.R_MU * M.C_FP)[0] - base_g,
    M.optimal_design(M.BASE_PD, M.AMP_MU, M.C_FP, M.R_MU * 0.6 * M.C_FP)[0] - base_g,
    M.optimal_design(M.BASE_PD, M.AMP_MU, M.C_FP * 0.7, M.R_MU * M.C_FP)[0] - base_g))