17  数据清洗

18 数据清洗

真实世界的数据从来都不是”干净”的。缺失值、异常值、格式不一致、重复记录……这些问题在生态学数据中尤为常见。数据清洗是数据分析的第一步,也是最重要的一步–垃圾进,垃圾出(Garbage In, Garbage Out)。

library(tidyverse)
library(lubridate)
library(naniar)

18.1 示例数据

我们使用一个基于真实研究参数模拟的土壤调查数据集来演示数据清洗的各个步骤。

18.1.1 数据来源

数据来源:Li, D., Wen, L., Yang, L., Luo, P., Xiao, W., Chen, H., … & Wang, K. (2017). Dynamics of soil organic carbon and nitrogen following agricultural abandonment in a karst region. Journal of Geophysical Research: Biogeosciences, 122(11), 2849-2864. DOI: 10.1002/2016JG003683

研究背景:在中国西南喀斯特地区,利用”空间换时间”方法,研究了白云岩和石灰岩两种岩性下弃耕地演替过程中土壤碳氮含量的变化。共125个样地,涵盖耕地、草地、灌丛和次生林四种土地利用类型。

数据集基本信息:

项目	内容
研究区域	中国西南喀斯特地区(广西、贵州)
样地数量	125个样地
岩性类型	白云岩 / 石灰岩
土地利用	耕地、草地、灌丛、次生林
土层深度	0-10 cm, 10-20 cm, 20-30 cm
测定指标	土壤有机碳(SOC)、全氮(TN)、pH、可交换钙

18.1.2 模拟数据生成

以下数据基于该研究的实际参数范围模拟,用于演示数据清洗流程:

set.seed(2027)

# 模拟125个样地的土壤数据
# 基于Li et al. (2017)研究的实际参数范围

n_sites <- 125

# 土地利用类型和岩性组合
land_use <- rep(c("耕地", "草地", "灌丛", "次生林"), length.out = n_sites)
lithology <- rep(c("白云岩", "石灰岩"), length.out = n_sites)

# 根据土地利用类型设置SOC (g/kg)的基础值
# Li et al. (2017)发现:耕地最低,森林最高
soc_base <- case_when(
  land_use == "耕地" ~ 8.5,
  land_use == "草地" ~ 12.3,
  land_use == "灌丛" ~ 18.7,
  land_use == "次生林" ~ 25.4
)

# 根据岩性调整(石灰岩地区SOC更高)
soc_base <- ifelse(lithology == "石灰岩", soc_base * 1.15, soc_base)

# 生成SOC数据,添加随机变异和异常值
soc <- soc_base + rnorm(n_sites, 0, 4)
# 添加异常值(一些不合理的高值和负值)
soc[c(15, 45, 78)] <- c(150, -3.5, 85)
# 添加缺失值
soc[c(23, 67, 102)] <- NA

# TN (g/kg) - 与SOC正相关
tn_base <- soc_base / 10
tn <- tn_base + rnorm(n_sites, 0, 0.3)
tn[c(34, 89)] <- c(-0.5, 25)  # 异常值
tn[c(12, 56, 91)] <- NA  # 缺失值

# pH值
ph <- 5.5 + rnorm(n_sites, 0, 0.5)
ph[8] <- -0.8  # 不可能的负值
ph[67] <- 15.2  # 不可能的极高值
ph[c(3, 45, 88)] <- NA

# 可交换钙 (cmol/kg)
ca <- (soc_base / 3) + rnorm(n_sites, 0, 1.5)
ca[c(7, 34)] <- c(-2, 120)
ca[c(15, 78)] <- NA

# 土层深度
depth <- rep(c("0-10 cm", "10-20 cm", "20-30 cm"), length.out = n_sites)

# 采样日期(故意设置格式不一致)
date_raw <- c(
  rep("2024-03-15", 40),
  rep("2024/03/16", 35),
  rep("15-Mar-2024", 30),
  rep("2024-03-17", 20)
)

# 组合数据框
soil_raw <- tibble(
  site_id = paste0("K", sprintf("%03d", 1:n_sites)),
  land_use = land_use,
  lithology = lithology,
  depth = depth,
  soc = round(soc, 2),        # 土壤有机碳 g/kg
  tn = round(tn, 3),          # 全氮 g/kg
  ph = round(ph, 2),
  ca = round(ca, 2),         # 可交换钙 cmol/kg
  date_raw = date_raw[1:n_sites]
)

# 添加一些重复行
soil_raw <- bind_rows(soil_raw, soil_raw[c(1, 2, 2, 15), ])

cat("数据维度:", dim(soil_raw), "\n")
cat("原始数据前10行:\n")
head(soil_raw, 10)

18.2 第一步:了解你的数据

在动手清洗之前,先全面了解数据的基本情况:

# 数据维度
dim(soil_raw)

# 数据结构
str(soil_raw)

# 描述性统计
summary(soil_raw)

# 每列缺失值数量
colSums(is.na(soil_raw))

Note数据字典

数据字典是描述数据集中每个变量含义、取值范围和数据类型的元数据文档。在生态学研究中，一个完善的数据字典是数据共享和重复利用的基础，也是论文方法学部分的重要内容。

为什么要建立数据字典？ 当你把数据分享给合作者或将数据上传到公共数据库（如 GBIF）时，数据字典帮助他人理解每个变量的含义。没有数据字典，别人拿到数据时只能靠猜，不仅效率低下，还容易产生误解。

数据字典的核心要素

要素	说明	示例
变量名	字段的英文或拼音缩写	site_id, soc, tn
含义	变量的中文描述	土壤有机碳含量
单位	测量单位	g/kg, cmol/kg
合理范围	正常取值范围	pH: 3.5-9.5
数据类型	数值/分类/日期	numeric, character
缺失值标记	NA 如何定义	-999 表示未测定

在 R 中管理数据字典

library(tidyverse)
library(readxl)

