名师讲堂｜使用 Python 测算各区县鸟类丰度指数

今天给大家分享使用 Python 测算各区县鸟类相对丰度指数的方法。配套的数据来自中国观鸟记录中心（RStata 数据中心爬取），核心方法是复现 Liang et al. (2020, PNAS) 提出的「努力量校正后的相对鸟类丰度」指数。本次以 2015–2020 年 的观鸟数据为例进行演示——样本窗口是出于示例目的选定的，并非参照某篇论文来确定。

在 R 中我们直接用 fixest 包；本讲义改用 Python 实现：用 pandas 做数据整理、用 pyfixest 的 fepois 跑带大量固定效应的泊松回归（对应 R 的 fixest::fepois），用 matplotlib 出图，再用 reticulate 在 R Markdown 里管理一个专属 Python 虚拟环境来运行这些代码。

完整的观鸟数据可以从下面的链接获取：

1980~2025 年观鸟记录、经纬度及其所处的省市区县数据(2026年3月爬取)：https://rstata.duanshu.com/#/brief/course/0005692119984a0b90677782a8f0a333

下面先讲清楚这个指标到底是怎么算出来的，再逐行讲解计算代码。

一、指标来源与计算过程

1.1 论文背景：从空气污染管制到鸟类保护

Liang et al. (2020) 发表在 PNAS 的论文 *”Conservation cobenefits from air pollution regulation: Evidence from birds”* 想回答一个问题：美国《清洁空气法》削减空气污染，是否顺带保护了鸟类？

论文用的是 eBird 公民科学观鸟数据，覆盖全美各县。它的核心因变量不是「鸟类丰富度（species richness）」，而是 「努力量校正后的相对鸟类丰度」。这个指数把「观察者有多努力」「当时好不好看见」这些干扰剥离掉，只留下真正反映鸟类数量的信号，再用来和臭氧、PM2.5 做回归。

一句话：论文要的是「相对」丰度，不是「绝对」有多少只鸟——它只能告诉你 A 地比 B 地、今年比去年多还是少。

1.2 两步法（Methods: Bird Abundance Estimation）

1.3 本讲义如何用中国观鸟数据复现

数据来自 RStata 数据中心「1980～2025 年观鸟记录、经纬度及其所处的省市区县数据」，这里选择 2015～2020 年的数据样本演示。包含两份核心 dta：

• 观鸟记录表（report 级）：约等于论文的 “checklist”；
• 鸟种观测统计报告（report×物种级）：约等于 eBird Reference Dataset；

关键：两份数据里 鸟种数量 含义不同！report 级是「物种数」，report×物种级是「个体数」。论文 Eq.1 左边的 #birds observed = 明细表按 reportId 求和的个体数。

论文分类群划分（Rosenberg et al. 2019 口径，映射到中文目名）如下：

• 水禽 waterfowl：雁形目；
• 涉禽 shorebirds：鸻形目；
• 水鸟 waterbirds：鹈形目、鹤形目、䴙䴘目、鲣鸟目、鹳形目、潜鸟目、鹱形目、鹲形目、红鹳目；
• 其余目 = 陆鸟 landbirds。

本讲义以 2015–2020 年的数据为例进行演示。更早年份（2015 之前）数据极稀疏，县×月×年固定效应难以识别，故不纳入本示例窗口。

二、使用 reticulate 创建与管理 Python 虚拟环境

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间，避免与系统 Python（如 Anaconda）发生版本冲突。本文档的 setup chunk 已经抢先锁定了 .venv 的 Python；下面再展示安装与验证过程。

重要说明（避免”已初始化”报错）：reticulate 在 R 会话（或 knit 过程）中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错。因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径，这是 reticulate 选取 Python 的最高优先级入口。

安装 reticulate（仅首次）

# 设置 CRAN 镜像（knit 时 R 处于非交互模式，不会自动选择镜像）
options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

# 仅在尚未安装时才安装，避免每次 knit 都重装
if (!requireNamespace("reticulate", quietly = TRUE)) {
  install.packages("reticulate")
  message("reticulate 安装完成！")
} else {
  message("reticulate 已安装，版本：", packageVersion("reticulate"))
}

在虚拟环境中安装 Python 包（仅首次）

# 检查关键包是否已安装，缺失的才安装
py_pkgs <- c("numpy", "pandas", "pyarrow", "matplotlib",
             "pyreadstat", "pyfixest", "scipy")
