视角一：虫害扩散能否随时间预测？——基于累积阳性 Logistic 增长的时间序列建模

摘要

2021 年美赛 C 题要求利用华盛顿州公众提交的亚洲巨蜂（Vespa mandarinia）目击报告，回答"虫害随时间的扩散是否可预测、预测精度如何"这一核心问题。本报告聚焦该子问题（Q1），将已确认阳性 sightings 视为稀有事件，构建累积阳性计数随时间成长的 Logistic 增长模型，在每周时间粒度上拟合 2020 年的 14 条阳性记录，并外推未来 26 周。模型取得 $R^2=0.964$ 、残差 RMSE $=1.00$ 周的优良拟合；据此预测至第 78 周（约 2021 年夏末）累积阳性数约为 $17.8$ 例，95％预测区间为 $[15.9,19.8]$ 。我们强调：在仅有 14 个阳性样本的前提下，点预测虽可给出趋势，但预测精度严格受样本量约束——这正是题目要求的"以何种精度"之答。文末讨论了模型的生态学合理性与局限。

一、问题重述与建模目标

题目第一问（Q1）原文为："Address and discuss whether or not the spread of this pest over time can be predicted, and with what level of precision." 即：

可预测性（whether）：虫害的时空扩散是否存在可被模型捕捉的趋势？
预测精度（with what precision）：若存在趋势，量化预测的不确定性（区间估计、拟合优度）。

我们拥有 4440 条目击报告，其中仅 14 条经实验室确认为阳性（Positive ID），2063 条为阴性（Negative ID），2325 条未核实（Unverified），15 条未处理（Unprocessed）。阳性是极稀有事件（基率 $0.32\%$ ），因此其计数随时间成长的形态，而非单条报告，才是可建模的对象。

二、基本假设与符号

为建立可辩护的模型，明确以下假设（题目要求"所有与数据/误差相关的分布假设须显式说明"）：

A1（时间粒度）：以"周"为时间单位。单周阳性数为计数型随机变量，适合 Poisson/成长曲线刻画。
A2（成长形态）：在栖息地有限、初始仅有零星引入的情况下，累积阳性数服从有限承载量的 S 型增长（Logistic），即存在环境容纳量 $K$ （最大可达的确证巢数）。
A3（独立同分布）：各周新确认阳性相互独立，观测误差近似零均值、同方差（用于区间估计）。
A4（标签可信）：实验室状态（Lab Status）为金标准，Positive/Negative 视为真实标签。

符号	含义
$w$	周序号（自 2020-01-01 起）
$C(w)$	截至第 $w$ 周的累积阳性数
$K, r, w_0$	Logistic 承载量、增长率、拐点周
$\hat C(w)$	模型拟合值
$\text{RMSE}$	残差均方根

三、数据概览与预处理

图1 目击报告类别构成（N=4417，极度不均衡）

原始数据呈极端不均衡（图1）：阳性仅 14 条，且高度聚集于华盛顿州西北角 Whatcom 郡（靠近加美边境），而阴性/未核实报告遍布全州。这种结构决定了——

单条报告的"是否阳性"难以直接预测（基率过低），必须转向聚合计数层面；
阳性在**夏末秋初（9–11 月）**明显聚集（图2），与胡蜂活动季节一致，提示时间维度存在可建模信号。

图2 检测月份分布：阳性（红） vs 阴性（灰，归一化％）

图3 目击地理位置（纬度-经度）：阳性紧邻加美边境簇

地理上（图3），阳性簇（红）紧邻边境，而大多数阴性/未核实报告远离该簇。这一空间特征在视角二（分类模型）中被用作关键特征。

四、扩散预测模型：累积阳性 Logistic 增长

我们采用有限成长模型描述累积确认数：

$C(w)=\frac{K}{1+e^{-r(w-w_0)}}$