# 读取数据字典
dict <- read_excel("data/data_dictionary.xlsx")

# 查看所有变量
dict

# 检验数据集是否符合数据字典
expert_ranges <- list(
  ph = c(3.5, 9.5),
  soc = c(0, 80),
  tn = c(0, 10),
  ca = c(0, 60)
)

# 检测超出合理范围的异常值
for (var_name in names(expert_ranges)) {
  range_vec <- expert_ranges[[var_name]]
  outliers <- soil |>
    filter(.data[[var_name]] < range_vec[1] | .data[[var_name]] > range_vec[2]) |>
    select(site_id, all_of(var_name))
  if (nrow(outliers) > 0) {
    cat(var_name, "有", nrow(outliers), "个值超出范围\n")
    print(outliers)
  }
}

生态学案例

某团队在提交数据到国家生态科学数据中心时，审核人员指出数据字典缺失，导致数据被退回。团队补充了完整的数据字典后，数据被接受并公开。现在该数据集已被引用 23 次，成为同类研究的基准数据。这个案例说明：数据字典不是可有可无的文档，而是数据共享的”通行证”。

扩展记录：2026-04-09 | 目标字数：800+

18.3 处理重复值

# 检查重复行
sum(duplicated(soil_raw))

# 查看重复的行
soil_raw[duplicated(soil_raw), ]

# 去除重复行
soil <- soil_raw |>
  distinct()

cat("去重前:", nrow(soil_raw), "行\n")
cat("去重后:", nrow(soil), "行\n")

18.4 处理异常值

Important为什么先处理异常值？

在数据清洗流程中，处理异常值和缺失值的顺序至关重要。许多初学者会直觉地先填充缺失值，但这是一个常见的错误。正确的顺序应该是：先移除/修正异常值，再填充缺失值。这个原则背后有深刻的统计学逻辑。

异常值对统计量的“污染”效应：

异常值会严重扰乱数据的中心趋势和离散程度。具体来说：

1. 均值偏移：一个极端异常值可以将均值拉向不合理的方向。例如，如果土壤pH值数据中有一个错误记录为150（正常范围为4-9），计算出的均值将远高于真实水平。

2. 标准差虚高：异常值会大幅增加数据的方差和标准差，让数据看起来比实际更加离散。这会影响后续的统计检验（如t检验、方差分析）的结果。

3. 相关性失真：在多变量分析中，异常值会扭曲变量间的相关关系，导致错误的回归系数或相关系数。

错误顺序的后果：

如果先填充缺失值，再处理异常值，会发生以下问题：

# 错误示例：先填充再处理异常值
soil_wrong <- soil |>
  mutate(soc = ifelse(is.na(soc), mean(soc, na.rm = TRUE), soc))  # 先填充
# 此时的均值已经被异常值（如150）“污染”
# 填充值不是真实的中心值，而是偏高的值

# 正确示例：先处理异常值再填充
soil_correct <- soil |>
  filter(soc > 0 & soc < 100) |>  # 先移除异常值
  mutate(soc = ifelse(is.na(soc), mean(soc, na.rm = TRUE), soc))  # 再填充
# 此时的均值是基于正常数据计算的，更接近真实值

生态学案例：

在一项植物群落调查中，某样地的物种丰富度被错误记录为999（实际应为9）。如果先用均值填充其他缺失值，计算出的均值为85种（包含异常值999），远高于真实的平均水平12种。这样填充的缺失值就会是85种，完全不合理。而如果先移除异常值999，再计算均值，得到的12种，用这个值填充缺失值才是合理的。

实际操作流程：

1. 第一步：识别异常值（箱线图、IQR方法、专业知识判断）
2. 第二步：处理异常值（删除、替换为缺失值、或修正）
3. 第三步：识别缺失值模式（随机缺失、系统性缺失）
4. 第四步：填充缺失值（均值、中位数、回归预测等）

例外情况：

如果你使用的是中位数填充而非均值填充，异常值对中位数的影响较小，顺序可以灵活调整。但为了保持流程的一致性和可解释性，仍然建议先处理异常值。

扩展记录：2026-04-09 | 目标字数：800+

18.4.1 识别异常值

识别异常值是数据清洗的第一步。在正式处理之前，首先要发现哪些值可能是异常的。常用的方法包括：基于统计学的 IQR 方法和 Z-score 方法，基于专业知识的合理范围判断，以及可视化方法。

箱线图法

箱线图是识别单变量异常值的直观工具。在 R 中使用 ggplot2 可以快速生成：

library(tidyverse)

# 箱线图识别异常值
soil |>
  select(where(is.numeric)) |>
  pivot_longer(everything(), names_to = "variable", values_to = "value") |>
  ggplot(aes(x = variable, y = value)) +
  geom_boxplot(fill = "lightblue", outlier.color = "red", outlier.shape = 1) +
  facet_wrap(~variable, scales = "free") +
  labs(title = "数值变量箱线图（红点为潜在异常值）") +
  theme_minimal()

箱线图中，超出箱体上下边缘（Q1-1.5xIQR 和 Q3+1.5xIQR）的点会被标记为异常值（红点）。

Z-score 法

Z-score 表示某个值偏离均值的标准差倍数。通常 |Z| > 3 的值被视为潜在异常值：

# Z-score 方法检测异常值
detect_outliers_zscore <- function(x, threshold = 3) {
  z_scores <- abs((x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE))
  return(z_scores > threshold)
}

# 对所有数值变量应用
soil_marked <- soil |>
  mutate(
    soc_z = abs((soc - mean(soc, na.rm = TRUE)) / sd(soc, na.rm = TRUE)),
    tn_z = abs((tn - mean(tn, na.rm = TRUE)) / sd(tn, na.rm = TRUE)),
    ph_z = abs((ph - mean(ph, na.rm = TRUE)) / sd(ph, na.rm = TRUE)),
    ca_z = abs((ca - mean(ca, na.rm = TRUE)) / sd(ca, na.rm = TRUE)),
    soc_outlier = soc_z > 3,
    tn_outlier = tn_z > 3,
    ph_outlier = ph_z > 3,
    ca_outlier = ca_z > 3
  )

