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

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

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

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

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

一、指标来源与计算过程

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

一句话：论文要的是「相对」丰度，不是「绝对」有多少只鸟——它只能告诉你 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 之前）数据极稀疏，县×月×年固定效应难以识别，故不纳入本示例窗口。

二、详细讲解计算代码

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

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

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

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

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

suppressMessages({ library(tidyverse); library(haven) })

dir_in <- "/Users/ac/Desktop/计算地区生物多样性/1980~2025年观鸟记录、经纬度及其所处的省市区县数据(2026年3月爬取)"
dir_out <- "/Users/ac/Desktop/使用 R 语言测算各区县鸟类丰度指数"
dir.create(dir_out, showWarnings = FALSE, recursive = TRUE)

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

YEAR_MIN <- 2015
YEAR_MAX <- 2020

# ------------------------------------------------------------------------------
# 1) 观鸟记录表：按年份筛选
# ------------------------------------------------------------------------------
cat("== 步骤1：读取观鸟记录表并按年份筛选 ==\n")
rep_all <- read_dta(f_rep) |>
  mutate(year = lubridate::year(观测时间_起始)) |>
  filter(year %in% YEAR_MIN:YEAR_MAX)

cat("筛选后报告数:", nrow(rep_all), "\n年份分布:\n")
print(rep_all |> count(year) |> arrange(year))

out_rep <- file.path(dir_out, "观鸟记录_2015_2020.dta")
write_dta(rep_all, out_rep)
cat("已写出 ->", out_rep, "\n")

# ------------------------------------------------------------------------------
# 2) 鸟种观测统计报告：按 reportId 关联筛选 （2015–2020 的报告才保留）
# ------------------------------------------------------------------------------
cat("\n== 步骤2：按 reportId 关联筛选鸟种观测统计报告 ==\n")
rep_ids <- rep_all |> distinct(reportId)
cat("2015–2020 唯一 reportId 数:", nrow(rep_ids), "\n")

sp_all <- read_dta(f_sp) |>
  inner_join(rep_ids, by = "reportId")

cat("鸟种报告筛选后行数:", nrow(sp_all), "\n")

out_sp <- file.path(dir_out, "鸟种观测统计报告_2015_2020.dta")
write_dta(sp_all, out_sp)
cat("已写出 ->", out_sp, "\n")

cat("\n样本筛选完成。两个样本 dta 已保存至:", dir_out, "\n")

筛选结果如下：

cat("筛选后报告数:", nrow(rep_all), "\n")
cat("鸟种报告筛选后行数:", nrow(sp_all), "\n")
print(rep_all |> count(year) |> arrange(year))

2.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
# 此处代码需要下载讲义材料查看～

# ------------------------------------------------------------------------------
# 0. 两份数据的角色对照（沿用复现脚本口径）
# ------------------------------------------------------------------------------
# 观鸟记录表(report级, 约 125 万行)  ≈ 论文里的 "checklist"
#   reportId          -> checklist 唯一 id
#   观测时间_起始/结束 -> hours（努力量）与 hour-of-day（可探测性）
#   县代码             -> 论文的 county c
#   鸟种数量           -> 该报告里的【物种数】(richness)
#   用户编号/地点编号   -> 可构造观察者、地点层面的努力量代理
#
# 鸟种观测统计报告(report×物种级, 约 2000 万行) ≈ eBird Reference Dataset
#   鸟种数量           -> 该物种在该报告中的【个体数】
#   鸟种目             -> 划分论文的 waterfowl / shorebirds / waterbirds / landbirds
#
# 因此论文 Eq.1 左边的 "#birds observed" = 明细表按 reportId 求和的个体数。

# 论文分类群划分（Rosenberg et al. 2019 口径，映射到中文目名）
grp_waterfowl  <- c("雁形目")
grp_shorebirds <- c("鸻形目")
grp_waterbirds <- c("鹈形目", "鹤形目", "䴙䴘目", "鲣鸟目", "鹳形目",
                    "潜鸟目", "鹱形目", "鹲形目", "红鹳目")
# 其余目 = 陆鸟 landbirds