installed <- py_list_packages(.venv_path)$package
need <- setdiff(py_pkgs, installed)

if (length(need) > 0) {
  virtualenv_install(.venv_path, packages = need)
  message("已安装缺失的包：", paste(need, collapse = ", "))
} else {
  message("所有 Python 包已就绪，无需安装")
}

说明：pyfixest 是 Python 里对应 R fixest 的包，提供 fepois（带高维固定效应的泊松回归）；pyreadstat 用于写回带中文变量标签的 Stata dta；matplotlib 负责出图。

验证激活状态

# 验证当前绑定的 Python 路径（应指向 .venv 目录）
py_config()

虚拟环境管理常用命令

# 查看所有已创建的虚拟环境
virtualenv_list()

# 删除虚拟环境（当不再需要时）
# virtualenv_remove(".venv")

三、详细讲解计算代码

整套计算分三步，对应三个 Python 脚本：01_筛选样本 → 02_测算丰度指数 → 03_对照刘钊论文。下面逐段给出完整代码并讲解。所有 Python 代码都在 .venv 中运行。

3.1 步骤 1：筛选 2015–2020 样本

观鸟记录表按 观测时间_起始 的年份直接筛选；鸟种观测统计报告本身无时间字段，需按 2015–2020 的 reportId 集合回连筛选。

1980~2025年观鸟记录、经纬度及其所处的省市区县数据(2026年3月爬取) 文件夹在附件中未提供，如需要可以从前述链接获取。

# ==============================================================================
# 步骤 1：从原始数据中筛选 2015–2020 年样本，导出到本项目文件夹
# 数据处理：微信公众号 RStata
# 数据来源：RStata 数据中心「1980–2025 年观鸟记录、经纬度及其所处的省市区县数据」
#   - 观鸟记录表(report 级)：直接按 观测时间_起始 的年份筛选
#   - 鸟种观测统计报告(report×物种级)：本身无时间字段，需按 2015–2020 的
#     reportId 集合回连筛选
# 产出：本项目文件夹下的两个样本 dta
#   - 观鸟记录_2015_2020.dta
#   - 鸟种观测统计报告_2015_2020.dta
# ==============================================================================

# 强制 stdout/stderr 走 ASCII：reticulate/knitr 捕获 Python 输出时硬编码 latin-1，
# 任何中文打印都会抛 UnicodeEncodeError 中断 knit。这里把输出中的非 ASCII 字符替换为
# '?'，作为兜底；正文讲解与 .py 源文件注释仍保留中文。
import sys as _sys
class _AsciiOut:
    def __init__(self, raw):
        self._raw = raw
    def write(self, s):
        if isinstance(s, str):
            s = s.encode("ascii", "replace").decode("ascii")
        return self._raw.write(s)
    def flush(self):
        return self._raw.flush()
    def __getattr__(self, n):
        return getattr(self._raw, n)
if not isinstance(_sys.stdout, _AsciiOut):
    _sys.stdout = _AsciiOut(_sys.stdout)
if not isinstance(_sys.stderr, _AsciiOut):
    _sys.stderr = _AsciiOut(_sys.stderr)

import os
import pandas as pd
import pyreadstat  # 写回 dta 时用 UTF-8，保留中文列名（pandas 原生 to_stata 默认 latin-1 会破坏中文列名）

dir_in = "/Users/ac/Desktop/计算地区生物多样性/1980~2025年观鸟记录、经纬度及其所处的省市区县数据(2026年3月爬取)"
dir_out = "/Users/ac/Desktop/使用 Python 测算各区县鸟类丰度指数"
os.makedirs(dir_out, exist_ok=True)

f_rep = os.path.join(dir_in, "1980～2025年观鸟记录、经纬度及其所处的省市区县数据(2026年3月爬取).dta")
f_sp = os.path.join(dir_in, "1980～2025年鸟种观测统计报告(2026年3月爬取).dta")

YEAR_MIN = 2015
YEAR_MAX = 2020

# ------------------------------------------------------------------------------
# 1) 观鸟记录表：按年份筛选
# ------------------------------------------------------------------------------
print("== Step 1: read checklist data and filter by year ==")
rep_all = pd.read_stata(f_rep)
# 观测时间_起始 在 read_stata 后通常为 datetime64，用 .dt.year 取年份
rep_all["year"] = pd.to_datetime(rep_all["观测时间_起始"]).dt.year
rep_all = rep_all[rep_all["year"].between(YEAR_MIN, YEAR_MAX)]