# 查看异常值汇总
outlier_summary <- soil_marked |>
  summarise(
    soc_outliers = sum(soc_outlier, na.rm = TRUE),
    tn_outliers = sum(tn_outlier, na.rm = TRUE),
    ph_outliers = sum(ph_outlier, na.rm = TRUE),
    ca_outliers = sum(ca_outlier, na.rm = TRUE)
  )
print(outlier_summary)

基于专业知识的合理范围

统计学方法只是辅助，最终判断需要结合专业知识。例如，土壤 pH 值在 3.5-9.5 范围之外几乎不可能：

# 专业知识判断
expert_ranges <- list(
  ph = c(3.5, 9.5),
  soc = c(0, 80),
  tn = c(0, 10),
  ca = c(0, 60)
)

# 检测超出专业知识范围的异常值
for (var_name in names(expert_ranges)) {
  range_vec <- expert_ranges[[var_name]]
  outliers <- soil |>
    filter(.data[[var_name]] < range_vec[1] | .data[[var_name]] > range_vec[2]) |>
    select(site_id, all_of(var_name))
  if (nrow(outliers) > 0) {
    cat(var_name, "超出范围 [", range_vec[1], "-", range_vec[2], "] 的样点：\n")
    print(outliers)
  }
}

生态学案例

在某植物多样性调查中，Z-score 方法识别出某样点物种数为 999，而正常范围在 10-80 之间。进一步调查发现，这是录入员手指滑键导致的输入错误。统计学异常检测帮助团队在数据使用前发现了这个错误，避免了错误的分析结果。

扩展记录：2026-04-09 | 目标字数：800+

18.4.2 IQR 方法

IQR（四分位距）方法是识别数值型变量异常值最常用的统计方法。它的核心思想是：用数据的中间 50% 来定义”正常范围”，超出这个范围的值被视为潜在异常值。

IQR 的计算公式

Q1 = 第25百分位数（下四分位数）
Q3 = 第75百分位数（上四分位数）
IQR = Q3 - Q1

异常值判定：
下界 = Q1 - 1.5 × IQR
上界 = Q3 + 1.5 × IQR

超出 [下界, 上界] 范围的值即为异常值。

在 R 中的实现

library(tidyverse)

# 定义 IQR 异常值检测函数
detect_outliers_iqr <- function(x, k = 1.5) {
  q1 <- quantile(x, probs = 0.25, na.rm = TRUE)
  q3 <- quantile(x, probs = 0.75, na.rm = TRUE)
  iqr <- q3 - q1
  lower <- q1 - k * iqr
  upper <- q3 + k * iqr
  return(x < lower | x > upper)
}

# 对 SOC 应用 IQR 方法
soc_outliers <- detect_outliers_iqr(soil$soc)
cat("检测到", sum(soc_outliers, na.rm = TRUE), "个 SOC 异常值\n")

# 查看异常值详情
soil |>
  filter(detect_outliers_iqr(soc)) |>
  select(site_id, soc, land_use)

k 值的调整

k = 1.5 是默认值，适用于大多数情况。但对于生态学数据，可以根据实际情况调整：

• k = 1.5：标准保守检测，检测约 0.7% 的数据为异常值
• k = 3：极端值检测，检测约 0.002% 的数据
• k = 1.0：更敏感的检测，适合噪声较大的野外数据

# 不同 k 值的比较
soil |>
  mutate(
    outlier_1.5 = detect_outliers_iqr(soc, k = 1.5),
    outlier_1.0 = detect_outliers_iqr(soc, k = 1.0),
    outlier_3.0 = detect_outliers_iqr(soc, k = 3.0)
  ) |>
  summarise(
    outliers_1.5 = sum(outlier_1.5, na.rm = TRUE),
    outliers_1.0 = sum(outlier_1.0, na.rm = TRUE),
    outliers_3.0 = sum(outlier_3.0, na.rm = TRUE)
  )

生态学案例

在某植物多样性调查中，研究者使用 IQR 方法（k=1.5）检测物种丰富度异常值。结果发现，同一个样点上记录的物种数为 999——这显然是录入错误（多打了一个9）。通过 IQR 方法，这个异常值被成功识别并修正为 99。这个案例说明：IQR 方法能有效检测录入错误导致的异常值。

注意事项

IQR 方法假设数据大致对称分布。对于严重偏态的数据（如很多零值的计数数据），可能需要先进行变换（如 log(x+1)）再应用 IQR 方法，或者使用基于分位数的其他方法。

扩展记录：2026-04-09 | 目标字数：800+

18.4.3 处理异常值

识别出异常值后，需要决定如何处理。处理方式取决于异常值的成因：如果是数据录入错误，应该修正；如果是真实的极端值，应该保留并记录；如果是测量误差，应该标记或删除。

常见的异常值处理策略

1. 修正：如果是明显的录入错误（如 999、-999），直接修正为正确值
2. 标记为缺失：保留原始值，转换为 NA，并在备注中说明原因
3. 截断：将极端值替换为合理边界值（不删除，但限制范围）
4. 删除行：如果异常值无法解释，且样本量充足，可以删除整行
5. 保留：如果是真实的生态现象（如喀斯特地区的高钙含量），应保留并讨论

在 R 中的实现

library(tidyverse)

# 策略1：修正明显的录入错误
soil_fixed <- soil |>
  mutate(dbh = if_else(dbh > 500, dbh / 10, dbh))

# 策略2：标记为缺失值（保留原始记录）
soil_flagged <- soil |>
  mutate(
    soc_flag = if_else(soc < 0 | soc > 100, "outlier", NA_character_),
    soc_clean = if_else(soc < 0 | soc > 100, NA_real_, soc)
  )

