公民科学物候数据中估计量选择造成的趋势偏倚：以 2015–2025 记录量爆发期为自然实验，比较首现日、分位数与 Weibull 校正估计量给出的昆虫出现期位移

1 · 研究问题

在 GBIF 上 2015–2025 年公民科学记录量增长约一个数量级的自然实验条件下，用"首次观测日（first observation date）""第 10/50 百分位日""Weibull 分布校正的起始日"三类物候估计量对同一批昆虫物种计算年际趋势，所得的趋势方向与幅度差异有多大？这些差异能否被各物种各年的记录数增长速率定量解释？

2 · 研究背景与空白

背景。 物候（phenology）位移是气候变化最直接的生物学指标之一。公民科学平台（iNaturalist、eBird 等）经 GBIF 汇聚了海量带日期与坐标的观测记录，成为物候研究的主力数据源。但这类数据没有标准化的采样设计：某物种某年的"首次被观测到的日期"同时取决于该物种何时真正出现，和当年有多少人在看——观测努力（observer effort）越大，首现日越早。这是一个经典的极值统计量偏倚问题。

已有工作到哪一步。 该偏倚在文献中被反复指出：加拿大 PlantWatch 研究（Scientific Reports, 2013）为控制观测者差异，只保留"至少 5 名观测者、至少 5 个地点"的首花日记录，报告 19 种植物 2001–2012 年首花日平均提前约 9 天；丹麦的对比研究显示定向公民科学项目（目标就是报告"首花"）记录到的花期显著早于标本馆与随机记录；Belitz 等（2020）提出并实现了基于 Weibull 分布的物候分位数估计量（R 包 phenesse），可给出对采样强度不敏感的起始日估计与置信区间。

空白在于：这些工作要么只做努力筛选、要么只提出更好的估计量，没有人把"估计量选择"当成自变量，在一个记录量剧烈变化的时间窗上量化它对趋势结论本身的影响。 换句话说，文献知道首现日有偏，但没有回答"这个偏差有多大、会不会把一个真实的物候提前趋势放大成两倍，或者把没有趋势的物种伪造出趋势"。2015–2025 恰好是 iNaturalist 记录量爆发的十年，构成一个天然放大器。这个空白适合学生课题：数据完全开放、算力可控、结论正负都成立。

3 · 可检验假设

• H1：以"首现日"估计的年际趋势斜率（天/年）系统性地比 Weibull 校正起始日估计的斜率更负（更"提前"），跨物种的中位差值 ≥ 1.0 天/年，且配对自助法 95% 置信区间不跨 0。
• H2（机制假设）：物种层面的斜率差值与该物种 2015–2025 记录数的年均增长率呈正相关（Spearman ρ ≥ 0.5，BH-FDR < 0.05）；若把每物种每年的记录数随机稀释到该物种十年内的最小年记录数，斜率差值应缩小 ≥ 60%。若稀释后差值不缩小，说明偏倚来自记录的空间/物种构成变化而非单纯数量，需另作解释。

4 · 量化验收标准

1. 方法学校验（硬门槛）：用自建管线复现两组已发表结果：(a) 在 Belitz 等（2020）phenesse 论文所用的模拟设定下重跑其 Weibull 估计量，要求覆盖率（nominal 95% CI 的实际覆盖率）与原文报告值偏差 ≤ 3 个百分点；(b) 选取 1 篇已发表的公民科学物候趋势论文（须公开物种清单、区域与年份范围），在同源数据上重算其报告的趋势斜率，要求 ≥ 70% 的物种斜率符号一致且中位斜率的绝对偏差 ≤ 0.5 天/年。任一组不达标，管线不可信，后续全部结论无效。
2. 统计口径预先写死：(a) 趋势用物种级线性回归（估计日 ~ 年）与跨物种混合效应模型（lme4，物种为随机截距与随机斜率）双口径都报；(b) 跨物种的多重检验一律报 BH-FDR，不报裸 p 值；(c) 所有估计量的不确定性用 ≥ 500 次自助法给 95% 置信区间，phenesse 的内置区间与自助区间都报；(d) 时间序列性质必须尊重——趋势检验不得做随机划分交叉验证，若做预测评估须按年份前后划分并额外报"随机划分 vs 时间划分"的差值以量化高估幅度；(e) 空间混杂必须控制——每物种限定在一个固定的地理网格（如 100 km × 100 km 或固定行政区）内，并报告各年记录的空间质心漂移，漂移超过 50 km 的物种-年标记并做敏感性分析。
3. 样本量与检验力预先估计：≥ 40 个物种，每物种每年 ≥ 30 条记录且 ≥ 8 个有效年份。开工前用模拟估计：在 n=40 物种、10 年序列下，能以 80% 检验力检出的最小斜率差值是多少；若大于 H1 的 1.0 天/年阈值，必须先扩大物种数或放宽年份要求。
4. 数据来源地位：GBIF 记录中的物种鉴定（含 iNaturalist 的 research-grade 标记）仅作为标签来源使用，不计入本项目的数据贡献；必须报告各物种的鉴定可靠性代理指标（research-grade 比例）并做敏感性分析。
5. 可复现性：GBIF 下载必须通过 rgbif::occ_download 生成带 DOI 的可引用下载，DOI 写入论文；全部脚本、物种清单、随机种子、区域定义开源；提供 10 物种的降规模复现包（≤ 1 小时可跑完）。