print("reports after filtering:", len(rep_all))
print(rep_all.groupby("year").size().sort_index())

out_rep = os.path.join(dir_out, "观鸟记录_2015_2020.dta")
pyreadstat.write_dta(rep_all, out_rep)
print("written -> sample checklist dta (2015-2020)")

# ------------------------------------------------------------------------------
# 2) 鸟种观测统计报告：按 reportId 关联筛选 （2015–2020 的报告才保留）
# ------------------------------------------------------------------------------
print("\n== Step 2: filter species-report by joining on reportId ==")
rep_ids = rep_all[["reportId"]].drop_duplicates()
print("unique reportId count (2015-2020):", len(rep_ids))

sp_all = pd.read_stata(f_sp).merge(rep_ids, on="reportId", how="inner")
print("species-report rows after filtering:", len(sp_all))

out_sp = os.path.join(dir_out, "鸟种观测统计报告_2015_2020.dta")
pyreadstat.write_dta(sp_all, out_sp)
print("written -> sample species-report dta (2015-2020)")

print("\nSample filtering done. Two sample dta files saved to the project folder.")

筛选结果如下：

print("reports after filtering:", len(rep_all))
print("species-report rows after filtering:", len(sp_all))
print(rep_all.groupby("year").size().sort_index())

3.2 步骤 2（上）：报告级计数与努力量变量

读取样本后，先把鸟种明细（report×物种级）按 reportId 汇总成「报告级」的总个体数 n_ind_all 与物种数 n_sp_all，并切出四个分类群；再把报告表整理出努力量（hours、n_observers）与地理时间维度，最后按 reportId 合并。

# ==============================================================================
# 步骤 2：基于 2015–2020 样本数据，测算各区县鸟类相对丰度指数
# 数据处理：微信公众号 RStata
# 复现 Liang, Rudik, Zou, Johnston, Rodewald & Kling (2020, PNAS)
#       "Conservation cobenefits from air pollution regulation: Evidence from birds"
# 的核心因变量：努力量校正后的「相对鸟类丰度」std(Γ̂)_cmy
#
# 此处代码需要下载讲义材料查看～

# 3. 合并 + 样本筛选
# ------------------------------------------------------------------------------
chk = checklist_counts.merge(
    rep_raw[["reportId", "date", "year", "month", "hour_day", "hours", "n_observers",
             "user_id", "site_id", "lon", "lat", "prov", "prov_code", "city",
             "city_code", "county", "county_code", "n_sp_report"]],
    on="reportId", how="inner")

print(chk.groupby("year").size())

# 更早年份的数据非常少，不值得使用和计算。
chk = chk[
    chk["county_code"].notna() & chk["year"].notna() &
    chk["year"].between(YEAR_MIN, YEAR_MAX) &
    chk["hours"].notna() & (chk["hours"] > 0) & (chk["hours"] <= HOURS_MAX) &
    (chk["n_ind_all"] >= 1)
].copy()

chk["cmy"] = (chk["county_code"].astype("int64").astype(str) + "_" +
              chk["year"].astype("int64").astype(str) + "_" +
              chk["month"].astype("int64").astype(str))

def season(m):
    if m in (3, 4, 5):
        return "spring"
    if m in (6, 7, 8):
        return "summer"
    if m in (9, 10, 11):
        return "autumn"
    return "winter"

chk["season"] = chk["month"].map(season)

print("sample:", len(chk), "checklists;",
      chk["county_code"].nunique(), "counties;",
      chk["cmy"].nunique(), "county-month-year units")

# 强制 stdout/stderr 走 ASCII：reticulate/knitr 捕获 Python 输出时硬编码 latin-1，
# 任何中文打印都会抛 UnicodeEncodeError 中断 knit。这里把输出中的非 ASCII 字符替换为
# '?'，作为兜底；正文讲解与 .py 源文件注释仍保留中文。
import sys as _sys
class _AsciiOut:
    def __init__(self, raw):
        self._raw = raw
    def write(self, s):
        if isinstance(s, str):
            s = s.encode("ascii", "replace").decode("ascii")
        return self._raw.write(s)
    def flush(self):
        return self._raw.flush()
    def __getattr__(self, n):
        return getattr(self._raw, n)