# 策略3：截断到合理范围
soil_trimmed <- soil |>
  mutate(soc_trimmed = pmax(0, pmin(80, soc)))

# 策略4：使用分组统计量处理
soil_grouped <- soil |>
  group_by(land_use) |>
  mutate(
    soc_z = (soc - mean(soc, na.rm = TRUE)) / sd(soc, na.rm = TRUE),
    soc_clean = if_else(abs(soc_z) > 3, NA_real_, soc)
  ) |>
  ungroup()

处理后的验证

# 查看处理结果
soil_comprehensive |>
  filter(!is.na(clean_flag)) |>
  select(site_id, land_use, soc, soc_clean, clean_flag)

# 对比处理前后的统计量
cat("\n处理前统计：\n")
print(summary(soil$soc))

cat("\n处理后统计：\n")
print(summary(soil_comprehensive$soc_clean))

生态学案例

在某土壤碳汇研究中，研究者发现部分样点的 SOC 值高达 150 g/kg，远超正常范围。通过追溯原始记录发现，这些样点的数据是将”15.0”误录为”150”。研究者将这类异常值标记为 NA，在报告中注明”3个样点因录入错误导致数据缺失”。这种透明的处理方式获得了审稿人的认可。

扩展记录：2026-04-09 | 目标字数：800+

18.4.4 喀斯特土壤特殊背景

根据 Li et al. (2017) 的研究: - 喀斯特地区土壤SOC一般在5-40 g/kg之间 - pH值通常在5.0-7.5之间(中性偏酸) - 可交换钙含量与岩性密切相关 - 负值或极端高值通常是测量或录入错误

# 根据专业知识设置合理范围,将超出范围的值替换为NA
soil_clean <- soil |>
  mutate(
    ph = if_else(ph < 3.5 | ph > 9.5, NA_real_, ph),        # pH合理范围
    soc = if_else(soc < 0 | soc > 80, NA_real_, soc),       # SOC合理范围(g/kg)
    tn = if_else(tn < 0 | tn > 10, NA_real_, tn),           # TN合理范围(g/kg)
    ca = if_else(ca < 0 | ca > 60, NA_real_, ca)             # Ca合理范围(cmol/kg)
  )

cat("异常值处理后的统计:\n")
summary(soil_clean |> select(soc, tn, ph, ca))

Note异常值 ≠ 错误值

这一节需要特别强调一个常见误区：异常值不等于错误值。在生态学研究中，异常值可能反映真实的生态现象，处理不当反而会掩盖重要发现。

三种不同性质的”异常值”

1. 数据录入错误：这是真正的错误，需要修正或删除
   • 示例：土壤有机碳 150 g/kg（应为 15.0 g/kg）
   • 特征：明显超出物理可能范围
2. 测量误差：设备故障或操作失误导致，需要评估影响
   • 示例：某次测量中仪器未校准，导致读数系统性偏高
   • 特征：可能呈现某种模式（如集中在某批次）
3. 真实的生态变异：这是宝贵的数据，反映了真实的生态现象
   • 示例：喀斯特石灰岩地区土壤钙含量极高（可达 60 cmol/kg）
   • 特征：有真实的生态学解释

如何区分？

library(tidyverse)

# 步骤1：检查是否超出物理可能范围
physical_impossible <- soil |>
  filter(soc < 0 | soc > 100)  # 物理上不可能

# 步骤2：检查是否超出合理生态学范围但可能有合理解释
unusual_ecological <- soil |>
  filter(soc > 80) |>
  mutate(possible_reason = case_when(
    lithology == "石灰岩" ~ "石灰岩地区SOC可能较高",
    land_use == "次生林" ~ "次生林SOC积累较多",
    TRUE ~ "需进一步调查"
  ))

cat("超出物理可能范围（疑似录入错误）：\n")
print(physical_impossible)
cat("\n超出合理范围但可能有生态学解释：\n")
print(unusual_ecological)

处理决策树

发现异常值后，按照以下流程判断：是否超出物理可能范围？如果是，则是录入错误；如果不是，则考虑是否有生态学合理解释；如果仍无法判断，进行敏感性分析。

敏感性分析

当你不确定某个值是否应该保留时，可以进行敏感性分析：分别在使用和不使用该值的情况下运行分析，比较结果是否一致。

# 敏感性分析示例
# 分析1：包含可疑值
result_with <- soil |>
  group_by(land_use) |>
  summarise(mean_soc = mean(soc, na.rm = TRUE))

# 分析2：排除可疑值
result_without <- soil |>
  filter(soc < 100) |>
  group_by(land_use) |>
  summarise(mean_soc = mean(soc, na.rm = TRUE))

# 比较结果
comparison <- result_with |>
  left_join(result_without, by = "land_use", suffix = c("_with", "_without")) |>
  mutate(diff = mean_soc_with - mean_soc_without)

cat("敏感性分析：\n")
print(comparison)
cat("\n如果差异小于5%，说明结果对异常值不敏感\n")

生态学案例

在某区域土壤碳储量研究中，统计方法检测到一个极端高值（85 g/kg）。研究者没有简单删除，而是查阅了该样点的原始记录：原来这是一个喀斯特石山区，土层薄但有机质含量高，这个值是真实的。保留该数据后，研究的结论更具代表性。审稿人对此给予高度评价。

扩展记录：2026-04-09 | 目标字数：800+

18.5 处理缺失值

18.5.1 识别缺失值

在处理异常值之后，下一步是系统性地识别数据中的缺失值。仅仅知道有多少缺失值是不够的，还需要了解缺失值的分布模式和成因。

识别缺失值的基本方法

library(tidyverse)