5 · 数据与工具

用途	来源 / 工具
观测记录	GBIF（https://www.gbif.org/ ），经 R 包 rgbif 的 occ_download 提交筛选式：basisOfRecord = HUMAN_OBSERVATION、hasCoordinate = TRUE、hasGeospatialIssue = FALSE、year = 2015,2025、限定目标分类群与地理范围。限定单一区域 + 单一目（如鳞翅目或蜻蜓目）后，下载体积通常在数百 MB 至 3 GB（DwC-A 压缩包），data.table::fread 分块读取即可，16 GB 内存充裕。全球全类群下载（数十 GB 至 TB 级）属超出范围
为什么选昆虫成虫出现期而非植物花期	昆虫成虫的观测日期本身即物候信号，无需依赖注释字段；而植物花期需要 iNaturalist 的 "Plant Phenology" 注释，该注释是否完整传递到 GBIF 的 DwC-A（reproductiveCondition / dynamicProperties 字段）需核实。若核实为可用，则可把植物花期作为第二个类群加入（见 go/no-go）
物候估计量	R 包 phenesse（Belitz et al. 2020）：weib_percentile() 与 weib_percentile_ci()。单次估计（n≈100 观测、iterations=500）约数秒；40 物种 × 11 年 × 3 分位点 × 500 次自助 ≈ 数小时至一夜，纯 CPU 可行。另实现朴素首现日与经验分位数作为对照
趋势建模	R lme4（混合效应）、boot、p.adjust(method="BH")
气候协变量（可选）	ERA5-Land 或 WorldClim 月均温（公开下载），仅用于说明性分析，不作为主结论
观测努力代理	同区域同期全部该分类群记录数（不限物种）作为努力代理；亦可用 eBird/iNaturalist 的观测者数。eBird Basic Dataset 需提交免费申请并等待审批，虽通常获批但存在时间不确定性，因此不设为主线依赖，仅列为可选补充
对照基准	Belitz 等 2020 的模拟覆盖率结果、所选复现论文的趋势斜率。仅用于校验与对比，不计入本项目的数据贡献
伦理	全部为已公开的观测记录，不涉及人体数据、不涉及活体动物操作、不涉及可识别个人信息（GBIF 记录的观测者字段在分析中一律不使用），无需伦理审批

6 · 方法路径

1. 装环境（R + rgbif + phenesse + lme4），跑通 phenesse 内置示例，先完成验收标准第 1 条的两组复现，把覆盖率表与斜率对比表写入 validation/。
2. 按预先写死的规则确定研究区域、分类群与物种清单（规则：每年 ≥ 30 条记录、≥ 8 个有效年份、research-grade 比例 ≥ 0.8、空间质心年漂移 < 50 km），提交带 DOI 的 GBIF 下载。规则先于数据探索确定，不得看到结果后再改。
3. 核实注释字段可用性（决定是否加入植物花期第二类群）与各物种的记录量增长曲线，出具"物种 × 年 × 记录数"矩阵与空间质心漂移表。
4. 对每个物种每年计算三类估计量（首现日、第 10/50 经验百分位、Weibull 校正第 10/50 百分位）及其自助区间。
5. 分别用三类估计量拟合年际趋势（物种级 OLS + 跨物种混合效应），得到三套斜率，做配对比较并给自助区间，检验 H1。
6. 检验 H2：把斜率差值对记录数年均增长率做回归；再做稀释实验（每物种每年下采样到最小年记录数），重算三套斜率与差值。
7. 独立交叉校验：换一套采样设计完全不同的数据源复算至少 10 个物种——首选具备标准化协议的公开数据（如英国 UKBMS 蝴蝶监测的公开年度指数、NEON 的地面节肢动物取样数据），检验"标准化数据上三类估计量差异是否消失"。若消失，则确证差异源于公民科学的采样不均而非估计量本身的数学性质。