if not isinstance(_sys.stdout, _AsciiOut):
    _sys.stdout = _AsciiOut(_sys.stdout)
if not isinstance(_sys.stderr, _AsciiOut):
    _sys.stderr = _AsciiOut(_sys.stderr)

print("sample:", len(chk), "checklists;",
      chk["county_code"].nunique(), "counties;",
      chk["cmy"].nunique(), "county-month-year units")
print("1 SD individuals per checklist (paper reports 98.4):", round(float(chk["n_ind_all"].std()), 1),
      "; median:", float(chk["n_ind_all"].median()))

3.3 步骤 2（中）：Poisson 努力量校正（复现 Eq.1）

这是整个方法的核心。fit_gamma() 把「对一个被解释变量跑一次固定效应泊松回归、取出 cmy 固定效应估计值 Γ^」封装成可复用函数，再用循环批量生成全样本与四个分类群的指数。

# ------------------------------------------------------------------------------
# 4. 第一步：Poisson 努力量校正（复现论文 Eq.1）
# ------------------------------------------------------------------------------
# 原式:
#   #birds_{cohdmy} = exp(β_d·hours + β_n·n_observers + ζ_h + Γ_{cmy} + e)
# 用 pyfixest.fepois 把 ζ_h 与 Γ_{cmy} 都放进固定效应，
# 再用 fit.fixef() 取出 C(cmy) 这一组估计值（dict: cmy -> Γ̂）。
# Γ̂ 只被识别到相差一个常数，但下一步做 z 标准化后不影响结论。
# 强制 stdout/stderr 走 ASCII：reticulate/knitr 捕获 Python 输出时硬编码 latin-1，
# 任何中文打印都会抛 UnicodeEncodeError 中断 knit。这里把输出中的非 ASCII 字符替换为
# '?'，作为兜底；正文讲解与 .py 源文件注释仍保留中文。
# 此处代码需要下载讲义材料查看～

# 强制 stdout/stderr 走 ASCII：reticulate/knitr 捕获 Python 输出时硬编码 latin-1，
# 任何中文打印都会抛 UnicodeEncodeError 中断 knit。这里把输出中的非 ASCII 字符替换为
# '?'，作为兜底；正文讲解与 .py 源文件注释仍保留中文。
import sys as _sys
class _AsciiOut:
    def __init__(self, raw):
        self._raw = raw
    def write(self, s):
        if isinstance(s, str):
            s = s.encode("ascii", "replace").decode("ascii")
        return self._raw.write(s)
    def flush(self):
        return self._raw.flush()
    def __getattr__(self, n):
        return getattr(self._raw, n)
if not isinstance(_sys.stdout, _AsciiOut):
    _sys.stdout = _AsciiOut(_sys.stdout)
if not isinstance(_sys.stderr, _AsciiOut):
    _sys.stderr = _AsciiOut(_sys.stderr)

print("=== Eq.1 Poisson (all species) coefficients ===")
print(m_all)

3.4 步骤 2（下）：面板合并、中文标签与诊断

此处内容需要下载讲义材料查看～

3.5 步骤 2（图）：三张诊断图

# ------------------------------------------------------------------------------
# 8. 图 1：全国相对丰度趋势（校正 vs 未校正）
# ------------------------------------------------------------------------------
tr_rows = []
for year, g in panel[panel["gstd_all"].notna()].groupby("year"):
    z1 = (np.log(g["raw_ind_mean"]) - np.log(g["raw_ind_mean"]).mean()) / np.log(g["raw_ind_mean"]).std(ddof=1)
    z2 = (np.log(g["sp_richness"]) - np.log(g["sp_richness"]).mean()) / np.log(g["sp_richness"]).std(ddof=1)
    tr_rows.append({
        "year": year,
        "努力量校正后 std(Γ̂)": wmean(g["gstd_all"], g["n_checklist"]),
        "未校正: log(平均个体数)": wmean(z1, g["n_checklist"]),
        "未校正: log(累计物种数)": wmean(z2, g["n_checklist"]),
    })
tr = pd.DataFrame(tr_rows)
tr_long = tr.melt(id_vars="year", var_name="variable", value_name="value")