# 方法1：逐列统计缺失值数量
soil |>
  summarise(across(everything(), ~sum(is.na(.)))) |>
  pivot_longer(everything(), names_to = "variable", values_to = "n_missing") |>
  filter(n_missing > 0) |>
  arrange(desc(n_missing))

# 方法2：计算缺失比例
soil |>
  summarise(across(everything(), ~mean(is.na(.)) * 100)) |>
  pivot_longer(everything(), names_to = "variable", values_to = "pct_missing") |>
  filter(pct_missing > 0)

# 方法3：使用 naniar 包
library(naniar)
miss_var_summary(soil)

为什么会缺失？

在检查缺失值时，需要思考：为什么这些值会缺失？

• 随机缺失（MAR）：与其他变量相关，但本身与值无关
• 完全随机缺失（MCAR）：与任何变量都无关，纯属偶然
• 非随机缺失（MNAR）：与缺失值本身相关

# 检查缺失是否与某些变量相关
soil |>
  mutate(
    soc_missing = is.na(soc),
    ph_group = ntile(ph, 3)
  ) |>
  group_by(ph_group) |>
  summarise(
    n_total = n(),
    n_soc_missing = sum(soc_missing),
    pct_missing = mean(soc_missing) * 100
  )

生态学案例

在某长期生态监测研究中，研究者发现 SOC 的缺失值集中在2019年的数据中。进一步调查发现，那一年使用的土壤碳分析仪出现了故障，部分样品无法测定。这就是典型的MCAR——缺失与任何生态学变量无关。通过检查缺失模式，研究者避免了错误的统计推断。

扩展记录：2026-04-09 | 目标字数：800+

18.5.2 可视化缺失模式

除了数值统计，可视化是理解缺失模式的重要工具。通过图形可以直观地发现缺失值的空间分布、时间趋势和变量间关联。

常用的缺失值可视化方法

library(tidyverse)
library(naniar)

# 1. 缺失值热力图
soil |>
  mutate(row_id = row_number()) |>
  pivot_longer(-c(site_id, land_use, lithology, depth, date_raw),
               names_to = "variable", values_to = "value") |>
  mutate(is_missing = is.na(value)) |>
  ggplot(aes(x = variable, y = row_id, fill = is_missing)) +
  geom_tile(color = "white", linewidth = 0.2) +
  scale_fill_manual(values = c("lightblue", "red"), labels = c("有值", "缺失")) +
  labs(title = "缺失值分布热力图", x = "变量", y = "样点编号", fill = "") +
  theme_minimal()

# 2. UpSet 图（适合多变量缺失分析）
gg_miss_upset(soil |> select(soc, tn, ph, ca))

# 3. 缺失与变量的关系
gg_miss_var(soil, facet = land_use) +
  labs(title = "不同土地利用类型的缺失情况")

生态学案例

在某森林动态监测研究中，研究者通过缺失值热力图发现：pH 和 SOC 的缺失值高度重叠——几乎所有缺失 pH 的样点同时也缺失 SOC。进一步调查发现，这两种指标来自同一次实验室分析。通过可视化快速定位问题，研究者采取了针对性的重新分析方法，避免了不必要的数据删除。

扩展记录：2026-04-09 | 目标字数：800+

18.5.3 缺失值模式分析

仅仅知道有多少缺失值是不够的，还需要了解缺失值的模式（pattern）——即哪些变量倾向于同时缺失、缺失是否与某些已观测变量有关。理解缺失模式对于选择正确的填充方法至关重要。

三种经典的缺失机制：

1. 完全随机缺失（MCAR, Missing Completely at Random）

   缺失与任何变量都无关，纯属偶然。例如，野外采样时某仪器突然故障，导致当天所有样本的pH值都缺失。在这种情况下，删除缺失记录不会引入偏倚，但会损失样本量。

2. 随机缺失（MAR, Missing at Random）

   缺失与已观测变量相关，但与未观测变量无关。例如，在本数据集中，石灰岩地区的样本可能更容易缺失可交换钙数据（因为测量方法不同）。这种情况下，可以利用已观测变量（岩性）来帮助填充。

3. 非随机缺失（MNAR, Missing not at Random）

   缺失与未观测变量本身相关。例如，土壤有机碳含量极高的样本可能在实验室分析时被遗漏（因为超出了标准曲线的范围）。这种情况下，简单的填充方法会引入系统性偏倚。

使用 naniar 包分析缺失模式：

# 使用 naniar 包分析缺失模式
library(naniar)

# 查看缺失值的汇总统计
miss_var_summary(soil_clean)

# 查看每个变量的缺失比例
miss_var_pct(soil_clean)

# 可视化缺失与某个变量的关系（例如与land_use的关系）
gg_miss_var(soil_clean, facet = land_use)

生态学案例：

在一项长期森林样地研究中，研究者发现有些样本的叶片氮含量数据缺失。进一步分析发现，这些缺失的样本主要集中在2019年的冬季采样期。查阅采样记录发现，那年冬天使用的叶氮分析仪出现了故障，导致部分样品无法测定。这就是典型的MCAR——缺失与任何生物学变量无关，只与采样技术问题相关。

但如果缺失的样本主要集中在高温季节，且都是落叶树种，则可能是MNAR——因为夏季落叶树种没有叶片可供采样。这种情况下，简单地用均值填充会低估氮含量的季节变异。

upset图解读：

naniar包的gg_miss_upset()函数可以展示多变量缺失的组合模式。每个垂直的条形图表示同时缺失的变量组合，连接的点表示这些变量同时缺失。例如，如果SOC和TN的点被连接起来，说明这两个变量倾向于同时缺失——这可能是因为它们来自同一次实验室分析。

扩展记录：2026-04-09 | 目标字数：800+

18.5.4 处理策略

根据缺失的原因和比例,选择不同的处理方式:

策略	适用场景	R 实现
删除行	缺失比例低(<5%),随机缺失	drop_na()
删除列	某列缺失比例过高(>50%)	select(-col)
均值/中位数填充	数值型,缺失比例适中	replace_na()
分组填充	缺失值与分组有关	group_by() + fill()
标记为单独类别	分类变量	replace_na(list(x = "未知"))

在本数据集中: - SOC、TN、pH、Ca的缺失都是随机缺失(MAR) - 缺失比例都在可接受范围内(<10%) - 由于异常值已处理,均值和中位数差异不大,使用中位数填充

# 方法1:删除含缺失值的行
soil_complete <- soil_clean |> drop_na()
cat("删除缺失后:", nrow(soil_complete), "行\n")

# 方法2:用中位数填充数值型缺失值
soil_filled <- soil_clean |>
  mutate(
    across(
      .cols = c(soc, tn, ph, ca),
      .fns = ~replace_na(.x, median(.x, na.rm = TRUE))
    )
  )

cat("填充后缺失值数量:", sum(is.na(soil_filled)), "\n")

Warning缺失值处理的注意事项

缺失值处理是数据清洗中最需要谨慎的环节之一。不当的处理方式可能引入系统性偏倚，导致错误的结论。

核心原则：理解缺失机制再处理

在选择缺失值处理方法之前，必须先理解数据为什么会缺失。不同的缺失机制需要不同的处理策略：

MCAR（完全随机缺失） - 定义：缺失与任何变量都无关，纯属偶然 - 例子：某天实验室突然停电，导致当天所有样品数据缺失 - 适合方法：列表删除（listwise deletion）对参数估计无偏，但会损失样本量

MAR（随机缺失） - 定义：缺失与已观测变量相关，但与未观测变量无关 - 例子：低收入调查对象更不愿意报告收入，但收入缺失与教育程度相关 - 适合方法：多重插补（Multiple Imputation）

MNAR（非随机缺失） - 定义：缺失与未观测变量本身相关 - 例子：重病患者更可能缺失体检数据 - 适合方法：选择模型（Selection Model）或模式混合模型（Pattern Mixture Model）

在 R 中实现多重插补

library(mice)
library(tidyverse)

# 检查缺失模式
md.pattern(soil)

# 可视化缺失关系
library(VIM)
aggr(soil, col = c("navyblue", "red"), numbers = TRUE)

# 多重插补（生成5个完整数据集）
imp_data <- mice(soil, m = 5, method = "pmm", seed = 2027)

# 查看插补结果
print(imp_data)

# 对每个插补数据集进行分析
fit <- with(imp_data, lm(soc ~ land_use + ph))

# 汇总结果（Pool 合并）
pooled <- pool(fit)
summary(pooled)

绝对不能做的事

• 用 0 填充缺失的浓度数据（0 通常不是合理的自然值）
• 用均值填充后不报告填充比例
• 删除所有含缺失值的行后不说明样本量变化

报告要求

在论文方法部分，必须明确报告：各变量的缺失比例、选择的处理方法及理由、如使用插补需报告插补次数和汇总方法、以及敏感性分析结果（如有）。

# 报告缺失情况的代码模板
missing_report <- soil |>
  summarise(across(everything(), ~sum(is.na(.)))) |>
  pivot_longer(everything(), names_to = "variable", values_to = "n_missing") |>
  mutate(
    pct_missing = round(n_missing / nrow(soil) * 100, 1),
    n_valid = nrow(soil) - n_missing
  ) |>
  filter(n_missing > 0) |>
  arrange(desc(pct_missing))

cat("缺失值报告：\n")
print(missing_report)

生态学案例

某 meta 分析研究在处理缺失值时，研究团队对每种常用的缺失值处理方法都进行了敏感性分析。结果发现，使用均值填充和多重插补的结论方向一致，但效应量大小有显著差异。最终论文选择了更保守的多重插补结果，并在讨论部分专门说明了这一选择的影响。这个透明的处理方式得到了审稿人的认可。

扩展记录：2026-04-09 | 目标字数：800+

18.6 统一数据格式

18.6.1 日期格式统一

在生态学数据收集中，日期格式不一致是常见的数据质量问题。不同采样员可能使用不同的日期格式，如”2024-03-15”、“2024/03/15”、“15-Mar-2024”等。日期格式不统一会严重影响后续的时间序列分析。

常见的日期格式问题

格式类型	示例	问题
年-月-日（ISO）	2024-03-15	标准格式
年/月/日	2024/03/15	与 ISO 类似
日-月-年	15-Mar-2024	需要特殊解析
月/日/年	03/15/2024	美国常用格式
中文格式	2024年3月15日	需要特殊处理

在 R 中统一日期格式

library(tidyverse)
library(lubridate)

# 原始日期数据（包含多种格式）
date_raw <- c("2024-03-15", "2024/03/16", "15-Mar-2024", "2024-03-17")

# parse_date_time 处理多种格式
date_parsed <- parse_date_time(date_raw, orders = c("ymd", "dby", "mdy"))
cat("解析结果：\n")
print(date_parsed)

# 提取日期成分
soil <- soil |>
  mutate(
    date_parsed = parse_date_time(date_raw, orders = c("ymd", "dby", "mdy")),
    year = year(date_parsed),
    month = month(date_parsed),
    day = day(date_parsed),
    doy = yday(date_parsed)
  )

生态学案例

某研究团队在整理10年野外观测数据时，发现不同年份的日期格式完全不同：2015年用”2015/03/15”，2016-2018年用”15-Mar-2015”，2019年用”Mar-15-2015”。他们用 lubridate 的 parse_date_time 函数批量解析，成功将所有日期统一为标准格式，最终用于时间序列分析，揭示了土壤碳储量10年间的变化趋势。

扩展记录：2026-04-09 | 目标字数：800+

18.6.2 分类变量标准化

分类变量（如土地利用类型、物种名称、采样区等）在不同来源的数据中可能使用不同的编码方式。例如，“次生林”可能被记录为”次生林”、“次生林分”、“SSF”、“Secondary forest”等。统一分类变量的编码是数据清洗的重要步骤。