7 · 新颖性边界

本课题不提出新的物候估计量（Weibull 校正估计量由 Belitz 等 2020 提出并实现），不声称发现新的物候位移现象，不声称"首现日对采样强度敏感"为本项目发现——这一点自 Miller-Rushing 等以来是物候统计学的常识，PlantWatch 与丹麦对比研究都已明确报告。

已有工作具体完成了什么：一类工作用观测者/地点数下限筛掉有偏记录（PlantWatch，2001–2012，19 种植物，报出约 9 天提前）；另一类工作提出更好的估计量并用模拟验证其覆盖率（phenesse，2020）；还有工作对比了不同数据源（标本馆 vs 定向公民科学 vs 随机记录）的花期分布差异。没有一项工作把"估计量选择"当作自变量，在记录量剧烈增长的真实时间窗上量化它对最终趋势结论的偏倚幅度。

本项目的贡献（且是主结论）：利用 2015–2025 公民科学记录量增长约一个数量级这一自然实验，给出"估计量选择造成的趋势偏倚"的定量刻度（天/年），并检验该偏倚是否可由记录量增长率解释（H2）与是否在标准化数据上消失（第 7 步对照）。这是评价维度的转换——已有工作评估的是估计量的统计性质（覆盖率、偏差），本项目评估的是该选择对最终气候生物学结论的偏置。

为什么有价值：大量近年物候论文直接使用 GBIF/iNaturalist 的首现日。若偏倚达到 1 天/年量级，在十年窗口上就是 10 天——与文献报告的真实位移同量级，足以改变结论的解释。

风险声明："三类估计量给出的趋势无显著差异"同样是有效结论，它会为大量已有物候文献提供一次重要的稳健性背书。但必须由第 3 条的检验力估计证明本设计有能力检出 1.0 天/年的差值，否则该结论无效。

8 · 决策门槛（go / no-go）

• 第 5 周末：确认满足筛选规则的物种数。若合格物种 < 25，立即降级：其一，把区域扩大（从单省/单州扩到国家或大区）并把年记录数下限降到 20；其二，把时间窗前移到 2010–2025 以纳入更多年份。两条降级都保留"三类估计量并列 + 偏倚量化 + 增长率解释"的主结论框架。
• 第 8 周末（硬门槛）：验收标准第 1 条必须通过。若 phenesse 覆盖率复现偏差 > 3 个百分点，先核对自助迭代次数与随机种子；若已发表趋势论文的斜率复现符号一致率 < 70%，须逐项排查区域定义与物种同义名（GBIF 的分类骨架会合并同义名，这是最常见的不一致来源）。六周内仍不通过则降级：把校验目标改为"在模拟数据上验证三类估计量的已知理论性质"（首现日随 n 增大而单调提前、Weibull 估计量渐近无偏），阈值改为模拟结果与理论预期一致。主结论框架不变。
• 第 12 周末：核实植物花期注释字段是否可用。若可用则加入第二类群（增强结论普适性）；若不可用则明确放弃植物分支，只做昆虫，并在论文中写明原因——不得含糊带过。
• 第 24 周末：完成第 7 步的标准化数据对照。若找不到与目标物种重叠的标准化监测数据，立即降级：改用"高观测努力年份的子集"作为伪标准化对照（在记录最密集的 3 年内做稀释，模拟低努力情形），并标注这是内部对照而非独立数据源。
• 第 36 周结果冻结，第 44 周英文 PPT 初稿。
• 需提前核实而非边做边发现：(a) GBIF 下载配额与 DwC-A 体积（先用 occ_count 预估）；(b) iNaturalist 物候注释在 GBIF 中的传递情况；(c) phenesse 在 n < 20 时是否会失败或给出发散区间（须实测并写入物种筛选规则）。三项在第 8 周前查清。
• 选择前提：适合对统计推断有兴趣、能接受"结论是一把尺子而不是一个发现"的学生。答辩时须能用英文清楚说明"为什么估计量的选择会改变生态学结论"，这是本课题的说服力核心。
• 预算裁剪顺序：先砍气候协变量分析，再砍稀释实验（H2 后半），再砍第 50 百分位（只保留第 10 百分位与首现日的对比），最后把物种数从 40 降到 25。砍到只剩"25 物种 × 首现日 vs Weibull 第 10 百分位 × 混合效应趋势 + 增长率回归"时，主结论仍成立。

---