fig1, ax = plt.subplots(figsize=(12, 6))
colors = {"努力量校正后 std(Γ̂)": "#D7263D", "未校正: log(平均个体数)": "#1B98E0",
          "未校正: log(累计物种数)": "#4C9F70"}
for var, grp in tr_long.groupby("variable"):
    ax.plot(grp["year"], grp["value"], label=var, color=colors.get(var, None),
            linewidth=1, marker="o", markersize=2.2)
ax.axhline(0, color="grey", linewidth=0.3)
ax.set_xticks(range(2015, 2021))
ax.set_xlabel("")
ax.set_ylabel("标准化指数（z 分数）")
add_title_and_subtitle(
    ax,
    title="中国县级相对鸟类丰度：努力量校正前后对比（2015–2020）",
    subtitle="复现 Liang et al. (2020, PNAS) Eq.1 的 Poisson 努力量校正；按县×月×年报告数加权",
    title_fontsize=28, subtitle_fontsize=20)
ax.legend(frameon=False)
fig1.savefig(os.path.join(dir_out, "fig1_校正前后趋势_2015_2020.png"), dpi=300, bbox_inches="tight")
add_caption(ax, "数据来源：1980–2025 年观鸟记录（RStata 数据中心，2026 年 3 月爬取）；复现 Liang et al. (2020, PNAS)", fig=fig1)
plt.close(fig1)

# ------------------------------------------------------------------------------
# 9. 图 2：为什么必须做努力量校正
# ------------------------------------------------------------------------------
e1 = chk.copy()
e1["h"] = np.minimum(np.round(e1["hours"] * 2) / 2, 8)
e1 = e1.groupby("h")["n_ind_all"].agg(["mean", "size"]).reset_index()
e1 = e1[e1["size"] >= 200]
e2 = chk.groupby("hour_day")["n_ind_all"].mean().rename("m").reset_index().sort_values("hour_day")

fig2, axes = plt.subplots(1, 2, figsize=(15, 7.5))
axes[0].plot(e1["h"], e1["mean"], color="#D7263D", marker="o")
add_title_and_subtitle(axes[0], title="A  努力量：时长",
                       subtitle="单次观测越久，报告内平均个体数越多",
                       title_fontsize=26, subtitle_fontsize=18)
axes[0].set_xlabel("单次观鸟时长（小时）")
axes[0].set_ylabel("报告内平均个体数")
axes[1].bar(e2["hour_day"], e2["m"], color="#1B98E0")
axes[1].set_xticks(range(0, 24, 3))
add_title_and_subtitle(axes[1], title="B  可探测性：时刻",
                       subtitle="白天（尤其上午）更易观测到鸟，计数偏高",
                       title_fontsize=26, subtitle_fontsize=18)
axes[1].set_xlabel("开始观测时刻")
axes[1].set_ylabel("报告内平均个体数")
fig2.suptitle("原始观鸟计数被「观察者努力量」和「可探测性」严重污染", fontsize=28)
fig2.supxlabel("这正是 Liang et al. (2020) 先用 Poisson 模型剥离 hours 与 hour-of-day 的原因", fontsize=20, color="grey")
fig2.tight_layout(rect=[0, 0.04, 1, 0.90])
fig2.savefig(os.path.join(dir_out, "fig2_努力量与可探测性_2015_2020.png"), dpi=300, bbox_inches="tight")
add_caption(axes[0], "数据来源：1980–2025 年观鸟记录（RStata 数据中心，2026 年 3 月爬取）；复现 Liang et al. (2020, PNAS)", fig=fig2)
plt.close(fig2)

# ------------------------------------------------------------------------------
# 10. 图 3：论文口径的分类群指数
# ------------------------------------------------------------------------------
# 此处代码需要下载讲义材料查看～

三张诊断图（图 1–3）

图 1 展示全国县级相对鸟类丰度在努力量校正前后的对比：

图 2 从「努力量（时长）」和「可探测性（观测时刻）」两个角度说明为什么必须做努力量校正：

图 3 展示论文口径的四个分类群（水禽 / 涉禽 / 水鸟 / 陆鸟）相对丰度指数：

3.6 步骤 3：与刘钊等（2025）论文趋势对照

刘钊等 (2025, 经济学（季刊）) 以「绿色金融改革创新试验区（2017）」为准自然实验，核心结论是试验区鸟类丰富度显著提升；其被解释变量明确「借鉴 Liang et al. (2020) 计算得到并标准化」。下面用我们产出的面板做方向性（naive DiD）检验，看趋势是否一致。