常见的分类变量问题

问题类型	示例	影响
同义不同名	“次生林” vs “次生林分”	被视为不同类别
拼写差异	“草地” vs “草 地”	空格导致匹配失败
大小写	“Forest” vs “forest”	被视为不同类别
缩写混用	“灌丛” vs “灌木丛”	需要统一标准

在 R 中标准化分类变量

library(tidyverse)

# 定义标准分类编码表
land_use_std <- tribble(
  ~raw, ~standard,
  "次生林", "次生林",
  "次生林分", "次生林",
  "SSF", "次生林",
  "Secondary forest", "次生林",
  "天然林", "天然林",
  "天然林分", "天然林"
)

# 应用标准化
soil <- soil |>
  left_join(land_use_std, by = c("land_use" = "raw")) |>
  mutate(land_use_std = if_else(is.na(standard), land_use, standard))

# 检查是否还有未匹配的类别
soil |>
  group_by(land_use_std) |>
  count()

生态学案例

某 meta 分析研究在合并不同数据集时，发现各国学者对”森林”类型的分类完全不同：巴西学者用”Floresta Secundaria”，中国学者用”次生林”，欧洲学者用”Secondary forest”。研究者建立了一个包含120个同义词的对照表，成功将所有数据统一到7个标准分类，极大地提高了 meta 分析的可靠性。

扩展记录：2026-04-09 | 目标字数：800+

18.7 数据类型转换

# 查看数据类型
str(soil_filled)

18.8 完整的数据清洗流水线

将以上步骤整合为一个完整的流水线:

# 定义清洗函数
clean_soil_data <- function(df) {
  df |>
    # 1. 去重
    distinct() |>
    # 2. 处理不合理的值(异常值 → NA)
    mutate(
      ph = if_else(ph < 3.5 | ph > 9.5, NA_real_, ph),
      soc = if_else(soc < 0 | soc > 80, NA_real_, soc),
      tn = if_else(tn < 0 | tn > 10, NA_real_, tn),
      ca = if_else(ca < 0 | ca > 60, NA_real_, ca)
    ) |>
    # 3. 统一日期格式
    mutate(
      date = parse_date_time(date_raw, orders = c("ymd", "y/m/d", "d-b-Y")) |> as.Date()
    ) |>
    # 4. 填充缺失值(中位数)
    mutate(
      across(
        .cols = c(soc, tn, ph, ca),
        .fns = ~replace_na(.x, median(.x, na.rm = TRUE))
      )
    ) |>
    # 5. 转换分类变量
    mutate(
      land_use = factor(land_use, levels = c("耕地", "草地", "灌丛", "次生林")),
      lithology = factor(lithology, levels = c("白云岩", "石灰岩")),
      depth = factor(depth, levels = c("0-10 cm", "10-20 cm", "20-30 cm"), ordered = TRUE)
    ) |>
    # 6. 选择和重命名列
    select(
      site_id, land_use, lithology, depth,
      soc, tn, ph, ca, date
    )
}

# 应用清洗流水线
soil_final <- clean_soil_data(soil_raw)

# 验证
cat("清洗后数据维度:", dim(soil_final), "\n")
summary(soil_final)
str(soil_final)

18.9 清洗前后对比报告

# 构建对比报告
report <- tibble(
  指标 = c("总行数", "重复行数", "SOC缺失", "TN缺失", "pH缺失", "Ca缺失",
           "SOC异常", "TN异常", "pH异常", "Ca异常", "日期格式数"),
  清洗前 = c(
    nrow(soil_raw),
    sum(duplicated(soil_raw)),
    sum(is.na(soil_raw$soc)),
    sum(is.na(soil_raw$tn)),
    sum(is.na(soil_raw$ph)),
    sum(is.na(soil_raw$ca)),
    sum(soil_raw$soc < 0 | soil_raw$soc > 80, na.rm = TRUE),
    sum(soil_raw$tn < 0 | soil_raw$tn > 10, na.rm = TRUE),
    sum(soil_raw$ph < 3.5 | soil_raw$ph > 9.5, na.rm = TRUE),
    sum(soil_raw$ca < 0 | soil_raw$ca > 60, na.rm = TRUE),
    length(unique(soil_raw$date_raw))
  ),
  清洗后 = c(
    nrow(soil_final),
    sum(duplicated(soil_final)),
    sum(is.na(soil_final$soc)),
    sum(is.na(soil_final$tn)),
    sum(is.na(soil_final$ph)),
    sum(is.na(soil_final$ca)),
    0, 0, 0, 0,
    1
  )
)

report

18.10 清洗后数据验证

18.10.1 按土地利用类型汇总

清洗后的数据，通常需要按分组变量进行汇总统计。土地利用类型是最常用的分组变量之一，因为它直接影响土壤理化性质。

按分组进行描述性统计

library(tidyverse)

# 按土地利用类型汇总数值变量
summary_by_group <- soil_clean |>
  group_by(land_use) |>
  summarise(
    across(
      .cols = c(soc, tn, ph, ca),
      .fns = list(
        n = ~sum(!is.na(.)),
        mean = ~mean(., na.rm = TRUE),
        sd = ~sd(., na.rm = TRUE)
      ),
      .names = "{col}_{fn}"
    )
  ) |>
  ungroup()

cat("各土地利用类型的 SOC 均值比较：\n")
print(summary_by_group |> select(land_use, soc_n, soc_mean, soc_sd))

组间差异检验

# 单因素方差分析（ANOVA）
aov_result <- aov(soc ~ land_use, data = soil_clean)
summary(aov_result)

# Kruskal-Wallis 检验（非参数替代）
kruskal.test(soc ~ land_use, data = soil_clean)

# 事后检验（Tukey HSD）
tukey_result <- TukeyHSD(aov_result)
print(tukey_result)