# ------------------------------------------------------------------------------
# 1. 物种明细 -> 报告级计数（总量 + 分类群）
# ------------------------------------------------------------------------------
sp_raw <- read_dta(f_sp) |>
  rename(sp_id = 鸟种编号, order = 鸟种目, n_ind = 鸟种数量)

unique(sp_raw$order)
# 也可以让豆包对这些进行分类，结果也是一样的。

# 个体数清洗
# 此处代码需要下载讲义材料查看～

# ------------------------------------------------------------------------------
# 2. 报告表：努力量 / 可探测性 / 地理时间
# ------------------------------------------------------------------------------
rep_raw <- read_dta(f_rep) |>
  distinct(reportId, .keep_all = TRUE) |>
  rename(user_id = 用户编号, site_id = 地点编号, lon = 经度_wgs84, lat = 纬度_wgs84,
         prov = 省, prov_code = 省代码, city = 市, city_code = 市代码,
         county = 县, county_code = 县代码, n_sp_report = 鸟种数量) |>
  mutate(
    date     = as.Date(观测时间_起始),
    year     = lubridate::year(观测时间_起始),
    month    = lubridate::month(观测时间_起始),
    hour_day = lubridate::hour(观测时间_起始),
    hours    = as.numeric(difftime(观测时间_结束, 观测时间_起始, units = "hours"))
  ) |>
  group_by(site_id, date) |>
  # 论文有 number of observers；本数据没有，用同一地点同一天提交的不同用户数近似团体规模
  mutate(n_observers = n_distinct(user_id)) |>
  ungroup()

rep_raw

# ------------------------------------------------------------------------------
# 3. 合并 + 样本筛选
# ------------------------------------------------------------------------------
# 此处代码需要下载讲义材料查看～

cat("样本:", nrow(chk), "份观鸟报告;",
    n_distinct(chk$county_code), "个县;",
    n_distinct(chk$cmy), "个县×月×年单元\n")
cat("1 SD 每份报告个体数 (论文报告 98.4):", round(sd(chk$n_ind_all), 1),
    "; 中位数:", median(chk$n_ind_all), "\n")

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

# ------------------------------------------------------------------------------
# 4. 第一步：Poisson 努力量校正（复现论文 Eq.1）
# ------------------------------------------------------------------------------
# 此处代码需要下载讲义材料查看～

# gamma_hat 对应论文 Eq.1 的 Γ̂_cmy；
# gamma_std 对应论文 Eq.2/Eq.3 拿去和 O3、PM2.5 回归的 std(Γ̂)_cmy——它就是论文因变量的现货版本；
# gamma_hat = 0 那行（第 1 行）也是正常的。固定效应是以某个参照为基准估计的，部分单元的估计值恰好落在 0 附近（显示精度截断成 0）。它不代表"没有鸟"，只是"该单元丰度等于模型暗含的基准水平"。

cat("=== Eq.1 Poisson (全部鸟种) 系数 ===\n")
print(summary(m_all))

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

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

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

# ------------------------------------------------------------------------------
# 8. 图 1：全国相对丰度趋势（校正 vs 未校正）
# ------------------------------------------------------------------------------
tr <- panel |>
  filter(!is.na(gstd_all)) |>
  group_by(year) |>
  summarise(
    `努力量校正后 std(Γ̂)` = weighted.mean(gstd_all, n_checklist),
    `未校正: log(平均个体数)` = weighted.mean(as.numeric(scale(log(raw_ind_mean))), n_checklist),
    `未校正: log(累计物种数)` = weighted.mean(as.numeric(scale(log(sp_richness))), n_checklist)
  ) |>
  pivot_longer(-year, names_to = "variable", values_to = "value") |>
  arrange(year)

p1 <- ggplot(tr, aes(year, value, colour = variable)) +
  geom_hline(yintercept = 0, linewidth = .3, colour = "grey60") +
  geom_line(linewidth = 1) + geom_point(size = 2.2) +
  scale_x_continuous(breaks = 2015:2020) +
  scale_colour_manual(values = c("#D7263D", "#1B98E0", "#4C9F70")) +
  labs(x = NULL, y = "标准化指数（z 分数）", colour = NULL,
       title = "中国县级相对鸟类丰度：努力量校正前后对比（2015–2020）",
       subtitle = "复现 Liang et al. (2020, PNAS) Eq.1 的 Poisson 努力量校正；按县×月×年报告数加权",
       caption = "数据：RStata 数据中心「1980–2025 年观鸟记录」（样本 2015–2020），数据处理：微信公众号 RStata")