# 对照刘钊等 (2025, 经济学季刊) 的趋势一致性检验
# 数据处理：微信公众号 RStata
# 论文：绿色金融改革创新试验区(2017)显著提升鸟类丰富度(借鉴 Liang et al.2020 方法)
# 本脚本用「02_测算鸟类丰度指数_2015_2020.py」产出的县×月×年面板（样本 2015–2020），
# 做方向性(naive DiD)与生态学合理性检验。
# 强制 stdout/stderr 走 ASCII：reticulate/knitr 捕获 Python 输出时硬编码 latin-1，
# 任何中文打印都会抛 UnicodeEncodeError 中断 knit。这里把输出中的非 ASCII 字符替换为
# '?'，作为兜底；正文讲解与 .py 源文件注释仍保留中文。
import sys as _sys
class _AsciiOut:
    def __init__(self, raw):
        self._raw = raw
    def write(self, s):
        if isinstance(s, str):
            s = s.encode("ascii", "replace").decode("ascii")
        return self._raw.write(s)
    def flush(self):
        return self._raw.flush()
    def __getattr__(self, n):
        return getattr(self._raw, n)
if not isinstance(_sys.stdout, _AsciiOut):
    _sys.stdout = _AsciiOut(_sys.stdout)
if not isinstance(_sys.stderr, _AsciiOut):
    _sys.stderr = _AsciiOut(_sys.stderr)

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# 统一绘图样式（中文字体 + 标题/副标题/数据来源说明），参考 python-china-map skill
_sys.path.insert(0, os.getcwd())
from plot_style import setup_chinese_font, add_title_and_subtitle, add_caption
setup_chinese_font()

def wmean(x, w):
    """加权平均（对应 R 的 weighted.mean，na.rm=TRUE 会先剔除缺失）"""
    x = np.asarray(x, dtype=float)
    w = np.asarray(w, dtype=float)
    mask = ~np.isnan(x)
    if mask.sum() == 0:
        return float("nan")
    return float(np.average(x[mask], weights=w[mask]))

dir_in = "/Users/ac/Desktop/使用 Python 测算各区县鸟类丰度指数"
dir_out = os.path.join(dir_in, "output")

panel = pd.read_parquet(os.path.join(dir_out, "panel_county_month_year_2015_2020.parquet"))
print("panel:", len(panel), "county-month-year units; years", int(panel["year"].min()), "-", int(panel["year"].max()))

## ---- 1. 试验区定义（首批准试验区，2017年6月）----
## 赣江新区→南昌/九江；贵安新区→贵阳；花都→广州。整市视为处理组（偏保守，稀释效应）。
pilot_cities = ["湖州市", "衢州市", "南昌市", "九江市",
                "广州市", "贵阳市", "兰州市", "哈密市"]
panel["treat"] = panel["city"].isin(pilot_cities)
panel["post"] = panel["year"] >= 2017

print("\nTreated-city units:", int(panel["treat"].sum()), " Control units:", int((~panel["treat"]).sum()))

## ---- 2. 朴素双重差分 (2015-2020, 政策前=2015-2016, 后=2017-2020) ----
dd = panel[panel["year"].between(2015, 2020)].copy()
dd["period"] = np.where(dd["year"] <= 2016, "pre", "post")
tab = dd.groupby(["treat", "period"]).apply(
    lambda g: pd.Series({"g": wmean(g["gstd_all"], g["n_checklist"])}), include_groups=False
).reset_index()
print(tab.to_string(index=False))

# 2×2 -> DiD = (post-pilot - pre-pilot) - (post-non - pre-non)
wide_tab = tab.pivot(index="treat", columns="period", values="g").reset_index()
g_post_pilot = wide_tab.loc[wide_tab["treat"], "post"].iloc[0]
g_pre_pilot = wide_tab.loc[wide_tab["treat"], "pre"].iloc[0]
g_post_non = wide_tab.loc[~wide_tab["treat"], "post"].iloc[0]
g_pre_non = wide_tab.loc[~wide_tab["treat"], "pre"].iloc[0]
naive_did = (g_post_pilot - g_pre_pilot) - (g_post_non - g_pre_non)
print(f"\nnaive DiD (abundance index) = {naive_did:.4f}")
print("Paper Table 1 DID coef beta1 is significantly positive -> if our index direction matches, it should be > 0.")
# 这个数据印证结论：试验区设立后鸟类相对丰度更高。