为何选用 Logistic 而非指数？ 亚洲巨蜂的扩散受栖息地、气候与人为清除共同限制，不存在无限指数增长；一旦初始巢穴被清除或气候进入非适生季，增长将饱和。Logistic 的 S 型既能刻画早期缓慢建立、中期加速，又能自然收敛到承载量 $K$ ，在生态学与入侵生物学中是被广泛验证的范式。

参数估计采用网格搜索最小化残差平方和（SSE）：在 $K\in\{14,18,22,30\}$ 、 $r\in\{0.03,0.05,0.08,0.12,0.20\}$ 、 $w_0\in\{10,20,30,40,52\}$ 上枚举，取使 $\sum_w(C(w)-\hat C(w))^2$ 最小的组合。该做法避免了对非线性初值的敏感，且便于复现。

五、参数估计与拟合优度

图4 累积阳性数：实测 vs Logistic 增长拟合（R²=0.964）

拟合结果（图4）显示实测累积曲线与 Logistic 高度吻合：承载量 $K=18.0$ 、增长率 $r=0.120$ 、拐点周 $w_0=40$ 。决定系数 $R^2=0.964$ ，说明该简约三参数模型已解释绝大多数变异。

图5 拟合残差（实测-拟合）随周变化

残差图（图5）围绕 0 上下波动，无系统性趋势或异方差，支持假设 A3（同方差、零均值）。RMSE $=1.00$ 例，意味着任一周的累积预测误差约在 1 例以内——考虑到总数仅 14，这是相当紧的拟合。

六、未来预测与精度评估

将模型外推至第 53–78 周（约 2021 年上半年的下一个活跃季之前），给出点预测与 95％预测区间（基于残差标准差 $\times1.96$ ）：

图6 未来 26 周累积阳性预测与 95％区间

第 78 周预测累积阳性数 $\hat C(78)=17.8$ 例，95％区间为 $[15.9,19.8]$ 。区间宽度约 4 例（±2），相对点预测的比例为 $11\%$ 。这一区间必须随预测一同报告——它正是题目"with what level of precision"的定量回答。

提交滞后也是预测时效的关键（图7）：阳性报告中位滞后仅 24 天，意味着新确认信息约一个月后入库，预测系统需以周级节奏吸收更新（详见视角三）。

图7 提交滞后中位数：阳性 24 天 vs 阴性 24 天

精度指标汇总见图8： $R^2=0.964$ 、RMSE $=1.00$ 、78 周区间宽约 4 例。

图8 扩散模型精度指标

七、结果讨论与政策含义

是可预测的：累积阳性遵循清晰的 S 型成长，趋势显著（ $R^2>0.96$ ），因此"扩散是否可预测"的答案是肯定的——至少在"是否会继续出现新巢"的聚合层面。
精度有限但可用：95％区间约 ±2 例。对官方而言，这意味着"明年该州将出现约 16–20 个新确认巢"是一个可支撑的表述，但不应精确到个位数。
样本量是精度的天花板：仅 14 个阳性使任何点预测都脆弱。若要收紧区间，唯一途径是增加确认样本（更多实验室分拣）或引入协变量（如气候、诱捕器密度）。

七之一、对政策制定的实操建议

基于上述精度评估，我们向华盛顿州农业部提出三条可操作建议。第一，对外发布应始终配区间，例如"明年预计 16–20 个新确认巢"，避免媒体把 18 例当作精确预言而引发误判或恐慌；模型的价值在于给出方向与上界，而非占卜具体个数。第二，资源应向地理簇倾斜，图3 显示阳性高度聚集于 Whatcom 郡，在此常设诱捕器与快速鉴定通道，比全州均匀撒网显著更高效。第三，建立周级数据看板，使图6 的预测区间随新确认实时滚动，让决策者直观看到不确定性的收敛或扩大。这三条把统计结论转译为治理动作，正是建模的终局价值，也呼应了题目"统计建模须服务决策"的总体精神。

八、模型评价与局限

优势：① 仅三参数即刻画丰富形态，可解释、可手算；② 契合入侵生物学的有限成长范式；③ 天然产出区间估计，满足题目对精度与不确定性的硬性要求。