# 可视化
library(ggpubr)
ggboxplot(soil_clean, x = "land_use", y = "soc",
          color = "land_use", palette = "jco",
          add = "jitter") +
  stat_compare_means(method = "kruskal.test") +
  labs(title = "不同土地利用类型的土壤有机碳含量比较",
       y = "SOC (g/kg)", x = "土地利用类型")

生态学案例

在某区域土地利用变化研究中，研究者比较了耕地、草地、灌丛和次生林的土壤有机碳含量。结果显示：次生林的 SOC 均值（25.4 g/kg）显著高于耕地（8.5 g/kg），这与已有研究结论一致。通过组间差异检验，研究者为”保护现有森林以维持土壤碳汇功能”这一建议提供了数据支持。

扩展记录：2026-04-09 | 目标字数：800+

18.10.2 验证清洗结果

数据清洗不是做完就结束了，还需要验证清洗结果是否符合预期。这包括：数据分布是否合理、与已有知识是否一致、是否引入了新问题。

与外部知识的对比

将清洗后的数据与已发表研究中的数据进行比较，验证是否符合基本规律：

library(tidyverse)

# 与 Li et al. (2017) 的研究结果对比
cat("清洗后数据与文献对比：\n\n")

# SOC 含量规律：耕地 < 草地 < 灌丛 < 次生林
observed_order <- soil_clean |>
  group_by(land_use) |>
  summarise(mean_soc = mean(soc, na.rm = TRUE)) |>
  arrange(mean_soc)

cat("各土地利用类型 SOC 均值排序：\n")
print(observed_order)

# 检查是否符合预期规律
expected_order <- c("耕地", "草地", "灌丛", "次生林")
actual_order <- observed_order$land_use
cat("\n是否符合文献规律（耕地 < 草地 < 灌丛 < 次生林）：",
    ifelse(identical(actual_order, expected_order), "是", "否"), "\n")

# 岩性与 SOC 关系
cat("\n不同岩性地区 SOC 均值：\n")
lithology_comparison <- soil_clean |>
  group_by(lithology) |>
  summarise(mean_soc = mean(soc, na.rm = TRUE))
print(lithology_comparison)

数据分布检验

清洗后的数据分布应该合理，不应该出现人为制造的”伪规律”：

library(tidyverse)

# 检查各变量的分布
soil_clean |>
  select(where(is.numeric)) |>
  pivot_longer(everything(), names_to = "variable", values_to = "value") |>
  ggplot(aes(x = value)) +
  geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7) +
  facet_wrap(~variable, scales = "free") +
  labs(title = "清洗后数值变量分布检验") +
  theme_minimal()

# 检验正态性
for (var_name in c("soc", "tn", "ph")) {
  p_val <- shapiro.test(soil_clean[[var_name]])$p.value
  cat(var_name, "正态性检验 p-value:", round(p_val, 4))
  cat(ifelse(p_val < 0.05, " 非正态", " 正态"), "\n")
}

变量间相关性检验

清洗后的变量间关系应该符合生态学常识：

library(corrplot)

# 计算相关矩阵
numeric_vars <- soil_clean |>
  select(soc, tn, ph, ca) |>
  cor(use = "pairwise.complete.obs")

# 可视化相关矩阵
corrplot(numeric_vars, method = "ellipse", type = "upper")

cat("SOC 与 TN 相关系数：", round(cor(soil_clean$soc, soil_clean$tn, use = "pairwise"), 3), "\n")
cat("SOC 与 Ca 相关系数：", round(cor(soil_clean$soc, soil_clean$ca, use = "pairwise"), 3), "\n")
cat("预期：SOC 与 TN、Ca 应呈正相关（有机质积累伴随氮和钙增加）\n")

完整性检验

清洗后的数据应该没有影响分析完整性的问题：

# 样本量检验
cat("清洗前后样本量：\n")
cat("原始数据：", nrow(soil), "行\n")
cat("清洗后数据：", nrow(soil_clean), "行\n")
cat("损失比例：", round((1 - nrow(soil_clean)/nrow(soil)) * 100, 1), "%\n")

# 各组样本量是否足够
g_by_group <- soil_clean |>
  group_by(land_use) |>
  count()
cat("\n各土地利用类型样本量：\n")
print(g_by_group)

# 最低样本量检验（每组至少5个）
min_n <- min(g_by_group$n)
cat("\n最小样本量：", min_n)
if (min_n >= 5) {
  cat("（满足基本统计检验要求）\n")
} else {
  cat("（样本量不足，建议重新采样或合并组别）\n")
}

生态学案例

某研究团队在完成数据清洗后，照例与已发表数据对比。结果发现清洗后 SOC 含量与预期相反（耕地 > 次生林）。进一步排查发现，清洗过程中误删了次生林的大部分数据，导致该组均值偏高。修正后，结论恢复正常。这个案例说明：验证步骤不是形式，而是发现错误的最后防线。

扩展记录：2026-04-09 | 目标字数：800+

18.11 课后练习

1. 下载 Li et al. (2017) 的原始数据(如有),或从以下来源获取真实数据:
   • 中国科学院广西环江喀斯特农田生态系统国家野外科学观测研究站
   • Dryad 数据平台(https://datadryad.org)
2. 检查数据中的缺失值、异常值和格式问题
3. 编写一个完整的数据清洗流水线
4. 对比清洗前后的数据质量(缺失值数量、异常值数量、数据类型)
5. 将清洗脚本和清洗报告提交到 GitHub

18.12 参考文献

Li, D., Wen, L., Yang, L., Luo, P., Xiao, W., Chen, H., … & Wang, K. (2017). Dynamics of soil organic carbon and nitrogen following agricultural abandonment in a karst region. Journal of Geophysical Research: Biogeosciences, 122(11), 2849-2864. DOI: 10.1002/2016JG003683