# 物种丰富度口径的朴素 DiD (对应论文表3中的“鸟种丰富度")
dd_sp = panel[panel["year"].between(2015, 2020)].copy()
dd_sp["period"] = np.where(dd_sp["year"] <= 2016, "pre", "post")
tab_sp = dd_sp.groupby(["treat", "period"]).apply(
    lambda g: pd.Series({"s": wmean(g["sp_richness"], g["n_checklist"])}), include_groups=False
).reset_index()
print("\nSpecies-richness 2x2 means:")
print(tab_sp.to_string(index=False))
wide_sp = tab_sp.pivot(index="treat", columns="period", values="s").reset_index()
s_post_pilot = wide_sp.loc[wide_sp["treat"], "post"].iloc[0]
s_pre_pilot = wide_sp.loc[wide_sp["treat"], "pre"].iloc[0]
s_post_non = wide_sp.loc[~wide_sp["treat"], "post"].iloc[0]
s_pre_non = wide_sp.loc[~wide_sp["treat"], "pre"].iloc[0]
naive_did_sp = (s_post_pilot - s_pre_pilot) - (s_post_non - s_pre_non)
print(f"\nnaive DiD (species richness) = {naive_did_sp:.4f}")
print("Note: this species count is raw (uncorrected); the paper usually also applies Liang correction for heterogeneity.")

## ---- 3. 事件研究式年度轨迹 (2015-2020) ----
# 此处代码需要下载讲义材料查看～

print("\n=== descriptive trend summary ===")
print("National annual raw abundance mean (gamma_hat, log scale):")
yr_rows = []
for year, g in panel.groupby("year"):
    yr_rows.append({"year": year, "raw": wmean(g["gamma_hat"], g["n_checklist"])})
print(pd.DataFrame(yr_rows).to_string(index=False))
print("\nGroup 2015 vs 2020 change (standardized index):")
cmp = panel[panel["year"].isin([2015, 2020])].groupby("year").apply(
    lambda g: pd.Series({
        "gstd_all": round(wmean(g["gstd_all"], g["n_checklist"]), 4),
        "gstd_waterbird": round(wmean(g["gstd_waterbird"], g["n_checklist"]), 4),
        "gstd_landbird": round(wmean(g["gstd_landbird"], g["n_checklist"]), 4),
    }), include_groups=False
).reset_index()
print(cmp.to_string(index=False))

与刘钊等（2025）论文趋势对照图

下图把我们的指数与刘钊等 (2025) 论文做趋势一致性对照（年度轨迹、季节模式、分类群趋势、区域基线）：

# 强制 stdout/stderr 走 ASCII：reticulate/knitr 捕获 Python 输出时硬编码 latin-1，
# 任何中文打印都会抛 UnicodeEncodeError 中断 knit。这里把输出中的非 ASCII 字符替换为
# '?'，作为兜底；正文讲解与 .py 源文件注释仍保留中文。
import sys as _sys
class _AsciiOut:
    def __init__(self, raw):
        self._raw = raw
    def write(self, s):
        if isinstance(s, str):
            s = s.encode("ascii", "replace").decode("ascii")
        return self._raw.write(s)
    def flush(self):
        return self._raw.flush()
    def __getattr__(self, n):
        return getattr(self._raw, n)
if not isinstance(_sys.stdout, _AsciiOut):
    _sys.stdout = _AsciiOut(_sys.stdout)
if not isinstance(_sys.stderr, _AsciiOut):
    _sys.stderr = _AsciiOut(_sys.stderr)

print(tab.to_string(index=False))
print(f"naive DiD (abundance index) = {naive_did:.4f}")
print(f"naive DiD (species richness) = {naive_did_sp:.4f}")

解读：丰度指数口径的朴素 DiD = r round(py$naive_did, 4) > 0，与刘钊等 (2025) 表 1 中 DID 系数 β1 显著为正的方向一致；而未校正物种数口径的 DiD 为负，是因为非试验区观测基数增长更快——这正说明必须用努力量校正后的指数做比较，不能拿原始物种数直接比。

四、小结

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Python 测算各区县鸟类丰度指数