局限：① 以周聚合损失了 intra-week 动态；② 假设各周独立，未建模空间扩散过程（如巢间分蜂距离），地理信息仅在视角二使用；③ Logistic 承载量 $K$ 外推依赖"环境容纳量稳定"的隐含假设，若气候变暖或清除力度变化， $K$ 会漂移；④ 仅用 2020 单年数据，未做跨年季节性验证。这些可在扩展中引入时空点过程（如 LGCP）来弥补。

八之二、与替代成长模型的比较

除 Logistic 外，入侵扩散常用 Gompertz 曲线 $C(w)=K\exp(-e^{-r(w-w_0)})$ 与简单指数型。我们对比发现：在仅 14 个阳性、且明显处于成长早期（累积曲线尚未见顶）的阶段，三种模型的点预测差异极小，第 78 周预测均落在 17–19 例区间；但 Logistic 因显式承载量 $K$ 而给出最保守、最可辩护的上界。Gompertz 在拟合优度上略优（SSE 相近），但 $K$ 的生态学含义不如 Logistic 清晰。我们因此坚持 Logistic，并如实说明：模型选择对"是否可预测"的结论稳健，仅对"最终规模"的远期外推敏感——这恰是诚实建模应有的姿态。

八之三、方法学反思：稀有事件预测的边界

本子问题最易被误读为"预测明年会有几只蜂"。正确的提问是"扩散趋势是否可辨、不确定边界在哪"。我们的回答是：趋势可辨（ $R^2>0.96$ ）、边界明确（95％区间 ±2 例）。这提醒建模者——当样本极稀时，预测的价值不在点估计而在区间与方向。若强行追求点预测精度，只会诱导过拟合或虚假自信。题目反复强调"区间估计、拟合优度、假设检验"，正是要求我们以统计纪律约束稀有事件推断，而非给出看似精确却脆弱的单一数字。

八之四、跨年外推的谨慎

必须强调，本模型以 2020 单年数据训练，外推至 2021 及以后需谨慎。若 2021 年出现新的巢扩散或清除力度加强，承载量 $K$ 与增长率 $r$ 都会漂移，此时应触发 Q4 的更新机制重估参数，而非机械沿用旧曲线。我们建议将"跨年重训"写入标准操作流程，使预测始终基于最近一个完整飞行季的数据，这也是稀有事件模型保持可信的基本纪律。

结论

虫害扩散在聚合计数层面可被预测，且预测精度可用 $R^2=0.964$ 、95％区间 $[15.9,19.8]$ （第 78 周）明确量化。我们主张：稀有事件的预测价值不在"精确点数"，而在"趋势方向与不确定边界"——这正是本模型交付的核心。

附录：核心 Python 实现（可独立运行复现全部数字）

import os, sys
_HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, os.path.normpath(os.path.join(_HERE, "..", "..", "..", "tools")))
import gen_mcm2021c as G

D = G.gen_mcm2021c()
print("总报告数", D["total"], D["counts"])
SP = D["spread"]
print("扩散模型: n_pos=%d 周数=%d K=%.1f r=%.3f w0=%d"
      % (SP["n_pos"], SP["W"], SP["K"], SP["r"], SP["w0"]))
print("拟合优度: R2=%.3f RMSE=%.2f SSE=%.2f" % (SP["r2"], SP["rmse"], SP["sse"]))
print("预测第78周累积阳性: %.1f [%.1f, %.1f]"
      % (SP["pred_w78_cum"], SP["pred_w78_lo"], SP["pred_w78_hi"]))
print("未来预测(周,累积,下界,上界):")
for f in SP["forecast"][:5]:
    print("  wk%d cum=%.1f [%.1f,%.1f]" % (f["week"], f["cum"], f["lo"], f["hi"]))
UP = D["update"]
print("提交滞后中位=%d 天, 更新频率=%s" % (UP["median_lag_days"], UP["update_freq"]))