ggsave(file.path(dir_out, "fig1_校正前后趋势_2015_2020.png"), p1,
       width = 12, height = 6, dpi = 300, bg = "white", device = png)

# ------------------------------------------------------------------------------
# 9. 图 2：为什么必须做努力量校正
# ------------------------------------------------------------------------------
# 此处代码需要下载讲义材料查看～

# ------------------------------------------------------------------------------
# 10. 图 3：论文口径的分类群指数
# ------------------------------------------------------------------------------
tg <- panel |>
  select(year, n_checklist, gstd_waterfowl, gstd_shorebird, gstd_waterbird, gstd_landbird) |>
  pivot_longer(-c(year, n_checklist), names_to = "variable", values_to = "value") |>
  filter(!is.na(value)) |>
  group_by(year, variable) |>
  summarise(v = weighted.mean(value, n_checklist), .groups = "drop") |>
  mutate(variable = factor(variable,
    levels = c("gstd_landbird", "gstd_waterfowl", "gstd_shorebird", "gstd_waterbird"),
    labels = c("陆鸟 landbirds", "雁鸭类 waterfowl", "鸻鹬鸥类 shorebirds", "其他水鸟 waterbirds")))

p3 <- ggplot(tg, aes(year, v, colour = variable)) +
  geom_hline(yintercept = 0, linewidth = .3, colour = "grey60") +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  scale_color_brewer(palette = "Set2") +
  scale_x_continuous(breaks = 2015:2020) +
  labs(x = NULL, y = "std(Γ̂)", colour = NULL,
       title = "论文口径的四个分类群相对丰度指数（2015–2020，可用本数据直接复现）",
       subtitle = "类群按「鸟种目」映射 Rosenberg et al. (2019) 的 waterfowl / shorebirds / waterbirds / landbirds")
ggsave(file.path(dir_out, "fig3_分类群指数_2015_2020.png"), p3,
       width = 12, height = 6, dpi = 300, bg = "white", device = png)

三张诊断图（图 1–3）

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

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

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

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

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

# 对照刘钊等 (2025, 经济学季刊) 的趋势一致性检验
# 数据处理：微信公众号 RStata
# 论文：绿色金融改革创新试验区(2017)显著提升鸟类丰富度(借鉴 Liang et al.2020 方法)
# 本脚本用「02_测算鸟类丰度指数_2015_2020.R」产出的县×月×年面板（样本 2015–2020），
# 做方向性(naive DiD)与生态学合理性检验。
suppressMessages({
  library(arrow); library(tidyverse); library(ggplot2); library(patchwork)
})

dir_in  <- "/Users/ac/Desktop/使用 R 语言测算各区县鸟类丰度指数"
dir_out <- file.path(dir_in, "output")

panel <- read_parquet(file.path(dir_out, "panel_county_month_year_2015_2020.parquet"))
cat("面板:", nrow(panel), "县×月×年单元; 年份", min(panel$year), "-", max(panel$year), "\n")

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

cat("\n处理组城市(县×月×年单元数):",
    sum(panel$treat), " 对照组:", sum(!panel$treat), "\n")

## ---- 2. 朴素双重差分 (2015-2020, 政策前=2015-2016, 后=2017-2020) ----
dd <- panel %>% filter(year %in% 2015:2020) %>%
  mutate(period = if_else(year <= 2016, "pre", "post"))
tab <- dd %>% group_by(treat, period) %>%
  summarise(g = weighted.mean(gstd_all, n_checklist, na.rm = TRUE), .groups = "drop")
print(tab)
did <- tab %>% pivot_wider(names_from = period, values_from = g) %>%
  mutate(diff = post - pre) %>% pull(diff)
naive_did <- did[2] - did[1]   # (post-pilot - pre-pilot) - (post-non - pre-non)
cat(sprintf("\n朴素 DiD (丰度指数) = %.4f\n", naive_did))
cat("论文表1 DID系数 β1 显著为正 → 若本指数方向一致，应 >0。\n")

# 这个数据印证结论：试验区设立后鸟类相对丰度更高。

# 物种丰富度口径的朴素 DiD (对应论文表3中的“鸟种丰富度")
dd_sp <- panel %>% filter(year %in% 2015:2020) %>%
  mutate(period = if_else(year <= 2016, "pre", "post"))
tab_sp <- dd_sp %>% group_by(treat, period) %>%
  summarise(s = weighted.mean(sp_richness, n_checklist, na.rm = TRUE), .groups = "drop")
cat("\n物种数(richness) 2×2 均值:\n"); print(tab_sp)
did_sp <- tab_sp %>% pivot_wider(names_from = period, values_from = s) %>%
  mutate(diff = post - pre) %>% pull(diff)
naive_did_sp <- did_sp[2] - did_sp[1]
cat(sprintf("\n朴素 DiD (物种数 richness) = %.4f\n", naive_did_sp))
cat("注意：本物种数为未做努力量校正原始计数，论文种类异质性一般也做 Liang 校正。\n")

# 此处代码需要下载讲义材料查看～

## ---- 6. 区域基线(描述性上下文) ----
region_map <- c(
  "北京市"="东部","天津市"="东部","河北省"="东部","上海市"="东部","江苏省"="东部",
  "浙江省"="东部","福建省"="东部","山东省"="东部","广东省"="东部","海南省"="东部",
  "山西省"="中部","安徽省"="中部","江西省"="中部","河南省"="中部","湖北省"="中部",
  "湖南省"="中部","广西壮族自治区"="中部",
  "内蒙古自治区"="西部","重庆市"="西部","四川省"="西部","贵州省"="西部","云南省"="西部",
  "西藏自治区"="西部","陕西省"="西部","甘肃省"="西部","青海省"="西部","宁夏回族自治区"="西部",
  "新疆维吾尔自治区"="西部",
  "辽宁省"="东北","吉林省"="东北","黑龙江省"="东北",
  "台湾省"="其他","香港特别行政区"="其他","澳门特别行政区"="其他")
panel$region <- recode(panel$prov, !!!region_map)
reg <- panel %>% group_by(region) %>%
  summarise(g = weighted.mean(gstd_all, n_checklist, na.rm = TRUE), .groups = "drop") %>%
  filter(region %in% c("东部","中部","西部","东北"))
p_reg <- ggplot(reg, aes(x = region, y = g, fill = region)) +
  geom_col() + scale_fill_brewer(palette = "Set2") +
  labs(x = "区域", y = "相对丰度指数均值", title = "图D 区域基线丰度(描述性)")

## ---- 输出 ----
fig <- (p_es + p_seas) / (p_tg + p_reg)
ggsave(file.path(dir_out, "fig4_与刘钊论文对照_2015_2020.png"), fig,
       width = 20, height = 12, dpi = 300, device = png)
cat("\n已保存 output/fig4_与刘钊论文对照_2015_2020.png\n")

cat("\n=== 描述性趋势汇总 ===\n")
cat("全国年度 原始丰度均值(gamma_hat, 对数尺度):\n")
print(panel %>% group_by(year) %>%
      summarise(raw = weighted.mean(gamma_hat, n_checklist, na.rm = TRUE)) %>%
      arrange(year), n = 20)
cat("\n分类群 2015 vs 2020 变化(标准化指数):\n")
print(panel %>% filter(year %in% c(2015, 2020)) %>%
      group_by(year) %>%
      summarise(across(c(gstd_all, gstd_waterbird, gstd_landbird),
                       ~ round(weighted.mean(., n_checklist, na.rm = TRUE), 4)), .groups="drop"))

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

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

print(tab)
cat(sprintf("朴素 DiD (丰度指数) = %.4f\n", naive_did))
cat(sprintf("朴素 DiD (物种数 richness) = %.4f\n", naive_did_sp))

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

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