6  R 语言基础

7 R 语言基础

本章介绍 R 语言的基本数据类型、数据读取和基础统计操作。如果你之前没有编程经验，不用担心——R 的语法直观易懂，我们会从最基础的概念开始。

7.1 R 语言简史与生态学应用

7.1.1 R 的诞生与发展

R 语言由新西兰奥克兰大学的 Ross Ihaka 和 Robert Gentleman 于 1993 年开发，名字取自两位创始人名字的首字母。R 是 S 语言的开源实现，而 S 语言由贝尔实验室开发，专门用于统计计算和图形展示。

发展里程碑：

• 1993：R 项目启动
• 1995：R 成为自由软件（GPL 协议）
• 1997：R Core Team 成立，负责维护和开发
• 2000：R 1.0.0 正式发布
• 2004：R 2.0.0 发布，引入懒惰加载（lazy loading）
• 2013：R 3.0.0 发布，支持长向量（> 2^31 元素）
• 2020：R 4.0.0 发布，改进字符串处理
• 2024：R 4.4.0 发布，持续优化性能

7.1.2 为什么生态学家选择 R？

1. 专为统计设计：R 的核心就是统计分析，内置大量统计函数
2. 丰富的生态学包：vegan（群落生态学）、ade4（多元分析）、biodiversityR（生物多样性）等
3. 强大的可视化：ggplot2 可以制作出版级图表
4. 可重复性研究：代码即分析流程，便于审查和重现
5. 活跃的社区：遇到问题可以在 Stack Overflow、RStudio Community 求助
6. 免费开源：无需购买许可证，适合学术研究

7.1.3 生态学中的 R 应用实例

# 群落生态学：物种多样性分析
library(vegan)
data(dune)  # 荷兰沙丘草地植被数据
diversity(dune, index = "shannon")  # Shannon 多样性指数

# 空间生态学：物种分布模型
library(dismo)
# 使用 MaxEnt 模型预测物种分布

# 系统发育分析：构建进化树
library(ape)
# 读取序列数据，构建系统发育树

# 景观生态学：景观格局分析
library(landscapemetrics)
# 计算景观指数（斑块数量、形状指数等）

NoteR vs Python：生态学家该选哪个？

R 的优势： - 统计分析更成熟（线性混合模型、广义可加模型等） - 生态学专用包更多（vegan、phyloseq、picante） - 数据可视化更优雅（ggplot2）

Python 的优势： - 机器学习更强大（scikit-learn、TensorFlow） - 通用编程能力更强（网络爬虫、自动化） - 地理空间分析更灵活（geopandas、rasterio）

建议：两者都学！R 用于统计分析和可视化，Python 用于数据采集和机器学习。

7.2 认识 RStudio

RStudio 界面分为四个区域：

• 左上：脚本编辑器（Source），编写和保存代码
• 左下：控制台（Console），直接运行代码
• 右上：环境（Environment），查看当前变量
• 右下：文件/图表/帮助（Files/Plots/Help）

Tip快捷键

• Ctrl + Enter：运行当前行代码
• Ctrl + Shift + Enter：运行整个脚本
• Ctrl + S：保存文件
• Tab：代码自动补全

7.3 基本数据类型

R 中有几种基本的数据类型：

7.3.1 数值型（numeric）

数值型是 R 中最基础、最常用的数据类型之一，主要用于保存可以参与数学运算的数值。它既可以表示整数，也可以表示带小数的实数。在 R 中，如果直接写 3.14、6.5、100 这类数字，通常都会被识别为数值型；如果想明确告诉 R 某个值是整数，可以在数字后面加上 L，例如 42L。学习数值型的意义，不只是会“存一个数字”，更重要的是理解：生态学研究中的绝大多数连续变量——例如土壤 pH、含水量、有机碳含量、植物高度、胸径、降雨量、温度——本质上都要以数值型形式进入分析流程。只有变量被正确识别为数值型，R 才能顺利执行平均值、标准差、回归分析、绘图等操作。如果数据读入后看起来像数字，但实际上被识别成字符型，那么很多统计函数会直接报错。因此，判断一个对象是不是数值型、在需要时进行类型转换，是 R 初学者必须掌握的基本功。

# 创建数值型对象
x <- 3.14        # 小数，默认是 numeric
y <- 42L         # 整数，使用 L 后缀
z <- c(5.2, 6.1, 7.3, 5.8)  # 数值向量

# 查看对象类型
class(x)
class(y)
class(z)

# 对数值型执行常见运算
mean(z)          # 计算平均值
sd(z)            # 计算标准差
round(z, 1)      # 保留 1 位小数

运行结果说明：class(x) 一般会返回 "numeric"，class(y) 通常会返回 "integer"，说明 R 会区分一般数值和显式整数。z 是由多个数值组成的向量，因此仍然属于数值型数据结构。mean(z) 会返回这组数值的平均值，sd(z) 会给出离散程度，round(z, 1) 则把结果保留到 1 位小数。这里最关键的认识是：只要对象是数值型，R 就能直接把统计函数逐步应用上去，这也是后续进行描述统计、假设检验和建模的前提。

从土壤生态学或马尾松研究的角度看，数值型变量几乎贯穿整个分析流程。例如，在马尾松人工林土壤研究中，我们常记录不同样地的土壤 pH、土壤含水率、全氮、有效磷、SOC（soil organic carbon，土壤有机碳）和微生物生物量碳等指标。这些指标都不是“标签”，而是可以比较大小、计算均值和建立相关关系的连续量。如果某次调查中 20 个样地的 pH 被错误读成字符型，那么你既无法计算平均 pH，也无法判断哪几个样地偏酸，更无法进一步分析 pH 与马尾松生长、凋落物分解或土壤酶活性的关系。相反，只要这些变量被正确保存为数值型，就可以马上做分布统计、箱线图、相关分析，甚至进一步构建线性模型，解释马尾松林地土壤环境差异。因此，数值型不仅是“一个数据类型”，更是生态学定量分析能够成立的技术基础。

7.3.2 字符型（character）

字符型用于存储文本信息，也就是不能直接参与数学运算、但对数据解释非常重要的内容。在 R 中，凡是放在单引号或双引号中的内容，通常都会被识别为字符型，例如物种名称、样地编号、采样者姓名、处理组名称、土层标签等。字符型并不表示“没用的数据”，恰恰相反，它在生态学研究中承担着数据标识、分类描述和结果解释的作用。比如“马尾松”“杉木”“对照组”“处理组”“0-10 cm”“2024-春季调查”这些值都应该用字符型来保存。初学者常犯的错误是把字符型当成数值型使用，或者误把本应是分类变量的内容写成数字代码却不做说明，结果后续分析时难以理解。掌握字符型的关键，不只是会写字符串，而是知道什么时候应该保留原始文本信息、什么时候需要把字符型进一步转换为因子或其他类型。

# 创建字符型对象
name <- "广西大学"
species <- '马尾松'     # 单引号和双引号都可以
plot_id <- c("P01", "P02", "P03")

# 查看类型和字符长度
class(name)
class(species)
nchar(name)            # 统计字符数

# 字符串拼接与输出
message_text <- paste("优势树种为", species)
message_text

运行结果说明：class(name) 和 class(species) 都会返回 "character"，说明这些对象被 R 当作文本处理。nchar(name) 会返回字符串包含的字符个数，常用于检查样地编码、物种名缩写或文本长度是否符合要求。paste() 会把不同对象拼接成一句完整文本，运行后 message_text 会得到类似“优势树种为 马尾松”的结果。这里要注意，字符型对象不能直接用于加减乘除运算，否则会报错；但它们非常适合承担“标签”和“说明”的功能。

在土壤生态学和马尾松研究中，字符型数据经常出现在元数据和分类信息中。例如，我们在调查马尾松林不同坡位的土壤性质时，可能会记录 species = "马尾松"、slope_position = "上坡位"、soil_layer = "0-10 cm"、treatment = "施氮" 等字段。这些内容本身不是连续数值，但它们决定了每一行数据的生态学含义。如果没有这些字符型信息，后续分析就会失去背景：你虽然能看到 pH 为 5.8、SOC 为 18.4，但不知道它来自哪一个样地、哪一层土壤、哪一种处理。尤其在马尾松林地长期定位监测中，样地编号、调查季节、群落类型和样品来源往往都需要用字符型精确记录，然后再与数值型指标联合分析。因此，字符型并不是“附属信息”，而是生态数据整理和结果复现的重要组成部分。一个规范的数据表，通常总是由字符型字段负责描述背景，由数值型字段负责承载测量值，两者缺一不可。

7.3.3 逻辑型（logical）

逻辑型（logical）是 R 中用于表示真假判断的数据类型，只有两个可能的取值：TRUE 和 FALSE（可简写为 T 和 F，但不推荐简写以避免与变量名混淆）。逻辑型数据在生态学数据分析中扮演着”筛选器”和”判断器”的角色，是数据清洗、条件筛选、逻辑判断的核心工具。

is_tree <- TRUE
is_herb <- FALSE
3 > 2              # 比较运算返回逻辑值

逻辑型的生成方式：

1. 直接赋值：is_native <- TRUE 表示某物种为本地种。

2. 比较运算：所有比较运算符（>、<、>=、<=、==、!=）的结果都是逻辑型。例如，soil_ph > 7 会对每个 pH 值判断是否大于 7，返回一个逻辑向量。

3. 逻辑运算：&（与）、|（或）、!（非）可以组合多个逻辑条件。例如，(elevation > 1000) & (slope < 30) 筛选出海拔高于 1000 米且坡度小于 30 度的样地。

4. 函数返回：许多 R 函数返回逻辑值，如 is.na(x) 判断是否为缺失值，is.numeric(x) 判断是否为数值型。

生态学应用场景：

• 数据筛选：data[data$dbh > 10, ] 筛选出胸径大于 10 cm 的树木记录。逻辑向量 data$dbh > 10 作为行索引，只保留 TRUE 对应的行。

• 条件分类：ifelse(soil_ph < 5.5, "酸性", "中性或碱性") 根据 pH 值将土壤分类。ifelse() 函数根据逻辑条件返回不同的值。

• 缺失值检测：is.na(tree_height) 返回一个逻辑向量，标记哪些树高数据缺失。sum(is.na(tree_height)) 可以统计缺失值数量（TRUE 在数值运算中被视为 1，FALSE 被视为 0）。

• 数据验证：all(biomass >= 0) 检查所有生物量值是否非负，any(is.na(species)) 检查是否存在物种名称缺失。

注意事项：

1. 逻辑运算符的区别：& 和 | 是向量化运算符，逐元素比较；&& 和 || 只比较第一个元素，常用于 if 语句中的单一条件判断。

2. 缺失值的逻辑运算：任何涉及 NA 的比较运算结果都是 NA，而非 FALSE。例如，NA > 5 返回 NA。因此，筛选数据时需要显式排除缺失值：data[!is.na(data$ph) & data$ph > 7, ]。

3. 逻辑型的数值转换：TRUE 转换为数值时为 1，FALSE 为 0。这一特性常用于统计：sum(soil_ph > 7) 统计碱性土壤样本数量，mean(is_native) 计算本地种比例。

实战示例：假设我们有一批马尾松林样地数据，需要筛选出”海拔 800-1200 米、坡度小于 25 度、土壤 pH 在 4.5-6.0 之间”的样地进行进一步分析。可以用逻辑运算组合多个条件：

selected <- (elevation >= 800) & (elevation <= 1200) &
            (slope < 25) &
            (soil_ph >= 4.5) & (soil_ph <= 6.0)
filtered_data <- data[selected, ]

这种基于逻辑型的数据筛选方式，比手动逐行检查数据表高效得多，且筛选逻辑完全可追溯、可重复。掌握逻辑型数据的使用，是从”手工整理数据”过渡到”编程式数据管理”的关键一步。

7.4 向量（vector）

向量是 R 中最基本的数据结构，用 c() 函数创建：

# 数值向量
heights <- c(12.5, 8.3, 15.7, 10.2, 9.8)

# 字符向量
species <- c("马尾松", "杉木", "桉树", "樟树", "楠木")

# 逻辑向量
is_native <- c(TRUE, TRUE, FALSE, TRUE, TRUE)

7.4.1 向量的深入理解

向量是 R 中最基础的数据结构，理解向量是掌握 R 的关键。

向量的特性：

1. 原子性（Atomic）：向量中所有元素必须是同一类型
2. 一维性：向量只有长度，没有维度
3. 向量化运算：运算自动应用到每个元素

# 类型强制转换（coercion）
mixed <- c(1, 2, "three")  # 数值被强制转换为字符
class(mixed)
mixed

# 类型优先级：logical < integer < numeric < character
c(TRUE, 1, 2.5, "text")  # 全部转为字符型

# 命名向量（named vector）
soil_ph <- c(plot1 = 6.5, plot2 = 5.8, plot3 = 7.2)
soil_ph
names(soil_ph)  # 提取名称
soil_ph["plot1"]  # 用名称索引

7.4.2 向量生成函数

向量生成函数是 R 初学者必须掌握的一组基础工具，因为在实际分析中，我们经常需要快速创建一串有规律的数据，而不是一个一个手动输入。所谓“向量生成函数”，本质上就是帮助我们批量构造向量的函数，包括生成顺序序列、按固定步长取值、重复元素以及构造随机数样本等。最常见的有 :、seq()、rep()、runif()、rnorm() 和 rpois()。这些函数之所以重要，是因为生态学数据分析并不总是从现成文件开始：很多时候你需要先构建测试数据、模拟样本、设定时间序列，或者创建样地编号、深度序列、处理水平序列等辅助变量。掌握它们后，你会发现很多本来看起来繁琐的准备工作，其实几行代码就能完成。更进一步说，向量生成函数不仅节省时间，也能减少人工录入错误，让数据结构更规范、可重复。

# 序列生成
1:10                       # 生成 1 到 10 的整数序列
seq(1, 10, by = 2)         # 指定步长为 2
seq(0, 1, length.out = 5)  # 指定输出长度为 5

# 重复生成
rep(1:3, times = 3)        # 整个向量重复 3 次
rep(1:3, each = 3)         # 每个元素分别重复 3 次
rep(c("A", "B"), times = c(2, 3))  # 不同元素重复不同次数

# 随机数生成
set.seed(123)              # 固定随机种子，保证结果可重复
runif(5, min = 0, max = 10)   # 生成 5 个 0~10 的均匀分布随机数
rnorm(5, mean = 10, sd = 2)   # 生成 5 个正态分布随机数
rpois(5, lambda = 3)          # 生成 5 个泊松分布计数值

运行结果说明：1:10 会直接返回从 1 到 10 的整数序列；seq(1, 10, by = 2) 会返回 1、3、5、7、9；seq(0, 1, length.out = 5) 会把 0 到 1 之间平均切成 5 个点。rep() 系列主要用于重复元素，常适合构造实验处理分组或重复测量标签。随机数函数则更适合做模拟：runif() 适合生成均匀分布样本，rnorm() 常用于模拟近似正态分布的生态指标，rpois() 常用于模拟个体数、事件次数等计数型变量。set.seed() 很关键，因为它能确保你和别人再次运行代码时得到相同结果，从而保证教学演示和科研复现实验的一致性。

在土壤生态学和马尾松研究中，向量生成函数的应用非常广。例如，如果你要建立一个马尾松林样地数据框，首先就需要创建样地编号 plot_id <- paste0("P", 1:20)，再用 rep() 生成 4 个坡位类型各重复 5 次的分组变量；如果要分析不同土层的养分变化，可以用 seq(0, 40, by = 10) 快速生成采样深度标签；如果要在课堂上演示土壤 pH 或有机碳含量的统计分析，又暂时没有真实数据，就可以用 rnorm() 模拟一组接近实测分布的数据；如果要模拟土壤动物个体数或根际微生物菌落数，则常会用到 rpois()。在马尾松人工林土壤养分研究中，这些函数不仅能帮助我们构造教学示例，还能用来测试代码、预演分析流程、检查模型是否能正确运行。换句话说，向量生成函数是把“分析思路”快速变成“可运行数据对象”的桥梁。

7.4.3 生态学实例：土壤 pH 值分析

土壤 pH 值分析是土壤生态学中极其常见的入门任务，因为 pH 会直接影响养分有效性、微生物活性、根系吸收以及凋落物分解速率。在马尾松林地中，土壤通常偏酸性，不同坡位、不同林龄、不同经营措施下，pH 往往存在明显差异。因此，学习如何用 R 对 pH 向量进行生成、汇总、分类和异常值识别，是把基础语法转化为生态分析能力的一个典型例子。这个案例的重点不是“记住几个函数”，而是理解：当一组生态测量值进入 R 后，我们如何快速把它变成可解释的统计信息。

# 模拟 30 个样地的土壤 pH 值
set.seed(2024)
soil_ph <- rnorm(30, mean = 6.2, sd = 0.8)

# 基本统计
cat("平均 pH:", round(mean(soil_ph), 2), "\n")
cat("标准差:", round(sd(soil_ph), 2), "\n")
cat("变异系数:", round(sd(soil_ph) / mean(soil_ph) * 100, 2), "%\n")

# 分类统计
acidic <- sum(soil_ph < 6.5)                    # 酸性土壤数量
neutral <- sum(soil_ph >= 6.5 & soil_ph <= 7.5) # 中性土壤数量
alkaline <- sum(soil_ph > 7.5)                  # 碱性土壤数量

cat("\n土壤类型分布：\n")
cat("酸性:", acidic, "个样地\n")
cat("中性:", neutral, "个样地\n")
cat("碱性:", alkaline, "个样地\n")

# 异常值检测：超出均值 ± 2 倍标准差
outliers <- soil_ph[abs(soil_ph - mean(soil_ph)) > 2 * sd(soil_ph)]
if (length(outliers) > 0) {
  cat("\n检测到异常值:", round(outliers, 2), "\n")
} else {
  cat("\n未检测到异常值\n")
}

运行结果说明：代码首先利用 rnorm() 生成 30 个近似正态分布的土壤 pH 值，然后分别计算平均值、标准差和变异系数。平均值反映总体酸碱水平，标准差表示样地间波动大小，变异系数便于比较不同指标的离散程度。后面的分类统计把连续 pH 值转换为“酸性—中性—碱性”三类，更适合生态解释。异常值检测则帮助研究者识别可能存在测量误差、样品污染或特殊环境条件的样地。如果输出结果显示大多数样地集中在酸性范围，这与许多南方马尾松林土壤偏酸的背景是一致的。

从生态学解释看，这类分析可以直接服务于马尾松人工林土壤质量评价。假设你比较 30 个马尾松样地的 pH，发现大部分样地低于 6.5，说明土壤总体偏酸，这可能与长期针叶凋落物积累、有机酸释放和强降雨淋溶有关。如果其中少数样地 pH 明显偏高，研究者就要进一步检查这些样地是否位于不同母质区、是否经过石灰处理，或者是否受到人为扰动。在土壤养分研究中，pH 还常与速效磷、交换性铝、微生物群落结构一起分析，因为酸化往往会影响养分可利用性和微生物过程。对于马尾松研究而言，pH 不只是一个“数字”，它还是理解林地土壤环境、解释林分生长差异和制定经营措施的重要基础指标。通过这个案例，学生应当学会把一个普通的数值向量，转化为具有生态含义的分析结果。

7.4.4 向量运算

R 的向量运算是自动逐元素进行的，这一特性称为”向量化运算”（vectorization），是 R 语言区别于传统编程语言的核心优势之一。向量化运算不仅代码简洁，而且执行效率高，避免了显式循环带来的性能损失。在生态学数据分析中，向量化运算使得对整列数据的批量处理变得极为便捷。

heights <- c(12.5, 8.3, 15.7, 10.2, 9.8)

# 基本运算
heights * 100          # 每个元素乘以 100（cm 转换）
heights + 1            # 每个元素加 1

# 统计函数
mean(heights)          # 平均值
sd(heights)            # 标准差
max(heights)           # 最大值
min(heights)           # 最小值
length(heights)        # 长度

向量化运算的原理：

当对向量执行算术运算时，R 会自动将运算应用到向量的每一个元素上。例如，heights * 100 等价于 c(12.5*100, 8.3*100, 15.7*100, 10.2*100, 9.8*100)，但前者的代码更简洁，执行速度也更快（因为 R 底层用 C 语言实现了向量化运算）。

向量间的运算：

两个等长向量可以直接进行算术运算，运算逐元素配对进行：

dbh <- c(15.2, 12.8, 18.5, 14.3, 16.7)  # 胸径 (cm)
height <- c(12.5, 8.3, 15.7, 10.2, 9.8)  # 树高 (m)

# 计算树木体积（简化公式：V = 0.0001 * dbh^2 * height）
volume <- 0.0001 * dbh^2 * height

如果两个向量长度不等，R 会自动”循环补齐”短向量（recycling），但这可能导致意外结果，应当谨慎使用。

常用统计函数：

R 提供了丰富的向量统计函数，这些函数是生态学数据探索性分析的基础工具：

函数	功能	生态学应用示例
mean(x)	平均值	样地平均物种丰富度
median(x)	中位数	树高中位数（对异常值不敏感）
sd(x)	标准差	土壤有机碳含量的变异程度
var(x)	方差	群落结构的空间异质性
min(x), max(x)	最小值、最大值	环境因子的取值范围
range(x)	范围（最小值和最大值）	海拔梯度范围
sum(x)	求和	样方内总生物量
length(x)	向量长度	样本数量
quantile(x)	分位数	计算 25%、50%、75% 分位数

向量化运算的生态学应用：

1. 单位转换：野外调查数据常需要单位转换。例如，树高从米转换为厘米（height_cm <- height_m * 100），土壤有机碳从 g/kg 转换为百分比（soc_percent <- soc_g_kg / 10）。

2. 标准化处理：将数据标准化为均值为 0、标准差为 1 的 z-score，便于比较不同量纲的变量：z_score <- (x - mean(x)) / sd(x)。

3. 对数转换：生物量、物种多度等数据常呈右偏分布，对数转换可改善正态性：log_biomass <- log(biomass + 1)（加 1 避免对数零值）。

4. 比率计算：计算碳氮比（cn_ratio <- carbon / nitrogen）、根冠比（root_shoot_ratio <- root_biomass / shoot_biomass）等生态学指标。

注意事项：

1. 缺失值处理：如果向量中包含 NA，大多数统计函数会返回 NA。需要添加 na.rm = TRUE 参数排除缺失值：mean(heights, na.rm = TRUE)。

2. 除零错误：向量运算中若出现除以零的情况，R 会返回 Inf（无穷大）或 NaN（非数值）。应在运算前检查分母是否为零。

3. 向量长度不匹配：两个向量运算时，若长度不等且短向量长度不是长向量长度的整数倍，R 会发出警告。应确保参与运算的向量长度一致或有意使用循环补齐。

实战示例：假设我们测量了 50 株马尾松的胸径和树高，需要计算每株树的材积并统计总材积。使用向量化运算可以一行代码完成：

dbh <- c(...)  # 50 个胸径值
height <- c(...)  # 50 个树高值
volume <- 0.0001 * dbh^2 * height  # 逐株计算材积
total_volume <- sum(volume)  # 总材积
mean_volume <- mean(volume)  # 平均材积

这种向量化的编程风格，不仅代码简洁易读，而且执行效率远高于逐个元素循环计算。掌握向量化运算，是从”面向过程编程”过渡到”面向数据编程”的关键一步。

7.4.5 向量索引

heights <- c(12.5, 8.3, 15.7, 10.2, 9.8)

heights[1]             # 第 1 个元素（R 从 1 开始计数！）
heights[c(1, 3)]       # 第 1 和第 3 个
heights[heights > 10]  # 筛选大于 10 的元素

7.5 因子（factor）

因子用于表示分类变量，在生态学数据中非常常见：

# 创建因子
habitat <- factor(c("森林", "草地", "湿地", "森林", "草地", "森林"))
habitat

# 查看水平
levels(habitat)

# 频次统计
table(habitat)

Note为什么需要因子？

在统计分析和绘图中，R 需要知道哪些变量是分类变量。因子的核心作用是告诉R：“这个变量有固定的几个类别”，而不是随意的文本。例如，在生态学研究中，“土地利用类型”可能包含”林地、草地、农田”三个类别；“采样地点”可能包含”A、B、C”三个地点。如果把这类变量存储为普通文本（character），R会把它们当作没有任何顺序或关系的字符串处理。但当你需要进行方差分析、绘制箱线图、或者建立回归模型时，R需要知道这些变量的分类结构——这时候就需要因子。

因子的实际用途：

因子在三种常见场景中必不可少。第一，分类统计——如果你想知道不同土地利用类型间的土壤有机碳差异，table() 或 group_by() 会自动按因子水平分组统计。第二，绘图——ggplot2的fill和color参数要求分类变量必须是因子，否则无法正确映射图例和颜色。第三，回归分析——在建立线性模型时，lm() 会将因子变量自动编码为虚拟变量（dummy variables），这是统计分析的基础。

library(tidyverse)

# 创建因子变量
site <- c("林地", "草地", "林地", "农田", "草地")
# 普通文本向量（不是因子）
site_char <- site
class(site_char)  # [1] "character"

# 转为因子
site_factor <- factor(site, levels = c("林地", "草地", "农田"))
class(site_factor)  # [1] "factor"
levels(site_factor)  # [1] "林地" "草地" "农田"

# 因子在数据分析中的实际应用
soil_data <- tibble(
  site = factor(c("林地", "草地", "林地", "农田", "草地", "林地"),
                 levels = c("林地", "草地", "农田")),
  oc = c(28.5, 15.3, 31.2, 8.7, 12.1, 26.8),  # 有机碳 g/kg
  ph = c(5.2, 5.8, 5.1, 6.2, 5.9, 5.3)
)

# 按因子水平分组统计
soil_data |>
  group_by(site) |>
  summarise(mean_oc = mean(oc), n = n())

# 因子绘图：ggplot2 自动识别因子
ggplot(soil_data, aes(x = site, y = oc, fill = site)) +
  geom_boxplot() +
  labs(title = "不同土地利用类型的土壤有机碳差异",
       x = "土地利用类型", y = "有机碳 (g/kg)") +
  theme_minimal()

# 在回归模型中，因子自动转为虚拟变量
model <- lm(oc ~ site, data = soil_data)
summary(model)
# R 自动将"林地"设为参照组，输出"草地"和"农田"相对于林地的回归系数

生态学案例：

在分析马尾松混交林土壤孔隙数据时，一位学生把所有变量都存成了character类型。结果在建立回归模型时，R把”林地”和”草地”当成了数值型变量，导致模型报错。改正方法是使用factor()将土地利用类型转为因子，并在模型中正确指定参照组。这个错误让论文数据重新分析浪费了一周时间。在数据导入阶段就明确因子结构，比事后修正要高效得多。

扩展记录: 2026-04-10 | 扩展者：Clawd | 目标字数：800+

7.5.1 有序因子（ordered factor）

有些分类变量有自然的顺序关系，比如”低 < 中 < 高”。使用有序因子可以让 R 正确处理这种顺序：

# 创建有序因子
soil_depth <- factor(
  c("0-10cm", "10-20cm", "20-30cm", "0-10cm", "10-20cm"),
  levels = c("0-10cm", "10-20cm", "20-30cm"),
  ordered = TRUE
)
soil_depth

# 有序因子支持比较运算
soil_depth[1] < soil_depth[3]  # TRUE：0-10cm < 20-30cm

7.6 数据框（data.frame）

数据框是 R 中最常用的数据结构，类似于 Excel 表格：

# 创建数据框
tree_data <- data.frame(
  species = c("马尾松", "杉木", "桉树", "樟树", "楠木"),
  height = c(12.5, 8.3, 15.7, 10.2, 9.8),
  dbh = c(25.3, 15.8, 32.1, 20.5, 18.7),    # 胸径 (cm)
  native = c(TRUE, TRUE, FALSE, TRUE, TRUE)
)

# 查看数据
tree_data
str(tree_data)     # 查看结构
summary(tree_data) # 描述性统计

7.6.1 数据框的深入理解

数据框是列表的特殊形式，每列可以是不同类型，但同一列内必须是相同类型。

# 数据框的本质
class(tree_data)  # "data.frame"
typeof(tree_data)  # "list"
is.list(tree_data)  # TRUE

# 数据框的属性
names(tree_data)    # 列名
rownames(tree_data) # 行名
dim(tree_data)      # 维度（行数，列数）
nrow(tree_data)     # 行数
ncol(tree_data)     # 列数

# 数据框与矩阵的区别
# 数据框：列可以是不同类型
# 矩阵：所有元素必须是同一类型

7.6.2 数据框的创建方式

数据框（data.frame）是 R 中最重要的数据结构，也是生态学数据分析的核心载体。理解数据框的多种创建方式，可以帮助你灵活地组织和准备数据。数据框本质上是一个特殊的列表，每一列可以是不同的数据类型（数值、字符、因子等），但同一列内的元素必须是相同类型，且所有列的长度必须相等。这种结构完美对应了生态学数据的典型形式：每一行代表一个观测单位（如一个样地、一棵树、一次测量），每一列代表一个变量（如物种、高度、pH 值）。掌握数据框的创建方法，不仅能让你从零开始构建数据集（如模拟数据、手动录入小样本数据），还能帮助你理解数据框的内部结构，为后续的数据操作和转换打下基础。在实际研究中，虽然大部分数据来自文件读取，但在编写函数、生成报告、构建示例时，手动创建数据框的能力仍然非常重要。

# 方式1：从向量创建（最基础的方法）
species <- c("马尾松", "杉木", "桉树")
height <- c(12.5, 8.3, 15.7)
df1 <- data.frame(species, height)
df1

# 方式2：直接在 data.frame() 中定义列
df2 <- data.frame(
  plot_id = 1:5,
  soil_ph = c(6.5, 5.8, 7.2, 6.1, 5.5),
  organic_matter = c(3.2, 2.8, 4.1, 3.5, 2.9)
)
df2

# 方式3：从矩阵转换
mat <- matrix(1:12, nrow = 4, ncol = 3)
df3 <- as.data.frame(mat)
colnames(df3) <- c("N", "P", "K")
rownames(df3) <- paste0("样地", 1:4)
df3

# 方式4：使用 tibble（tidyverse 推荐）
library(tibble)
df4 <- tibble(
  species = c("马尾松", "杉木", "桉树"),
  height = c(12.5, 8.3, 15.7)
)
df4  # tibble 打印更友好，自动显示数据类型

# 方式5：从列表创建
data_list <- list(
  plot_id = 1:3,
  species = c("马尾松", "杉木", "桉树"),
  dbh = c(25.3, 15.8, 32.1)
)
df5 <- as.data.frame(data_list)
df5

# 方式6：逐列添加（动态构建）
df6 <- data.frame(plot_id = 1:5)
df6$species <- c("马尾松", "杉木", "桉树", "樟树", "楠木")
df6$height <- c(12.5, 8.3, 15.7, 10.2, 9.8)
df6

# 方式7：使用 expand.grid 创建所有组合
df7 <- expand.grid(
  plot = 1:3,
  treatment = c("对照", "施氮", "施磷"),
  layer = c("0-10cm", "10-20cm")
)
df7  # 生成 3 × 3 × 2 = 18 行的完整组合

# 方式8：从其他数据框复制并修改
df8 <- df2
df8$CN_ratio <- df8$organic_matter / 0.5  # 添加新列
df8 <- df8[1:3, ]  # 只保留前 3 行
df8

# 方式9：使用 tribble（行式输入，适合小数据）
library(tibble)
df9 <- tribble(
  ~species,    ~height, ~dbh,
  "马尾松",    12.5,    25.3,
  "杉木",      8.3,     15.8,
  "桉树",      15.7,    32.1
)
df9

# 方式10：读取剪贴板（快速从 Excel 复制）
# df10 <- read.table("clipboard", header = TRUE, sep = "\t")

运行结果说明：data.frame() 会创建一个标准的数据框对象，默认会将字符向量转换为因子（R 4.0 之前），可以用 stringsAsFactors = FALSE 避免。tibble() 创建的对象打印时会显示数据类型（如 <chr>、<dbl>），且不会自动转换字符为因子，更符合现代数据分析习惯。expand.grid() 非常适合生成实验设计的所有处理组合。tribble() 的行式输入方式在手动录入小样本数据时非常直观。

生态学应用案例：在马尾松林土壤研究中，不同的创建方式适用于不同场景：

library(tibble)
library(dplyr)

# 场景1：手动录入样地基本信息
plot_info <- tibble(
  plot_id = paste0("P", 1:5),
  slope_position = c("上坡", "中坡", "下坡", "中坡", "上坡"),
  aspect = c("南", "东南", "东", "南", "西南"),
  elevation = c(450, 420, 380, 410, 460)
)

# 场景2：生成实验设计的所有处理组合
experimental_design <- expand.grid(
  plot = 1:4,
  treatment = c("对照", "低氮", "高氮"),
  sampling_time = c("春季", "夏季", "秋季", "冬季")
) %>%
  arrange(plot, treatment, sampling_time) %>%
  mutate(sample_id = paste(plot, treatment, sampling_time, sep = "_"))

head(experimental_design, 12)

# 场景3：从测定结果构建数据框
# 假设实验室返回了 3 个样地的测定结果
soil_results <- tibble(
  plot_id = c("P01", "P02", "P03"),
  pH = c(6.5, 5.8, 7.2),
  SOC = c(18.5, 15.2, 22.3),
  TN = c(1.2, 1.0, 1.5)
) %>%
  mutate(
    CN_ratio = SOC / TN,
    soil_type = case_when(
      pH < 6.0 ~ "酸性",
      pH < 7.0 ~ "弱酸性",
      TRUE ~ "中性"
    )
  )

print(soil_results)

# 场景4：合并多个来源的数据
# 野外调查数据
field_data <- tibble(
  plot_id = paste0("P", 1:3),
  canopy_cover = c(0.85, 0.78, 0.92),
  litter_depth = c(3.5, 2.8, 4.2)
)

# 实验室数据
lab_data <- tibble(
  plot_id = paste0("P", 1:3),
  pH = c(6.5, 5.8, 7.2),
  SOC = c(18.5, 15.2, 22.3)
)

# 合并
complete_data <- field_data %>%
  left_join(lab_data, by = "plot_id")

print(complete_data)

这些例子展示了数据框创建在实际研究中的多样化应用。掌握这些方法后，你可以灵活地组织数据、设计实验、整合多源信息，为后续的统计分析和可视化打下坚实基础。

7.6.3 访问数据框

数据框（data frame）是 R 中最常用的数据结构，它将多个向量组织成表格形式，每一列是一个变量，每一行是一个观测。访问数据框的数据有多种方式，掌握这些方法是数据操作的基础。

tree_data <- data.frame(
  species = c("马尾松", "杉木", "桉树", "樟树", "楠木"),
  height = c(12.5, 8.3, 15.7, 10.2, 9.8),
  dbh = c(25.3, 15.8, 32.1, 20.5, 18.7),
  native = c(TRUE, TRUE, FALSE, TRUE, TRUE)
)

tree_data$height       # 用 $ 访问列
tree_data[1, ]         # 第 1 行
tree_data[, 2]         # 第 2 列
tree_data[1, 2]        # 第 1 行第 2 列

数据框访问的三种主要方式：

1. 列访问（按变量名）：
   • tree_data$height：用 $ 符号提取单列，返回向量。这是最常用的方式，代码简洁且易读。
   • tree_data[["height"]]：用双方括号提取单列，返回向量。适合变量名存储在字符串中的情况（如 col_name <- "height"; tree_data[[col_name]]）。
   • tree_data["height"]：用单方括号提取单列，返回数据框（包含一列）。注意与双方括号的区别。
2. 行列索引（按位置）：
   • tree_data[1, ]：提取第 1 行，返回数据框（包含一行）。
   • tree_data[, 2]：提取第 2 列，返回向量。
   • tree_data[1, 2]：提取第 1 行第 2 列的单个元素。
   • tree_data[1:3, c(1, 3)]：提取第 1-3 行、第 1 和第 3 列，返回数据框。
3. 逻辑索引（按条件筛选）：
   • tree_data[tree_data$height > 10, ]：筛选树高大于 10 米的所有行。
   • tree_data[tree_data$native == TRUE, ]：筛选本地树种。
   • tree_data[tree_data$species %in% c("马尾松", "杉木"), ]：筛选特定物种。

生态学应用场景：

• 提取单个变量进行统计：mean(tree_data$height) 计算平均树高，sd(tree_data$dbh) 计算胸径标准差。

• 筛选特定条件的样本：large_trees <- tree_data[tree_data$dbh > 20, ] 筛选出胸径大于 20 cm 的树木，native_trees <- tree_data[tree_data$native == TRUE, ] 筛选本地树种。

• 提取多列进行分析：tree_data[, c("height", "dbh")] 提取树高和胸径两列，用于相关性分析或散点图绘制。

• 按行提取样本：tree_data[1:3, ] 提取前 3 行数据，用于查看数据结构或调试代码。

高级访问技巧：

1. 使用 subset() 函数：subset(tree_data, height > 10 & native == TRUE) 筛选树高大于 10 米的本地树种，语法更简洁。

2. 使用 dplyr 包的 filter() 和 select()：

   library(dplyr)
   tree_data %>%
     filter(height > 10) %>%
     select(species, height)

   这种管道式语法更符合数据分析的思维流程。

3. 使用 which() 函数获取行号：which(tree_data$height > 10) 返回满足条件的行号向量，可用于更复杂的索引操作。

注意事项：

1. 单方括号 vs 双方括号：tree_data["height"] 返回数据框，tree_data[["height"]] 返回向量。在需要向量的场合（如统计函数），应使用双方括号或 $。

2. 逻辑索引中的缺失值：如果逻辑条件中包含 NA，筛选结果会包含 NA 行。应使用 !is.na() 显式排除：tree_data[!is.na(tree_data$height) & tree_data$height > 10, ]。

3. 行列索引的顺序：tree_data[行索引, 列索引]，行在前，列在后。省略行索引表示所有行，省略列索引表示所有列。

4. 负索引：tree_data[-1, ] 排除第 1 行，tree_data[, -2] 排除第 2 列。负索引不能与正索引混用。

实战示例：假设我们有一批马尾松林样地调查数据，需要筛选出”树高大于 12 米、胸径大于 20 cm 的本地树种”进行进一步分析：

# 方法 1：基础 R 语法
selected_trees <- tree_data[tree_data$height > 12 &
                            tree_data$dbh > 20 &
                            tree_data$native == TRUE, ]

# 方法 2：dplyr 语法（推荐）
library(dplyr)
selected_trees <- tree_data %>%
  filter(height > 12, dbh > 20, native == TRUE)

# 进一步提取需要的列
selected_trees <- selected_trees %>%
  select(species, height, dbh)

掌握数据框的访问方法，是从”查看数据”过渡到”操作数据”的关键一步。在实际研究中，数据筛选、变量提取、条件查询是最常见的操作，熟练运用这些技巧可以大幅提高数据处理效率。

7.7 读取数据文件

实际研究中，数据通常存储在文件中。R 提供了多种读取方式。

Tip安装所需包

读取数据前，需要先安装相应的 R 包。只需安装一次，之后每次使用 library() 加载即可：

# 安装 tidyverse（包含 readr 等常用包）
install.packages("tidyverse")

# 安装 readxl（读取 Excel 文件）
install.packages("readxl")

7.8 R 包安装与管理

R 的强大之处在于丰富的扩展包（packages）。理解包管理是使用 R 的关键技能。

7.8.1 什么是 R 包？

R 包（package）是函数、数据集和文档的标准化集合，用于扩展 R 的基础功能。可以将 R 包理解为”工具箱”——基础 R 提供了基本工具，而各种扩展包则提供了专业领域的高级工具。

R 包的核心组成：

1. 函数库：实现特定功能的 R 函数（如 ggplot2 包的 ggplot() 函数用于绘图）
2. 数据集：示例数据或参考数据（如 datasets 包的 iris 数据集）
3. 文档：函数帮助文档、使用教程（vignettes）和引用信息
4. 依赖关系：包可能依赖其他包，安装时会自动处理依赖

包的来源与数量：

截至 2024 年，CRAN（Comprehensive R Archive Network，R 官方包仓库）上有超过 20,000 个包，涵盖统计分析、机器学习、生物信息学、空间分析、可视化等几乎所有数据科学领域。此外，还有 Bioconductor（生物信息学专用）、GitHub（开发版包）等其他来源。

生态学常用包示例：

包名	功能	典型应用
vegan	群落生态学分析	物种多样性指数、排序分析（NMDS、CCA）
ggplot2	数据可视化	绘制出版级图表
dplyr	数据处理	筛选、汇总、变换数据
sf	空间数据处理	地理信息系统（GIS）分析
lme4	混合效应模型	嵌套设计、重复测量数据分析

为什么需要安装包？

R 默认只加载少量核心包（如 base、stats、graphics），这样可以保持 R 启动速度快、内存占用小。当需要使用特定功能时，才安装和加载相应的包。这种设计类似于手机应用商店——系统自带基础功能，专业功能通过安装 App 获得。

包的版本管理：

R 包会持续更新以修复 bug 和添加新功能。使用 packageVersion("包名") 可以查看已安装包的版本。某些分析需要特定版本的包以确保可重复性，此时可以使用 remotes::install_version() 安装指定版本。

生态学案例：某研究团队在分析土壤微生物群落数据时，使用了 vegan 包进行 NMDS 排序分析，ggplot2 包绘制排序图，dplyr 包处理环境因子数据。这三个包协同工作，大幅简化了分析流程。如果没有这些包，研究者需要从零编写数百行代码实现相同功能。

注意事项：

• 包名区分大小写（ggplot2 不等于 Ggplot2）
• 安装包需要网络连接（从 CRAN 下载）
• 某些包依赖系统库（如 sf 包需要 GDAL），安装前需检查系统环境

扩展记录: 2026-04-11 | 扩展者：Clawd | 目标字数：800+

7.8.2 包的安装

R 包的安装是扩展 R 功能的第一步,也是每个 R 用户必须掌握的基本技能。所谓”安装包”,就是把别人编写好的函数、数据和文档下载到你的计算机上,让 R 能够识别和使用这些扩展功能。R 包的安装只需要执行一次,之后每次使用时只需要用 library() 加载即可,不需要重复安装。理解包安装的重要性在于:R 的基础功能虽然强大,但生态学研究中需要的专业分析工具(如群落排序、系统发育分析、空间统计等)都封装在各种专业包中。掌握包的安装方法,意味着你可以站在前人的肩膀上,直接使用经过验证的分析方法,而不需要从零开始编写复杂算法。R 包主要来自三个来源:CRAN(官方仓库,最常用)、Bioconductor(生物信息学专用)和 GitHub(开发版本)。不同来源的包有不同的安装方法,但 CRAN 包的安装最简单,也是初学者最常接触的。

# 方法1:从 CRAN 安装单个包(最常用)
install.packages("ggplot2")

# 方法2:安装多个包(一次性安装,节省时间)
install.packages(c("dplyr", "tidyr", "readr"))

# 方法3:安装 tidyverse(包含多个常用包的集合)
install.packages("tidyverse")

# 方法4:从 GitHub 安装开发版(需要先安装 devtools)
install.packages("devtools")
devtools::install_github("username/packagename")

# 方法5:从 Bioconductor 安装生物信息学包
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("phyloseq")

# 方法6:指定 CRAN 镜像(加速下载)
install.packages("vegan", repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")

# 方法7:安装本地包文件(.tar.gz 或 .zip)
install.packages("path/to/package.tar.gz", repos = NULL, type = "source")

# 检查包是否已安装
if (!require("ggplot2")) {
  install.packages("ggplot2")
}

运行结果说明:install.packages() 会从 CRAN 下载包并安装到你的 R 库目录(通常是 ~/R/library/)。安装过程中会显示下载进度和编译信息,成功后会提示”package ‘xxx’ successfully unpacked”。如果包有依赖项(即需要其他包才能运行),R 会自动安装这些依赖包。devtools::install_github() 用于安装 GitHub 上的开发版包,这些包可能包含最新功能但稳定性较差。BiocManager::install() 专门用于安装 Bioconductor 的生物信息学包,这些包在基因组学、微生物组学研究中非常重要。

生态学应用案例:在马尾松林土壤生态学研究中,你可能需要安装以下包来完成不同的分析任务:

# 群落生态学分析
install.packages("vegan")      # 多样性指数、排序分析
install.packages("biodiversityR")  # 生物多样性可视化

# 土壤微生物组分析
BiocManager::install("phyloseq")  # 微生物群落分析
BiocManager::install("DESeq2")    # 差异丰度分析

# 空间分析
install.packages("sf")         # 空间数据处理
install.packages("terra")      # 栅格数据分析

# 数据处理与可视化
install.packages("tidyverse")  # 数据科学工具集
install.packages("ggplot2")    # 高级绘图

# 统计建模
install.packages("lme4")       # 线性混合效应模型
install.packages("mgcv")       # 广义可加模型

这些包覆盖了从数据读取、清洗、分析到可视化的完整流程。安装完成后,你就可以用 library(vegan) 加载包并使用其中的函数。建议在项目开始时一次性安装所有需要的包,避免分析过程中频繁中断安装。

7.8.3 包的加载与使用

# 加载包（每次会话都需要）
library(ggplot2)

# 或者使用 require()（返回逻辑值）
if (!require(ggplot2)) {
  install.packages("ggplot2")
  library(ggplot2)
}

# 使用包中的函数而不加载整个包
dplyr::filter(data, condition)

# 查看已加载的包
search()

# 卸载包
detach("package:ggplot2", unload = TRUE)

7.8.4 包的管理

# 查看已安装的包
installed.packages()

# 更新所有包
update.packages()

# 更新特定包
update.packages("ggplot2")

# 删除包
remove.packages("packagename")

# 查看包的版本
packageVersion("ggplot2")

# 查看包的帮助文档
help(package = "ggplot2")

# 查看包中的函数列表
ls("package:ggplot2")

7.8.5 常用生态学 R 包

生态学研究涉及的数据类型和分析方法非常多样,从群落组成、物种多样性、系统发育关系,到空间分布、功能性状、生态系统过程,每个领域都有专门的 R 包提供支持。了解并掌握这些常用包,是开展现代生态学研究的必备技能。这些包不仅提供了经过同行评审的标准分析方法,还包含了丰富的示例数据和详细的文档,可以帮助研究者快速上手。更重要的是,使用这些广泛认可的包进行分析,可以提高研究的可重复性和可比性,便于与其他研究者交流和合作。在马尾松林土壤生态学研究中,我们通常需要整合群落生态学(植物和微生物群落)、土壤化学、空间分析等多个方面的工具,因此熟悉不同类型包的功能和适用场景非常重要。

# === 群落生态学 ===
install.packages("vegan")      # 群落分析核心包:多样性指数、排序分析(NMDS、PCA、CCA)
install.packages("ade4")       # 多元分析:对应分析、主成分分析
install.packages("biodiversityR")  # 生物多样性分析与可视化

# vegan 包使用示例
library(vegan)
data(dune)  # 荷兰沙丘草地植被数据
diversity(dune, index = "shannon")  # 计算 Shannon 多样性指数
nmds <- metaMDS(dune)  # 非度量多维尺度分析
plot(nmds)

# === 系统发育分析 ===
install.packages("ape")        # 系统发育树构建与操作
install.packages("phytools")   # 系统发育工具:祖先状态重建、性状进化
install.packages("picante")    # 系统发育群落生态学:系统发育多样性

# ape 包使用示例
library(ape)
tree <- read.tree(text = "((A:1,B:1):1,(C:1,D:1):1);")
plot(tree)

# === 空间生态学 ===
install.packages("sp")         # 空间数据类(旧版,但仍广泛使用)
install.packages("sf")         # 简单要素(Simple Features,推荐使用)
install.packages("raster")     # 栅格数据处理(旧版)
install.packages("terra")      # 栅格数据处理(新版,推荐使用)
install.packages("spdep")      # 空间依赖性分析:Moran's I、空间回归

# sf 包使用示例
library(sf)
# 创建空间点数据
points <- st_as_sf(data.frame(
  lon = c(108.3, 108.4, 108.5),
  lat = c(22.8, 22.9, 23.0),
  species = c("马尾松", "杉木", "桉树")
), coords = c("lon", "lat"), crs = 4326)

# === 功能性状分析 ===
install.packages("FD")         # 功能多样性指数
install.packages("TPD")        # 性状概率密度
install.packages("fundiversity")  # 功能多样性计算

# === 数据处理与可视化 ===
install.packages("tidyverse")  # 数据科学工具集(包含 dplyr、ggplot2、tidyr 等)
install.packages("ggplot2")    # 高级数据可视化
install.packages("dplyr")      # 数据处理:筛选、分组、汇总
install.packages("tidyr")      # 数据整理:宽长格式转换

# === 统计建模 ===
install.packages("lme4")       # 线性混合效应模型(LMM、GLMM)
install.packages("nlme")       # 非线性混合效应模型
install.packages("mgcv")       # 广义可加模型(GAM)
install.packages("glmmTMB")    # 广义线性混合模型(支持更多分布)

# === 微生物组分析(需要 Bioconductor) ===
BiocManager::install("phyloseq")  # 微生物群落分析核心包
BiocManager::install("DESeq2")    # 差异丰度分析
BiocManager::install("microbiome")  # 微生物组数据处理

运行结果说明:这些包安装后,可以通过 library() 加载使用。例如 library(vegan) 会加载 vegan 包,之后就可以使用 diversity()、metaMDS() 等函数。每个包都有详细的帮助文档,可以通过 ?diversity 或 help(package = "vegan") 查看。

生态学应用案例:在马尾松人工林土壤研究中,典型的分析流程可能需要以下包的组合:

# 1. 数据读取与清洗
library(readxl)  # 读取 Excel 数据
library(dplyr)   # 数据筛选、分组

# 2. 土壤理化性质分析
library(ggplot2)  # 绘制土壤 pH、有机碳分布图
library(corrplot)  # 绘制土壤指标相关矩阵

# 3. 植物群落分析
library(vegan)  # 计算物种多样性、群落排序

# 4. 土壤微生物群落分析
library(phyloseq)  # 处理 16S rRNA 测序数据
library(DESeq2)    # 比较不同处理间微生物差异

# 5. 空间分析
library(sf)     # 处理样地空间坐标
library(spdep)  # 分析土壤性质的空间自相关

# 6. 统计建模
library(lme4)   # 构建混合效应模型,考虑样地随机效应

这个工作流展示了如何整合多个包完成复杂的生态学分析。建议在项目开始时创建一个”包加载脚本”,统一管理所有需要的包,确保分析的可重复性。

7.8.6 包冲突处理

不同包可能有同名函数，导致冲突：

# 查看函数来自哪个包
conflicts()

# 明确指定使用哪个包的函数
dplyr::filter(data, condition)  # 使用 dplyr 的 filter
stats::filter(data, condition)  # 使用 stats 的 filter

# 使用 conflicted 包管理冲突
install.packages("conflicted")
library(conflicted)
conflict_prefer("filter", "dplyr")  # 优先使用 dplyr 的 filter

7.8.7 生态学实例：安装并使用 vegan 包

# 安装 vegan 包
install.packages("vegan")

# 加载包
library(vegan)

# 使用内置数据集
data(dune)  # 荷兰沙丘草地植被数据
head(dune)

# 计算 Shannon 多样性指数
shannon <- diversity(dune, index = "shannon")
shannon

# 进行 NMDS 排序
nmds <- metaMDS(dune)
plot(nmds)

7.8.8 读取 CSV 文件

CSV (Comma-Separated Values, 逗号分隔值) 文件是生态学研究中最常用的数据交换格式之一。CSV 文件本质上是纯文本文件,每行代表一条记录,字段之间用逗号分隔,第一行通常是列名。CSV 格式的优势在于:跨平台兼容性好(Windows、Mac、Linux 都支持)、文件体积小、可以用任何文本编辑器打开查看、几乎所有统计软件都能读取。在生态学数据管理中,野外调查数据、实验室测定结果、文献整理的数据集通常都保存为 CSV 格式。R 提供了多种读取 CSV 文件的方法,基础 R 的 read.csv() 功能完善但速度较慢,tidyverse 的 read_csv() 速度更快且能自动识别数据类型,是目前推荐的方法。掌握 CSV 文件读取不仅是数据分析的起点,也是实现数据共享和分析可重复性的基础。

# 方法1:基础 R 的 read.csv()(传统方法)
data <- read.csv("data/species_survey.csv")
# 常用参数:
# header = TRUE: 第一行是列名(默认)
# sep = ",": 分隔符(默认逗号)
# stringsAsFactors = FALSE: 不自动转换为因子(R 4.0+ 默认)
data <- read.csv("data/species_survey.csv",
                 header = TRUE,
                 stringsAsFactors = FALSE)

# 方法2:tidyverse 的 read_csv()(推荐)
library(readr)
data <- read_csv("data/species_survey.csv")
# 优势:
# - 速度快(比 read.csv 快 10 倍)
# - 自动识别数据类型
# - 返回 tibble(改进的数据框)
# - 显示读取进度和列类型

# 方法3:指定列类型(避免自动识别错误)
data <- read_csv("data/species_survey.csv",
                 col_types = cols(
                   plot_id = col_character(),
                   species = col_character(),
                   dbh = col_double(),
                   height = col_double(),
                   date = col_date(format = "%Y-%m-%d")
                 ))

# 方法4:跳过前几行(处理带注释的文件)
data <- read_csv("data/species_survey.csv", skip = 2)

# 方法5:只读取部分行(预览大文件)
data_preview <- read_csv("data/species_survey.csv", n_max = 10)

# 方法6:处理中文编码问题
data <- read_csv("data/species_survey.csv", locale = locale(encoding = "UTF-8"))
# 或者
data <- read_csv("data/species_survey.csv", locale = locale(encoding = "GBK"))

# 方法7:读取 URL 上的 CSV 文件
url <- "https://example.com/data.csv"
data <- read_csv(url)

# 读取后检查数据
head(data)      # 查看前 6 行
str(data)       # 查看数据结构
summary(data)   # 描述性统计

运行结果说明:read.csv() 会返回一个 data.frame 对象,read_csv() 返回 tibble 对象(tibble 是改进的 data.frame,打印时更友好)。读取成功后,R 会显示数据的维度(行数和列数)。read_csv() 还会显示每列的数据类型(如 col_character()、col_double())。如果文件路径错误,会报错”cannot open file”。如果编码不匹配,中文可能显示为乱码,需要调整 locale 参数。

生态学应用案例:在马尾松林土壤调查中,野外数据通常记录在 Excel 中,然后导出为 CSV 格式供 R 分析。假设你有一个 soil_survey.csv 文件,包含样地编号、坡位、土层、pH、有机碳等字段:

library(readr)

# 读取土壤调查数据
soil_data <- read_csv("data/soil_survey.csv",
                      col_types = cols(
                        plot_id = col_character(),
                        slope_position = col_factor(levels = c("上坡", "中坡", "下坡")),
                        soil_layer = col_character(),
                        pH = col_double(),
                        SOC = col_double(),
                        TN = col_double()
                      ))

# 检查数据
head(soil_data)
summary(soil_data)

# 数据清洗:移除缺失值
soil_data_clean <- soil_data %>%
  filter(!is.na(pH), !is.na(SOC))

# 按坡位分组统计
library(dplyr)
soil_summary <- soil_data_clean %>%
  group_by(slope_position) %>%
  summarise(
    n = n(),
    mean_pH = mean(pH),
    mean_SOC = mean(SOC),
    sd_SOC = sd(SOC)
  )

print(soil_summary)

这个流程展示了从 CSV 文件读取数据到初步分析的完整过程。通过指定列类型,可以确保 slope_position 被正确识别为有序因子,pH 和 SOC 被识别为数值型,避免后续分析出错。在实际研究中,建议将原始数据保存为 CSV 格式并妥善备份,所有数据清洗和分析步骤都用 R 脚本记录,确保分析的可重复性。

7.8.9 读取 Excel 文件

Excel 文件（.xlsx 或 .xls）是生态学研究中极其常见的数据格式，尤其在实验室数据记录、野外调查数据整理和跨学科协作中广泛使用。虽然 CSV 格式更适合程序化处理和长期存档，但 Excel 在数据录入、初步整理和与非编程人员协作时仍然不可替代。

Excel 文件的优势与局限：

优势： 1. 多工作表支持：一个文件可以包含多个工作表，便于组织不同类型的数据（如”样地信息”、“土壤数据”、“植被数据”分别存储在不同工作表中）。 2. 格式和公式：支持单元格格式、条件格式、公式计算，便于数据录入时的即时检查和计算。 3. 直观易用：图形界面友好，非编程人员也能轻松使用。 4. 广泛兼容：几乎所有数据分析软件都支持 Excel 格式。

局限： 1. 文件格式复杂：Excel 文件是二进制格式，不是纯文本，难以用版本控制系统（如 Git）追踪变化。 2. 跨平台兼容性：不同版本的 Excel 可能导致格式丢失或数据错误。 3. 数据类型自动转换：Excel 会自动将某些字符串转换为日期或数值（如基因名称 “SEPT2” 被转换为日期），导致数据损坏。 4. 文件体积大：相比 CSV，Excel 文件体积更大，不适合存储大规模数据。

推荐工作流程：用 Excel 录入和初步整理数据 → 在 R 中读取并进行正式分析 → 将清洗后的数据保存为 CSV 格式以便长期保存和共享。

# 方法1:使用 readxl 包(推荐,无需外部依赖)
library(readxl)
data <- read_excel("data/soil_data.xlsx", sheet = 1)
# 常用参数:
# sheet: 指定工作表(可以用数字或名称)
# range: 指定读取范围(如 "A1:D10")
# col_names: 是否第一行是列名(默认 TRUE)
# skip: 跳过前几行

# 方法2:读取指定工作表(用名称)
data <- read_excel("data/soil_data.xlsx", sheet = "土壤理化性质")

# 方法3:读取指定范围
data <- read_excel("data/soil_data.xlsx",
                   sheet = 1,
                   range = "A1:F50")

# 方法4:跳过前几行(处理带标题的表格)
data <- read_excel("data/soil_data.xlsx", skip = 3)

# 方法5:指定列类型
data <- read_excel("data/soil_data.xlsx",
                   col_types = c("text", "text", "numeric", "numeric", "date"))

# 方法6:查看 Excel 文件中的所有工作表
excel_sheets("data/soil_data.xlsx")

# 方法7:批量读取所有工作表
library(purrr)
all_sheets <- excel_sheets("data/soil_data.xlsx")
data_list <- map(all_sheets, ~read_excel("data/soil_data.xlsx", sheet = .x))
names(data_list) <- all_sheets

# 方法8:使用 openxlsx 包(可以写入 Excel)
library(openxlsx)
data <- read.xlsx("data/soil_data.xlsx", sheet = 1)

# 读取后检查数据
head(data)
str(data)
summary(data)

readxl 包的核心优势：

1. 无外部依赖：不需要安装 Java 或其他外部软件，纯 R 实现，安装和使用都非常简单。
2. 速度快：底层用 C++ 实现，读取速度远快于其他 Excel 读取包。
3. 稳定可靠：由 RStudio 团队维护，兼容性好，不易出错。
4. 自动类型推断：能够智能识别数值、日期、文本等数据类型，减少手动指定的工作量。

常见问题与解决方案：

1. 工作表名称包含中文或特殊字符：直接用中文名称即可，read_excel() 支持 UTF-8 编码。如果出现乱码，检查 Excel 文件的保存编码。

2. Excel 文件包含合并单元格：readxl 会将合并单元格的值填充到第一个单元格，其余单元格为 NA。建议在 Excel 中取消合并单元格后再读取。

3. Excel 文件包含公式：readxl 读取的是公式计算后的值，而非公式本身。如果需要公式，应在 Excel 中将公式转换为数值后再读取。

4. 日期格式问题：Excel 的日期存储为数值（距离 1900-01-01 的天数），readxl 会自动转换为 R 的日期格式。如果转换错误，可以手动指定 col_types = "date"。

5. 读取速度慢：如果 Excel 文件很大（> 10 MB），建议先在 Excel 中另存为 CSV 格式，然后用 read_csv() 读取，速度会快得多。

生态学应用案例：

假设我们有一个马尾松林土壤调查的 Excel 文件，包含三个工作表：“样地信息”、“土壤理化性质”、“微生物数据”。我们需要读取所有工作表并合并分析：

library(readxl)
library(dplyr)

# 查看所有工作表名称
sheets <- excel_sheets("data/pine_forest_soil.xlsx")
print(sheets)  # "样地信息" "土壤理化性质" "微生物数据"

# 读取各工作表
plot_info <- read_excel("data/pine_forest_soil.xlsx", sheet = "样地信息")
soil_chem <- read_excel("data/pine_forest_soil.xlsx", sheet = "土壤理化性质")
microbe <- read_excel("data/pine_forest_soil.xlsx", sheet = "微生物数据")

# 合并数据（按样地编号）
complete_data <- plot_info %>%
  left_join(soil_chem, by = "plot_id") %>%
  left_join(microbe, by = "plot_id")

# 保存为 CSV 格式（便于后续分析和共享）
write_csv(complete_data, "data/pine_forest_soil_clean.csv")

最佳实践建议：

1. Excel 用于数据录入，CSV 用于数据分析：Excel 适合数据录入阶段，但正式分析前应转换为 CSV 格式，避免格式问题和版本兼容性问题。

2. 避免在 Excel 中使用复杂格式：不要使用合并单元格、多级表头、嵌入图表等复杂格式，这些格式在 R 中读取时会出现问题。

3. 每个工作表只存储一种数据类型：不要在同一个工作表中混合存储不同结构的数据（如上半部分是样地信息，下半部分是物种列表），这会导致读取困难。

4. 使用标准的列名：列名应简洁、无空格、无特殊字符，使用下划线分隔单词（如 plot_id、soil_ph），避免使用中文列名（虽然 R 支持，但在某些情况下可能出现编码问题）。

5. 定期备份原始 Excel 文件：在 R 中读取和清洗数据后，原始 Excel 文件应妥善保存，不要覆盖或删除，以便日后追溯数据来源。

8 常见问题处理:

9 1. 日期格式问题

data\(date <- as.Date(data\)date, origin = “1899-12-30”)

10 2. 数值被识别为字符(因为包含非数值字符)

data\(pH <- as.numeric(data\)pH)

11 3. 空单元格被识别为 NA

data <- data %>% filter(!is.na(plot_id))

**运行结果说明**:`read_excel()` 会返回一个 tibble 对象,自动识别数据类型。如果 Excel 文件包含多个工作表,默认读取第一个。`excel_sheets()` 会返回所有工作表的名称向量。读取时如果遇到合并单元格,R 会将合并区域的值填充到第一个单元格,其余单元格为 NA。日期在 Excel 中存储为数字(距离 1900-01-01 的天数),读取后可能需要转换为 R 的 Date 类型。

**生态学应用案例**:在马尾松林土壤研究中,实验室测定数据通常由技术人员录入 Excel,一个文件可能包含多个工作表(如"土壤理化性质"、"土壤酶活性"、"微生物生物量")。以下是典型的读取和整合流程:
```r
library(readxl)
library(dplyr)

# 1. 查看文件中的所有工作表
sheets <- excel_sheets("data/马尾松林土壤数据.xlsx")
print(sheets)  # 输出: "土壤理化性质" "土壤酶活性" "微生物生物量"

# 2. 读取各个工作表
soil_chem <- read_excel("data/马尾松林土壤数据.xlsx",
                        sheet = "土壤理化性质",
                        skip = 1)  # 跳过标题行

soil_enzyme <- read_excel("data/马尾松林土壤数据.xlsx",
                          sheet = "土壤酶活性",
                          skip = 1)

soil_microbe <- read_excel("data/马尾松林土壤数据.xlsx",
                           sheet = "微生物生物量",
                           skip = 1)

# 3. 检查数据结构
str(soil_chem)

# 4. 数据清洗:统一列名、移除空行
soil_chem <- soil_chem %>%
  rename(plot_id = 样地编号, pH = pH值, SOC = 有机碳) %>%
  filter(!is.na(plot_id))

# 5. 合并多个工作表(按样地编号)
soil_all <- soil_chem %>%
  left_join(soil_enzyme, by = "plot_id") %>%
  left_join(soil_microbe, by = "plot_id")

# 6. 保存为 CSV 格式(便于后续分析和共享)
write_csv(soil_all, "data/soil_data_clean.csv")

这个流程展示了如何从一个包含多个工作表的 Excel 文件中读取数据,进行清洗和整合,最后保存为标准的 CSV 格式。在实际研究中,建议将原始 Excel 文件作为”数据录入版本”保留,将清洗后的 CSV 文件作为”分析版本”使用,所有数据处理步骤都用 R 脚本记录,确保分析的透明性和可重复性。

11.0.1 读取其他格式

除了 CSV 和 Excel，生态学研究中还会遇到其他数据格式。例如，Tab 分隔的文本文件（.txt 或 .tsv）常用于基因组学和生物信息学数据；R 原生格式（.rds 和 .RData）适合保存中间结果和复杂对象；SPSS、SAS、Stata 等统计软件的数据文件也可能需要读取。掌握这些格式的读取方法，可以让你更灵活地处理不同来源的数据，避免手动转换格式带来的错误和时间浪费。

# Tab 分隔文件
library(readr)
data <- read_tsv("data/results.txt")

# R 原生格式（推荐用于保存中间结果）
saveRDS(data, "data/processed.rds")
data <- readRDS("data/processed.rds")

# 保存和加载多个对象
save(data1, data2, data3, file = "data/workspace.RData")
load("data/workspace.RData")

# SPSS 文件
library(haven)
data <- read_sav("data/survey.sav")

# Stata 文件
data <- read_dta("data/panel.dta")

# SAS 文件
data <- read_sas("data/experiment.sas7bdat")

# JSON 文件（常用于网络数据）
library(jsonlite)
data <- fromJSON("data/api_response.json")

运行结果说明：read_tsv() 与 read_csv() 类似，只是分隔符从逗号变成了制表符。saveRDS() 和 readRDS() 是 R 原生的序列化方法，可以保存任何 R 对象（包括模型、图形、列表等），而且保存的文件通常比 CSV 小得多。haven 包专门用于读取 SPSS、Stata、SAS 等统计软件的数据文件，这在需要复现他人分析或整合历史数据时非常有用。

生态学应用案例：在马尾松林长期定位监测中，我们可能需要整合多年的历史数据。早期的数据可能保存在 SPSS 或 Excel 中，近期的数据可能是 CSV 格式，而中间处理步骤生成的复杂对象（如群落多样性分析的 vegan 对象、系统发育树、空间自相关模型等）则适合用 .rds 格式保存。使用 haven 包读取历史 SPSS 数据，用 readr 读取近期 CSV 数据，用 readRDS() 加载中间结果，我们可以构建一个完整的数据分析流程，避免重复计算，提高分析效率。此外，如果需要从在线数据库（如 GBIF、WorldClim）获取物种分布或气候数据，这些数据通常以 JSON 或 GeoJSON 格式提供，掌握 jsonlite 包的使用可以让你直接在 R 中处理这些数据，无需手动下载和转换。因此，熟悉多种数据格式的读取方法，是开展现代生态学研究、实现数据整合与可重复性分析的重要技能。

11.0.2 文件路径

• 使用 / 而不是 \ 作为路径分隔符
• 使用相对路径而不是绝对路径，这样别人也能运行你的代码
• 把数据文件放在项目的 data/ 文件夹中 :::

11.1 基本统计函数

# 生成示例数据
set.seed(42)  # 设置随机种子，确保可重复
soil_ph <- rnorm(50, mean = 6.5, sd = 0.8)  # 50 个土壤 pH 值

# 描述性统计
mean(soil_ph)          # 平均值
median(soil_ph)        # 中位数
sd(soil_ph)            # 标准差
var(soil_ph)           # 方差
range(soil_ph)         # 范围
quantile(soil_ph)      # 四分位数

# 综合统计
summary(soil_ph)

11.2 基本绘图

R 内置了基本绘图功能，后续课程会详细讲解 ggplot2：

set.seed(42)
soil_ph <- rnorm(50, mean = 6.5, sd = 0.8)

# 直方图
hist(soil_ph, main = "土壤 pH 分布", xlab = "pH", col = "lightblue")

# 箱线图
boxplot(soil_ph, main = "土壤 pH", ylab = "pH")

11.2.1 获取帮助

?mean              # 查看函数帮助
help(read.csv)     # 同上
??regression       # 搜索关键词

11.3 实践：生态学数据集探索

让我们用一个真实的生态学数据集来练习。R 内置了 iris（鸢尾花）数据集：

# 加载数据
data(iris)

# 查看前 6 行
head(iris)

# 数据结构
str(iris)

# 描述性统计
summary(iris)

# 按物种分组统计花瓣长度
# tapply(数据, 分组变量, 函数) —— 对每组分别计算
tapply(iris$Petal.Length, iris$Species, mean)

# 简单散点图
plot(iris$Sepal.Length, iris$Petal.Length,
     col = as.numeric(iris$Species),
     pch = 16,
     xlab = "花萼长度 (cm)",
     ylab = "花瓣长度 (cm)",
     main = "鸢尾花花萼与花瓣长度关系")
legend("topleft", levels(iris$Species),
       col = 1:3, pch = 16, cex = 0.8)

11.4 列表（list）

列表是 R 中最灵活的数据结构，可以包含不同类型的元素：

# 创建列表
field_survey <- list(
  site = "广西大学林场",
  date = "2024-03-15",
  species_count = 15,
  dominant_species = c("马尾松", "杉木", "桉树"),
  coordinates = c(lat = 22.8, lon = 108.3),
  has_water = TRUE
)

# 访问列表元素
field_survey$site              # 用 $ 访问
field_survey[[1]]              # 用 [[ ]] 访问
field_survey[["species_count"]]  # 用名称访问

# 查看列表结构
str(field_survey)

11.4.1 列表的实际应用

在生态学研究中，列表常用于存储复杂的分析结果：

# 线性回归返回的就是列表
model <- lm(Petal.Length ~ Sepal.Length, data = iris)
class(model)
names(model)  # 查看列表包含哪些元素

# 提取回归系数
model$coefficients

# 提取拟合值
head(model$fitted.values)

11.5 矩阵（matrix）

矩阵是二维数组，所有元素必须是同一类型：

# 创建矩阵
soil_nutrients <- matrix(
  c(12.5, 8.3, 15.7,   # 第一列：氮含量
    2.1, 1.8, 2.5,     # 第二列：磷含量
    8.5, 6.2, 9.1),    # 第三列：钾含量
  nrow = 3,
  ncol = 3,
  byrow = FALSE
)

# 添加行名和列名
rownames(soil_nutrients) <- c("样地1", "样地2", "样地3")
colnames(soil_nutrients) <- c("N", "P", "K")
soil_nutrients

# 矩阵运算
soil_nutrients * 10        # 每个元素乘以 10
rowSums(soil_nutrients)    # 每行求和
colMeans(soil_nutrients)   # 每列求平均

11.5.1 矩阵的深入理解

矩阵在生态学中广泛应用于群落数据分析（物种-样地矩阵）和多元统计。

# 创建物种-样地矩阵（species-site matrix）
# 行：样地，列：物种，值：个体数
set.seed(2024)
community_matrix <- matrix(
  rpois(20, lambda = 5),  # 20个随机计数值
  nrow = 4,
  ncol = 5
)
rownames(community_matrix) <- paste0("样地", 1:4)
colnames(community_matrix) <- c("马尾松", "杉木", "桉树", "樟树", "楠木")
community_matrix

# 矩阵转置（transpose）
t(community_matrix)  # 转换为物种-样地矩阵

# 矩阵乘法
# 计算样地间的相似性矩阵（基于物种组成）
similarity <- community_matrix %*% t(community_matrix)
similarity

# 矩阵的行列操作
rbind(community_matrix, 样地5 = c(3, 7, 2, 5, 4))  # 添加行
cbind(community_matrix, 枫香 = c(6, 4, 8, 5))      # 添加列

11.5.2 生态学实例：物种丰富度计算

# 计算每个样地的物种丰富度（species richness）
species_richness <- apply(community_matrix, 1, function(x) sum(x > 0))
cat("各样地物种丰富度：\n")
print(species_richness)

# 计算每个物种的出现频率
species_frequency <- apply(community_matrix, 2, function(x) sum(x > 0) / nrow(community_matrix))
cat("\n各物种出现频率：\n")
print(round(species_frequency, 2))

# 计算样地的总个体数
total_abundance <- rowSums(community_matrix)
cat("\n各样地总个体数：\n")
print(total_abundance)

11.6 条件语句

条件语句让程序根据不同情况执行不同操作：

# if-else 语句
ph_value <- 5.5

if (ph_value < 5.5) {
  soil_type <- "强酸性"
} else if (ph_value < 6.5) {
  soil_type <- "酸性"
} else if (ph_value < 7.5) {
  soil_type <- "中性"
} else {
  soil_type <- "碱性"
}

soil_type

11.6.1 向量化条件判断

对于向量，使用 ifelse() 更高效：

ph_values <- c(5.2, 6.8, 7.1, 4.9, 6.3)

# ifelse(条件, 真时的值, 假时的值)
soil_types <- ifelse(ph_values < 6.5, "酸性", "中性或碱性")
soil_types

# 多条件判断用 case_when（需要 dplyr 包）
library(dplyr)
soil_types <- case_when(
  ph_values < 5.5 ~ "强酸性",
  ph_values < 6.5 ~ "酸性",
  ph_values < 7.5 ~ "中性",
  TRUE ~ "碱性"
)
soil_types

11.7 循环

循环用于重复执行相同操作：

11.7.1 for 循环

for 循环是编程中最基础、最常用的控制结构之一,用于对一组元素逐个执行相同的操作。在 R 中,for 循环的基本语法是 for (变量 in 序列) { 操作 },其中”变量”会依次取”序列”中的每个值,然后执行大括号内的操作。for 循环的核心价值在于”批量处理”:当你需要对多个样地、多个物种、多个文件执行相同的分析步骤时,for 循环可以让你避免重复编写代码,提高效率并减少错误。在生态学数据分析中,for 循环常用于批量读取数据文件、逐个样地计算指标、逐层分析土壤剖面数据等场景。理解 for 循环的关键是认识到:循环变量在每次迭代时都会更新,你可以用它来访问不同的数据、构造不同的文件名,或者累积计算结果。虽然 R 提倡向量化编程以提高性能,但 for 循环在逻辑清晰、易于理解方面仍有独特优势,尤其适合初学者和复杂的多步骤操作。

# 基本 for 循环:遍历元素
species <- c("马尾松", "杉木", "桉树")

for (sp in species) {
  print(paste("正在处理:", sp))
}

# 用索引循环:更灵活的访问方式
for (i in 1:length(species)) {
  print(paste("第", i, "个物种是:", species[i]))
}

# 累积计算:计算累积和
numbers <- c(5, 10, 15, 20)
cumulative_sum <- 0
for (num in numbers) {
  cumulative_sum <- cumulative_sum + num
  cat("当前累积和:", cumulative_sum, "\n")
}

# 嵌套循环:处理二维数据
plots <- c("P1", "P2")
layers <- c("0-10cm", "10-20cm")
for (plot in plots) {
  for (layer in layers) {
    cat("样地", plot, "土层", layer, "\n")
  }
}

运行结果说明:第一个循环会依次输出”正在处理: 马尾松”、“正在处理: 杉木”、“正在处理: 桉树”。第二个循环使用索引 i,可以同时访问元素位置和元素值,这在需要知道”第几个”时非常有用。累积计算示例展示了如何在循环中更新变量,最终 cumulative_sum 会得到 50。嵌套循环会生成所有样地-土层的组合,共输出 4 行。

生态学应用案例:在马尾松林土壤研究中,假设你有 20 个样地,每个样地采集了 4 层土壤(0-10cm、10-20cm、20-30cm、30-40cm),每层测定了 pH、有机碳、全氮等指标。如果要批量计算每层土壤的碳氮比,可以用 for 循环遍历所有样地和土层:

for (plot_id in 1:20) {
  for (layer in 1:4) {
    cn_ratio <- soil_data[plot_id, layer, "SOC"] / soil_data[plot_id, layer, "TN"]
    cat("样地", plot_id, "第", layer, "层 C/N =", round(cn_ratio, 2), "\n")
  }
}

这种嵌套循环结构清晰地表达了”对每个样地的每一层”的逻辑,便于理解和调试。虽然向量化方法更快,但在数据结构复杂、需要条件判断或中间输出时,for 循环往往更直观。

11.7.2 while 循环

while 循环是另一种重要的循环结构,与 for 循环的主要区别在于:for 循环通常用于”已知次数”的重复操作,而 while 循环用于”满足条件就继续”的场景。while 循环的语法是 while (条件) { 操作 },只要条件为真(TRUE),就会一直执行大括号内的操作,直到条件变为假(FALSE)才停止。这种特性使得 while 循环特别适合处理”不确定何时结束”的任务,例如迭代算法收敛、文件逐行读取直到结束、模拟过程直到满足某个阈值等。在生态学研究中,while 循环常用于模拟种群动态(直到种群稳定)、迭代优化模型参数(直到误差足够小)、或者处理不规则数据(逐行读取直到遇到特定标记)。使用 while 循环时需要特别注意:必须确保循环条件最终会变为假,否则会陷入”无限循环”,导致程序卡死。因此,循环体内通常需要包含能够改变条件的语句,例如计数器递增、累积误差更新等。

# 基本 while 循环:计数到 5
count <- 1
while (count <= 5) {
  print(paste("计数:", count))
  count <- count + 1
}

# 条件判断:累积和直到超过阈值
set.seed(2024)
total <- 0
iteration <- 0
threshold <- 50
while (total < threshold) {
  value <- rnorm(1, mean = 10, sd = 3)
  total <- total + value
  iteration <- iteration + 1
  cat("第", iteration, "次迭代,当前累积:", round(total, 2), "\n")
}
cat("达到阈值,共迭代", iteration, "次\n")

# 实用案例:迭代计算直到收敛
x <- 100  # 初始值
tolerance <- 0.01  # 收敛容差
max_iter <- 100
iter <- 0
converged <- FALSE

while (!converged && iter < max_iter) {
  x_new <- sqrt(x)  # 迭代公式
  if (abs(x_new - x) < tolerance) {
    converged <- TRUE
  }
  x <- x_new
  iter <- iter + 1
}

if (converged) {
  cat("收敛于:", round(x, 4), ",迭代", iter, "次\n")
} else {
  cat("未收敛,达到最大迭代次数\n")
}

运行结果说明:第一个循环会输出 1 到 5 的计数。第二个示例模拟了”累积随机值直到超过 50”的过程,每次迭代生成一个随机数并累加,输出显示了达到阈值所需的迭代次数(每次运行结果可能不同)。第三个示例展示了数值迭代收敛的典型模式:不断更新 x 的值,直到新旧值之差小于容差,或者达到最大迭代次数。这种”双重退出条件”是 while 循环的常见模式,既保证收敛时及时停止,又防止无限循环。

生态学应用案例:在马尾松林土壤碳氮循环模拟中,while 循环常用于模拟动态过程。例如,模拟土壤有机碳分解直到达到平衡状态:

SOC <- 50  # 初始土壤有机碳 (g/kg)
decomp_rate <- 0.05  # 分解速率
input_rate <- 2.0  # 凋落物输入速率
tolerance <- 0.1
year <- 0

while (abs(SOC * decomp_rate - input_rate) > tolerance) {
  SOC <- SOC - SOC * decomp_rate + input_rate
  year <- year + 1
  if (year %% 10 == 0) cat("第", year, "年, SOC =", round(SOC, 2), "\n")
  if (year > 500) break  # 防止无限循环
}
cat("平衡状态: SOC =", round(SOC, 2), "g/kg, 用时", year, "年\n")

这个例子模拟了土壤碳库在分解和输入作用下逐渐达到平衡的过程,while 循环会一直运行直到分解量与输入量接近相等。这种动态模拟在生态系统建模中非常常见,while 循环提供了清晰的逻辑表达。

Warning循环效率

在 R 中，循环通常比向量化操作慢。能用向量化就不要用循环：

# 慢：用循环
x <- 1:1000
result <- numeric(1000)
for (i in 1:1000) {
  result[i] <- x[i] * 2
}

# 快：向量化
result <- x * 2

11.8 自定义函数

函数让你把重复的代码封装起来，提高代码复用性：

# 定义函数
calculate_cv <- function(x) {
  # 计算变异系数 (Coefficient of Variation)
  cv <- (sd(x) / mean(x)) * 100
  return(cv)
}

# 使用函数
heights <- c(12.5, 8.3, 15.7, 10.2, 9.8)
calculate_cv(heights)

11.9 函数编写基础与向量化编程

函数是 R 编程的核心，理解函数编写和向量化编程能大幅提升代码效率。

11.9.1 函数的基本结构

函数是组织代码的基本单元，将一段可以重复使用的逻辑封装起来，给它一个名字，之后通过函数名调用。良好的函数设计可以让代码更易读、易维护、易测试。在 R 中，函数也是一种对象，可以作为参数传递给其他函数，这也是 R 作为函数式编程语言的强大之处。在生态学数据分析中，我们通常将数据清洗、指标计算、模型拟合等操作封装成函数，这样不仅代码更简洁，也便于在不同项目间复用。

# 函数的基本结构
my_function <- function(arguments) {
  result <- some_calculation
  return(result)
}

# 实用函数：计算标准误（Standard Error）
se <- function(x, na.rm = FALSE) {
  if (na.rm) x <- x[!is.na(x)]
  n <- length(x)
  if (n == 0) return(NA)
  sd(x) / sqrt(n)
}

# 测试函数
data <- c(12.5, 8.3, 15.7, 10.2, 9.8, NA)
se(data, na.rm = TRUE)

# 返回多个值：返回列表
summary_stats <- function(x, na.rm = FALSE) {
  if (na.rm) x <- x[!is.na(x)]
  list(
    n = length(x),
    mean = mean(x),
    sd = sd(x),
    se = sd(x) / sqrt(length(x)),
    min = min(x),
    max = max(x)
  )
}
summary_stats(c(12.5, 8.3, 15.7, 10.2, 9.8))

运行结果说明：se() 函数计算标准误，接受两个参数：数值向量 x 和逻辑值 na.rm（是否移除缺失值）。return() 语句用于显式返回结果，但也可以省略——R 会自动返回最后一个表达式的值。summary_stats() 函数返回一个命名列表，包含样本量、平均值、标准差、标准误、最小值和最大值，这种”返回多个值”的模式在数据分析中非常实用。

生态学应用案例：在马尾松林土壤研究中，我们经常需要计算多个指标（pH、有机碳、全氮、有效磷等）的基本统计量，并为每个样地或每种处理生成统计报告。如果每次都手动编写统计代码，不仅繁琐，还容易出错。更好的做法是定义一个 calc_soil_stats() 函数，输入土壤指标向量，输出包含所有基本统计量的列表。然后，可以用 lapply() 批量计算所有指标：r soil_vars <- list(pH = soil_ph, SOC = soil_carbon, TN = total_nitrogen, AP = available_phosphorus) stats_list <- lapply(soil_vars, summary_stats)这样一行代码就能生成完整的统计报告。更进一步，我们可以将函数保存为独立的 .R 文件，在不同项目间共享复用，这也是构建个人 R 包的基础。

11.9.2 函数参数的类型

R 函数的参数可以分为必需参数（没有默认值，调用时必须提供）和可选参数（有默认值，可以省略）。理解参数类型是编写灵活函数的基础。R 还支持特殊的 ... 参数（ellipsis 或 dots），表示可以接受任意数量的额外参数，常用于将参数转发给其他函数，或者实现”传递未知参数”的功能。

# 必需参数 vs 默认参数
required_func <- function(x, y) {
  x + y
}
# required_func(1)  # 错误：缺少参数 y

default_func <- function(x, y = 10) {
  x + y
}
default_func(5)
default_func(5, 20)

# ... 参数
varargs_func <- function(x, ...) {
  args <- list(...)
  cat("主参数 x =", x, "
")
  cat("额外参数数量:", length(args), "
")
}
varargs_func(1, 2, 3, name = "马尾松")

# 参数验证
safe_divide <- function(a, b) {
  if (!is.numeric(a) || !is.numeric(b)) stop("两个参数都必须为数值型")
  if (b == 0) {
    warning("除数为零，返回 Inf")
    return(Inf)
  }
  a / b
}
safe_divide(10, 2)
safe_divide(10, 0)

运行结果说明：required_func() 的两个参数都没有默认值。default_func() 的 y 有默认值 10，因此可以只提供 x。varargs_func() 演示了 ... 的用法：它接受任意数量的额外参数并收集成列表。参数验证示例展示了如何检查输入的合法性，避免无效运算。

生态学应用案例：在马尾松林土壤数据分析中，假设你编写了一个 calc_diversity() 函数，计算 Shannon、Simpson、Pielou 等多个指数。为了让它更灵活，可以设计成这样：r calc_diversity <- function(community_matrix, index = "all", na.rm = FALSE, ...) { # index: 指定计算哪个指数，或者 "all" 计算全部 # na.rm: 是否移除缺失值 # ...: 传递给各指数计算函数的额外参数 }用户可以这样调用：calc_diversity(df)（计算所有指数）、calc_diversity(df, index = "shannon")（只计算 Shannon）、calc_diversity(df, na.rm = TRUE)（移除缺失值）。这种设计让函数既简单易用（默认参数满足大多数情况），又灵活强大。

11.9.3 函数的作用域（Scope）

作用域（Scope）决定了在哪里可以访问哪个变量。R 中有全局变量和局部变量之分：在函数内部创建的变量是局部变量，只能在函数内部使用；函数外部创建的变量是全局变量，在整个代码中都可以访问。理解作用域的规则，可以帮助你避免意想不到的错误，比如在函数内部意外修改了全局变量。R 还支持通过 <<- 运算符在函数内部修改全局变量，但这种用法需要谨慎。

# 全局变量 vs 局部变量
global_var <- 10

scope_demo <- function() {
  local_var <- 20
  cat("函数内部：global_var =", global_var, "
")
  cat("函数内部：local_var =", local_var, "
")
}
scope_demo()
cat("函数外部：global_var =", global_var, "
")

# <<- 修改全局变量（谨慎使用）
counter <- 0
increment <- function() {
  counter <<- counter + 1
  return(counter)
}
cat("计数器:", increment(), "
")
cat("计数器:", increment(), "
")
cat("全局变量 counter =", counter, "
")

运行结果说明：global_var 在函数外部创建，在函数内部可以读取；local_var 在函数内部创建，只能在函数内部访问。<<- 运算符将 counter 的修改”穿透”到全局环境，因此每次调用 increment() 都会增加全局计数器的值。<<- 虽然强大，但会破坏函数的纯粹性，因此应该避免滥用。

生态学应用案例：在马尾松林长期监测数据的分析中，我们可能需要在多个函数间共享样地元数据。一个好的做法是：将这些元数据保存在一个全局列表中，然后在各个函数中直接访问，而不需要每次都作为参数传递。但是要注意，如果多个函数都修改全局变量，代码会变得难以调试。更好的做法是通过参数传递，或者使用 R6 等面向对象框架来实现有封装的状态管理。

11.9.4 向量化编程

向量化是 R 的核心优势，避免使用循环能大幅提升性能。

# ❌ 不好的实践：使用循环
x <- 1:1000
result <- numeric(1000)
system.time({
  for (i in 1:1000) {
    result[i] <- x[i]^2
  }
})

# ✅ 好的实践：向量化
system.time({
  result <- x^2
})

# 向量化函数示例
# 计算多个样地的 Shannon 多样性指数
shannon_diversity <- function(abundance_matrix) {
  # abundance_matrix: 行为样地，列为物种
  apply(abundance_matrix, 1, function(row) {
    p <- row / sum(row)  # 相对丰度
    p <- p[p > 0]        # 移除 0 值
    -sum(p * log(p))     # Shannon 指数
  })
}

# 测试
set.seed(2024)
community <- matrix(rpois(20, lambda = 5), nrow = 4, ncol = 5)
shannon_diversity(community)

11.9.5 生态学实例：批量计算土壤指标

# 定义函数：计算土壤质量指数
soil_quality_index <- function(ph, organic_matter, nitrogen) {
  # 标准化各指标（0-1）
  ph_score <- ifelse(ph >= 6 & ph <= 7, 1, 
                     ifelse(ph < 6, (ph - 4) / 2, (8 - ph) / 2))
  om_score <- pmin(organic_matter / 5, 1)  # 有机质 > 5% 得满分
  n_score <- pmin(nitrogen / 2, 1)          # 氮含量 > 2% 得满分
  
  # 综合指数（加权平均）
  sqi <- 0.4 * ph_score + 0.3 * om_score + 0.3 * n_score
  
  # 返回详细结果
  data.frame(
    ph = ph,
    organic_matter = organic_matter,
    nitrogen = nitrogen,
    ph_score = round(ph_score, 2),
    om_score = round(om_score, 2),
    n_score = round(n_score, 2),
    sqi = round(sqi, 2),
    quality = case_when(
      sqi >= 0.8 ~ "优",
      sqi >= 0.6 ~ "良",
      sqi >= 0.4 ~ "中",
      TRUE ~ "差"
    )
  )
}

# 批量计算多个样地
set.seed(2024)
soil_data <- data.frame(
  plot_id = 1:10,
  ph = rnorm(10, mean = 6.5, sd = 0.8),
  organic_matter = rnorm(10, mean = 3.5, sd = 1.2),
  nitrogen = rnorm(10, mean = 1.5, sd = 0.5)
)

# 使用向量化函数
soil_results <- soil_quality_index(
  soil_data$ph,
  soil_data$organic_matter,
  soil_data$nitrogen
)

head(soil_results)

11.9.6 匿名函数（Anonymous Functions）

匿名函数,顾名思义,就是没有名字的函数。在 R 编程中,我们通常用 function_name <- function(x) { ... } 的方式定义一个有名字的函数,以便后续多次调用。但有时候,我们只需要在某个特定位置使用一次函数,不需要给它起名字并保存到环境中,这时就可以使用匿名函数。匿名函数最常见的应用场景是与 apply 系列函数、lapply、sapply、map 等高阶函数配合使用,用于对数据的每个元素执行临时的、一次性的操作。匿名函数的优势在于代码简洁、逻辑集中,不会污染全局环境(不会创建不必要的函数对象)。在生态学数据分析中,匿名函数常用于批量计算指标、数据转换、条件筛选等场景。R 4.1+ 还引入了更简洁的匿名函数语法 \(x) x^2,进一步提升了代码的可读性。理解匿名函数是掌握函数式编程思想的重要一步,它让你能够更灵活地处理复杂的数据操作。

# 方法1:传统的命名函数
square <- function(x) x^2
sapply(1:5, square)

# 方法2:匿名函数(不需要命名)
sapply(1:5, function(x) x^2)

# 方法3:R 4.1+ 简化语法(推荐)
sapply(1:5, \(x) x^2)

# 实际应用1:计算多个统计量
data <- c(12.5, 8.3, 15.7, 10.2, 9.8)
lapply(list(mean, sd, median, IQR), function(f) f(data))

# 实际应用2:对列表的每个元素执行复杂操作
soil_data_list <- list(
  plot1 = c(6.5, 5.8, 6.2),
  plot2 = c(7.1, 6.9, 7.3),
  plot3 = c(5.5, 5.2, 5.8)
)

# 计算每个样地的 pH 统计量
lapply(soil_data_list, function(x) {
  c(mean = mean(x),
    sd = sd(x),
    cv = sd(x) / mean(x) * 100)
})

# 实际应用3:条件转换
species_list <- c("Pinus massoniana", "Cunninghamia lanceolata", "Eucalyptus")
sapply(species_list, function(sp) {
  if (grepl("Pinus", sp)) "松属"
  else if (grepl("Cunninghamia", sp)) "杉属"
  else "其他"
})

# 实际应用4:与 dplyr 配合使用
library(dplyr)
iris %>%
  group_by(Species) %>%
  summarise(across(starts_with("Sepal"),
                   list(mean = \(x) mean(x),
                        cv = \(x) sd(x) / mean(x) * 100)))

# 实际应用5:嵌套匿名函数
data_matrix <- matrix(1:12, nrow = 3, ncol = 4)
apply(data_matrix, 1, function(row) {
  sapply(row, function(x) x^2)
})

# 对比:命名函数 vs 匿名函数
# 命名函数:适合多次使用
calculate_cv <- function(x) sd(x) / mean(x) * 100
cv1 <- calculate_cv(c(10, 12, 15))
cv2 <- calculate_cv(c(20, 22, 25))

# 匿名函数:适合一次性使用
lapply(list(c(10, 12, 15), c(20, 22, 25)),
       \(x) sd(x) / mean(x) * 100)

运行结果说明:sapply(1:5, function(x) x^2) 会对 1 到 5 的每个数字求平方,返回向量 c(1, 4, 9, 16, 25)。lapply() 返回列表,sapply() 尝试简化为向量或矩阵。匿名函数可以包含多行代码和复杂逻辑,只要用大括号 {} 包裹即可。R 4.1+ 的 \(x) 语法与 function(x) 完全等价,只是更简洁。

生态学应用案例:在马尾松林土壤研究中,假设你有多个样地的土壤剖面数据,需要对每个样地的每一层计算碳氮比,并判断是否符合”碳氮比 > 15”的标准:

library(purrr)

# 土壤剖面数据(列表格式)
soil_profiles <- list(
  plot1 = data.frame(
    layer = c("0-10", "10-20", "20-30"),
    SOC = c(18.5, 15.2, 12.8),
    TN = c(1.2, 1.0, 0.9)
  ),
  plot2 = data.frame(
    layer = c("0-10", "10-20", "20-30"),
    SOC = c(22.3, 18.7, 15.5),
    TN = c(1.5, 1.3, 1.1)
  )
)

# 使用匿名函数批量计算碳氮比并分类
results <- map(soil_profiles, \(df) {
  df %>%
    mutate(
      CN_ratio = SOC / TN,
      CN_class = ifelse(CN_ratio > 15, "高碳氮比", "低碳氮比")
    )
})

# 提取每个样地的平均碳氮比
mean_CN <- map_dbl(results, \(df) mean(df$CN_ratio))
print(mean_CN)

# 统计各样地高碳氮比土层的数量
high_CN_count <- map_int(results, \(df) sum(df$CN_class == "高碳氮比"))
print(high_CN_count)

这个例子展示了匿名函数在处理嵌套数据结构时的强大能力。通过 map() 和匿名函数的组合,我们可以用简洁的代码对每个样地的数据框执行相同的操作,避免了繁琐的 for 循环。匿名函数让代码更加紧凑,逻辑更加清晰,是现代 R 编程的重要技巧。

11.9.7 函数式编程：purrr 包

library(purrr)

# map 系列函数：更安全的 apply 替代品
data_list <- list(
  plot1 = c(12.5, 8.3, 15.7),
  plot2 = c(10.2, 9.8, 11.3),
  plot3 = c(14.1, 13.5, 12.8)
)

# map_dbl：返回数值向量
map_dbl(data_list, mean)

# map_df：返回数据框
map_df(data_list, ~data.frame(mean = mean(.x), sd = sd(.x)), .id = "plot")

# 嵌套数据处理
iris %>%
  split(.$Species) %>%
  map_df(~summarise(.x, 
                    mean_sepal = mean(Sepal.Length),
                    mean_petal = mean(Petal.Length)),
         .id = "Species")

11.10 调试技巧

编写代码时难免遇到错误，掌握调试技巧能快速定位和解决问题。

11.10.1 常见错误类型

在 R 编程中，新手经常遇到各种错误。理解这些错误的类型和原因，可以帮助你更快地定位和解决问题。常见的错误类型包括：语法错误（缺少括号、引号不匹配）、对象未找到（变量名拼写错误或未赋值）、类型错误（对错误类型的对象执行操作）、索引越界（访问不存在的元素）、参数错误（缺少必需参数或参数类型不对）等。

# 1. 语法错误
x <- c(1, 2, 3  # 缺少右括号

# 2. 对象未找到
print(undefined_variable)

# 3. 类型错误
"text" + 5

# 4. 索引越界
x <- c(1, 2, 3)
x[5]   # 返回 NA
x[0]   # 返回 numeric(0)

# 5. 参数错误
mean("text")

# 6. 缺失值未处理
x <- c(1, 2, NA, 4)
mean(x)
# 需要：mean(x, na.rm = TRUE)

# 7. 赋值错误：<- vs ==
if (x = 5) print("错误")  # 应该是 x == 5

运行结果说明：这些错误覆盖了 R 编程中最常见的问题。特别要注意：R 的索引从 1 开始（不是 0），且 x[5] 不会报错，而是返回 NA。在 if 条件中，必须使用 == 进行比较，使用 = 会变成赋值。

生态学应用案例：在马尾松林土壤数据分析中，这些常见错误经常出现。例如：- 把土壤 pH 的缺失值记录为 “NA”（字符），导致后续计算失败 - 假设有 20 个样地，但代码中写成了 for (i in 1:21)，试图访问不存在的第 21 个样地 - 计算平均胸径时，没有设置 na.rm = TRUE，导致结果为 NA

避免这些错误的方法包括：1) 养成在代码开头检查数据结构的习惯（用 str()）；2) 对缺失值使用 na.rm = TRUE 处理；3) 在循环前检查向量长度。

11.10.2 调试函数

1. print() 和 cat() 调试

# 最简单的调试方法：打印中间结果
my_function <- function(x) {
  cat("输入值:", x, "\n")
  result <- x^2
  cat("计算结果:", result, "\n")
  return(result)
}

my_function(5)

2. browser() 交互式调试

# browser() 会暂停函数执行，进入交互式调试模式
debug_function <- function(x, y) {
  z <- x + y
  browser()  # 在这里暂停
  result <- z * 2
  return(result)
}

debug_function(3, 5)
# 进入调试模式后可以：
# - 输入变量名查看值
# - n: 执行下一行
# - c: 继续执行到结束
# - Q: 退出调试

3. debug() 和 undebug()

# 对整个函数启用调试
my_function <- function(x) {
  result <- x^2
  return(result)
}

debug(my_function)    # 启用调试
my_function(5)        # 每次调用都会进入调试模式
undebug(my_function)  # 关闭调试

4. traceback() 追踪错误

# 当函数报错后，使用 traceback() 查看调用栈
func1 <- function(x) func2(x)
func2 <- function(x) func3(x)
func3 <- function(x) x + "text"  # 这里会报错

func1(5)  # 报错
traceback()  # 查看错误发生在哪个函数

5. tryCatch() 错误处理

# 捕获错误并优雅处理
safe_mean <- function(x) {
  tryCatch(
    {
      result <- mean(x)
      return(result)
    },
    error = function(e) {
      cat("错误:", e$message, "\n")
      return(NA)
    },
    warning = function(w) {
      cat("警告:", w$message, "\n")
      return(mean(x, na.rm = TRUE))
    }
  )
}

safe_mean(c(1, 2, 3))      # 正常
safe_mean("text")          # 捕获错误
safe_mean(c(1, 2, NA))     # 捕获警告

6. stopifnot() 和 stop() 参数检查

# 使用 stopifnot() 检查前置条件
calculate_diversity <- function(abundance) {
  stopifnot(
    is.numeric(abundance),           # 必须是数值型
    all(abundance >= 0),             # 必须非负
    length(abundance) > 0            # 不能为空
  )
  
  p <- abundance / sum(abundance)
  -sum(p * log(p + 1e-10))  # 避免 log(0)
}

# 测试
calculate_diversity(c(10, 5, 3))  # 正常
# calculate_diversity(c(-1, 5, 3))  # 报错：abundance 必须非负

# 使用 stop() 自定义错误信息
calculate_diversity2 <- function(abundance) {
  if (!is.numeric(abundance)) {
    stop("abundance 必须是数值型向量")
  }
  if (any(abundance < 0)) {
    stop("abundance 不能包含负值")
  }
  
  p <- abundance / sum(abundance)
  -sum(p * log(p + 1e-10))
}

7. warning() 发出警告

# 使用 warning() 提示潜在问题
calculate_cv <- function(x) {
  if (any(is.na(x))) {
    warning("数据包含 NA 值，已自动移除")
    x <- x[!is.na(x)]
  }
  
  if (length(x) < 3) {
    warning("样本量过小 (n < 3)，结果可能不可靠")
  }
  
  cv <- (sd(x) / mean(x)) * 100
  return(cv)
}

calculate_cv(c(1, 2, NA))  # 发出警告但继续执行

11.10.3 生态学实例：调试数据处理流程

调试（Debugging）是找出并修复代码错误的过程。常见的调试策略包括：使用 print() 或 cat() 输出中间变量、检查数据结构（str()）、使用 browser() 进入交互式调试器、使用 warning() 和 message() 输出提示信息、以及使用 tryCatch() 捕获错误并继续执行。

process_soil_data <- function(data) {
  # 步骤 1：检查输入
  if (!is.data.frame(data)) stop("输入必须是数据框")
  cat("原始数据行数:", nrow(data), "
")

  # 步骤 2：检查必需列
  required_cols <- c("plot_id", "species", "dbh", "height")
  missing_cols <- setdiff(required_cols, names(data))
  if (length(missing_cols) > 0) {
    stop(paste("缺少必需列:", paste(missing_cols, collapse = ", ")))
  }

  # 步骤 3：异常值检测
  outliers_idx <- abs(data$dbh - mean(data$dbh)) > 3 * sd(data$dbh)
  n_outliers <- sum(outliers_idx, na.rm = TRUE)
  if (n_outliers > 0) {
    warning(paste("检测到", n_outliers, "个异常值"))
  }

  # 步骤 4：计算汇总统计
  result <- aggregate(cbind(dbh, height) ~ species, data = data, FUN = mean)
  return(result)
}

# 测试
set.seed(2024)
test_data <- data.frame(
  plot_id = paste0("P", 1:20),
  species = rep(c("马尾松", "杉木"), each = 10),
  dbh = c(rnorm(10, 25, 5), rnorm(10, 18, 4)),
  height = c(rnorm(10, 15, 3), rnorm(10, 12, 2))
)
test_data$dbh[3] <- 150  # 插入异常值

tryCatch({
  result <- process_soil_data(test_data)
  print(result)
}, error = function(e) cat("错误:", conditionMessage(e), "
"))

运行结果说明：process_soil_data() 函数展示了规范的调试流程：在每个关键步骤输出中间结果，便于追踪错误发生的位置。用 tryCatch() 包装调用，可以捕获任何未处理的错误并给出友好的提示。

生态学应用案例：在马尾松林土壤研究中，野外观测数据经常存在各种问题：样地编号录入错误、测量值超出合理范围（如胸径 150cm 超过测量工具上限）、缺失值用特殊符号标记等。编写健壮的数据处理函数，可以在数据分析的早期发现这些问题，避免错误结论。

11.10.4 调试最佳实践

# ✅ 好的实践：
# 1. 使用有意义的变量名
species_richness <- calculate_richness(community_data)  # 清晰

# 2. 添加注释说明复杂逻辑
# 计算 Shannon 多样性指数: H' = -Σ(pi * ln(pi))
shannon <- -sum(p * log(p))

# 3. 分步骤调试复杂操作
step1 <- filter(data, condition)
cat("Step 1:", nrow(step1), "rows\n")
step2 <- mutate(step1, new_var = calculation)
cat("Step 2:", nrow(step2), "rows\n")

# 4. 使用 assertthat 包进行断言
library(assertthat)
assert_that(is.numeric(x))
assert_that(all(x >= 0))

# ❌ 不好的实践：
# 1. 忽略警告信息
suppressWarnings(mean(c(1, 2, "text")))  # 隐藏了真正的问题

# 2. 过度使用 try() 而不处理错误
try(risky_operation(), silent = TRUE)  # 错误被忽略

# 3. 不检查输入参数
my_function <- function(x) {
  x^2  # 如果 x 不是数值型会报错
}

11.11 apply 函数族

apply 函数族是 R 中处理数据的强大工具，可以替代大部分循环：

11.11.1 apply()：对矩阵或数组操作

# 创建示例矩阵
soil_data <- matrix(rnorm(15, mean = 10, sd = 2), nrow = 5, ncol = 3)
colnames(soil_data) <- c("N", "P", "K")
rownames(soil_data) <- paste0("样地", 1:5)
soil_data

# apply(数据, 维度, 函数)
# 维度：1 = 行，2 = 列
apply(soil_data, 1, mean)  # 每行的平均值
apply(soil_data, 2, sd)    # 每列的标准差

11.11.2 lapply()：对列表操作，返回列表

# 创建列表
measurements <- list(
  plot1 = c(12.5, 8.3, 15.7),
  plot2 = c(10.2, 9.8, 11.3),
  plot3 = c(14.1, 13.5, 12.8)
)

# 对每个元素计算平均值
lapply(measurements, mean)

# 对每个元素计算长度
lapply(measurements, length)

11.11.3 sapply()：简化版 lapply()，返回向量或矩阵

sapply() 是 lapply() 的简化版本，会自动尝试将结果简化（simplify）为向量或矩阵，而不是返回列表。这使得结果更易于阅读和后续处理，特别适合需要数值型输出的场景。

sapply() 与 lapply() 的核心区别

两者的主要区别在于返回值类型：

• lapply()：始终返回列表（list），无论函数返回什么
• sapply()：尝试简化结果为向量或矩阵，简化失败时退化为列表

# 示例数据：三个样地的土壤 pH 测量值
measurements <- list(
  plot1 = c(6.2, 6.5, 6.3, 6.4),
  plot2 = c(7.1, 7.3, 7.2),
  plot3 = c(5.8, 5.9, 6.0, 6.1, 5.7)
)

# lapply() 返回列表
result_lapply <- lapply(measurements, mean)
class(result_lapply)  # "list"
result_lapply

# sapply() 返回命名向量
result_sapply <- sapply(measurements, mean)
class(result_sapply)  # "numeric"
result_sapply

# sapply() 的结果可以直接用于数值运算
result_sapply + 1  # 所有均值加 1

返回值类型的自动判断

sapply() 根据函数返回值的长度和类型自动决定输出格式：

1. 返回向量（函数返回单个值）

# 计算每个样地的平均值 → 返回向量
sapply(measurements, mean)

# 计算每个样地的样本量 → 返回向量
sapply(measurements, length)

# 计算每个样地的标准差 → 返回向量
sapply(measurements, sd)

2. 返回矩阵（函数返回多个值）

# 计算每个样地的汇总统计 → 返回矩阵
result_matrix <- sapply(measurements, summary)
class(result_matrix)  # "matrix"
result_matrix

# 每列对应一个样地，每行对应一个统计量
# 可以方便地提取特定统计量
result_matrix["Mean", ]  # 提取所有样地的均值
result_matrix[, "plot1"]  # 提取 plot1 的所有统计量

3. 返回列表（函数返回不规则结果）

# 当结果长度不一致时，sapply() 退化为列表
irregular_function <- function(x) {
  if (length(x) > 3) return(c(mean(x), sd(x)))
  else return(mean(x))
}

result_irregular <- sapply(measurements, irregular_function)
class(result_irregular)  # "list"（因为结果长度不一致）

simplify 参数控制

sapply() 的 simplify 参数可以控制是否简化结果：

# simplify = TRUE（默认）：尝试简化
sapply(measurements, mean, simplify = TRUE)  # 返回向量

# simplify = FALSE：强制返回列表（等同于 lapply）
sapply(measurements, mean, simplify = FALSE)  # 返回列表

# simplify = "array"：尝试简化为数组
sapply(measurements, summary, simplify = "array")  # 返回矩阵

常见使用场景

场景 1：批量计算统计量

# 对 iris 数据集的所有数值列计算均值
iris_numeric <- iris[, 1:4]
sapply(iris_numeric, mean)

# 对比 lapply() 的结果（列表形式，不便于查看）
lapply(iris_numeric, mean)

场景 2：批量应用自定义函数

# 计算变异系数（CV = sd / mean * 100）
calculate_cv <- function(x) {
  sd(x) / mean(x) * 100
}

sapply(iris_numeric, calculate_cv)

场景 3：提取列表元素的特定属性

# 创建多个线性模型
models <- list(
  model1 = lm(Sepal.Length ~ Sepal.Width, data = iris),
  model2 = lm(Petal.Length ~ Petal.Width, data = iris),
  model3 = lm(Sepal.Length ~ Petal.Length, data = iris)
)

# 提取每个模型的 R²
sapply(models, function(m) summary(m)$r.squared)

# 提取每个模型的 AIC
sapply(models, AIC)

sapply() 的常见错误

错误 1：假设结果总是向量

# 当函数返回多个值时，结果是矩阵而非向量
result <- sapply(measurements, summary)
# 错误用法：result[1]（这会提取矩阵的第一个元素，而非第一个样地）
# 正确用法：result[, 1]（提取第一列，即第一个样地的所有统计量）

错误 2：忽略 NA 值

# 数据中包含 NA
measurements_na <- list(
  plot1 = c(6.2, 6.5, NA, 6.4),
  plot2 = c(7.1, 7.3, 7.2),
  plot3 = c(5.8, NA, 6.0, 6.1, 5.7)
)

# 错误：未处理 NA，结果包含 NA
sapply(measurements_na, mean)

# 正确：使用 na.rm = TRUE
sapply(measurements_na, mean, na.rm = TRUE)

错误 3：对数据框使用 sapply()

# sapply() 会将数据框的每一列作为独立元素处理
# 如果列类型不同（如数值列和字符列混合），结果可能不符合预期
mixed_df <- data.frame(
  numeric_col = 1:5,
  character_col = letters[1:5]
)

# 结果会将字符列转换为因子水平
sapply(mixed_df, class)

# 更安全的做法：先筛选数值列
sapply(mixed_df[sapply(mixed_df, is.numeric)], mean)

性能对比

对于大规模数据，sapply() 和 lapply() 的性能几乎相同（简化步骤的开销可忽略）。选择哪个主要取决于输出格式的需求：

# 性能测试（需要 microbenchmark 包）
library(microbenchmark)

large_list <- replicate(1000, rnorm(100), simplify = FALSE)

microbenchmark(
  lapply = lapply(large_list, mean),
  sapply = sapply(large_list, mean),
  times = 100
)
# 结果显示两者性能相近

生态学案例：某研究团队在分析 50 个森林样地的土壤养分数据时，需要计算每个样地的 pH、有机碳、全氮的均值和标准差。最初使用 lapply() 返回列表，导致后续提取结果时代码冗长。改用 sapply() 后，直接得到矩阵格式的结果，大幅简化了数据整理流程。该案例展示了 sapply() 在批量统计分析中的便利性。

何时使用 sapply()

场景	推荐函数	原因
需要向量或矩阵输出	sapply()	结果更简洁，便于后续处理
需要保持列表结构	lapply()	避免意外简化
结果长度不确定	lapply()	sapply() 可能产生意外的列表
需要明确控制输出类型	vapply()	更安全，需显式指定返回类型

扩展记录: 2026-04-11 | 扩展者：Clawd | 目标字数：800+

11.11.4 tapply()：分组操作

# 按物种分组计算平均花瓣长度
tapply(iris$Petal.Length, iris$Species, mean)

# 按物种分组计算标准差
tapply(iris$Petal.Length, iris$Species, sd)

# 自定义函数
tapply(iris$Petal.Length, iris$Species, function(x) {
  c(mean = mean(x), sd = sd(x), cv = sd(x)/mean(x)*100)
})

11.12 缺失值处理

真实数据常包含缺失值（NA），需要正确处理：

# 创建包含 NA 的数据
soil_ph <- c(6.5, 5.8, NA, 7.2, 6.1, NA, 5.5)

# 检测 NA
is.na(soil_ph)
sum(is.na(soil_ph))  # 统计 NA 个数

# 移除 NA
na.omit(soil_ph)

# 大多数函数需要 na.rm = TRUE 参数
mean(soil_ph)              # 返回 NA
mean(soil_ph, na.rm = TRUE)  # 正确计算

# 替换 NA
soil_ph[is.na(soil_ph)] <- mean(soil_ph, na.rm = TRUE)
soil_ph

11.12.1 完整案例分析

# 创建包含缺失值的数据框
tree_survey <- data.frame(
  plot = 1:10,
  species = c("马尾松", "杉木", "桉树", "樟树", "楠木",
              "马尾松", "杉木", "桉树", "樟树", "楠木"),
  height = c(12.5, 8.3, NA, 10.2, 9.8, 13.1, NA, 15.7, 11.3, 10.5),
  dbh = c(25.3, 15.8, 32.1, NA, 18.7, 26.5, 16.2, NA, 21.8, 19.3)
)

# 查看缺失值模式
summary(tree_survey)

# 完整案例分析（删除任何包含 NA 的行）
complete_data <- na.omit(tree_survey)
nrow(complete_data)

# 按列删除 NA
tree_survey$height[is.na(tree_survey$height)] <-
  mean(tree_survey$height, na.rm = TRUE)

11.13 字符串处理

生态学数据中经常需要处理物种名、样地编号等文本信息：

# 基础字符串操作
species_name <- "Pinus massoniana"

# 字符串长度
nchar(species_name)

# 转换大小写
toupper(species_name)
tolower(species_name)

# 字符串拼接
paste("物种:", species_name)
paste0("plot_", 1:5)  # 无分隔符拼接

# 字符串分割
strsplit(species_name, " ")

# 字符串替换
gsub("massoniana", "elliottii", species_name)

11.13.1 stringr 包（推荐）

library(stringr)

# 更直观的函数名
str_length(species_name)
str_to_upper(species_name)
str_split(species_name, " ")
str_replace(species_name, "massoniana", "elliottii")

# 检测模式
species_list <- c("Pinus massoniana", "Cunninghamia lanceolata",
                  "Eucalyptus robusta", "Pinus elliottii")
str_detect(species_list, "Pinus")  # 哪些包含 "Pinus"
str_subset(species_list, "Pinus")  # 提取包含 "Pinus" 的

11.14 日期和时间

处理野外调查日期、实验时间等：

# 创建日期
survey_date <- as.Date("2024-03-15")
class(survey_date)

# 日期运算
survey_date + 7  # 7 天后
survey_date - as.Date("2024-01-01")  # 距离年初多少天

# 提取日期成分
format(survey_date, "%Y")  # 年
format(survey_date, "%m")  # 月
format(survey_date, "%d")  # 日
format(survey_date, "%Y年%m月%d日")  # 自定义格式

11.14.1 lubridate 包（推荐）

library(lubridate)

# 更方便的日期创建
ymd("2024-03-15")
mdy("03/15/2024")
dmy("15-03-2024")

# 提取成分
date <- ymd("2024-03-15")
year(date)
month(date)
day(date)
wday(date, label = TRUE)  # 星期几

# 日期序列
seq(ymd("2024-01-01"), ymd("2024-12-31"), by = "month")

11.15 数据重塑

生态学数据常需要在宽格式和长格式之间转换：

# 宽格式数据
wide_data <- data.frame(
  plot = 1:3,
  pine = c(10, 15, 12),
  fir = c(8, 12, 10),
  eucalyptus = c(5, 8, 6)
)
wide_data

# 转换为长格式（tidyr 包）
library(tidyr)
long_data <- pivot_longer(
  wide_data,
  cols = c(pine, fir, eucalyptus),
  names_to = "species",
  values_to = "count"
)
long_data

# 长格式转回宽格式
pivot_wider(
  long_data,
  names_from = species,
  values_from = count
)

11.16 数据合并

合并不同来源的数据：

# 创建两个数据框
tree_height <- data.frame(
  plot = 1:5,
  species = c("马尾松", "杉木", "桉树", "樟树", "楠木"),
  height = c(12.5, 8.3, 15.7, 10.2, 9.8)
)

tree_diameter <- data.frame(
  plot = 1:5,
  dbh = c(25.3, 15.8, 32.1, 20.5, 18.7),
  age = c(15, 10, 12, 18, 14)
)

# 按列合并（cbind）
cbind(tree_height, tree_diameter[, -1])  # 去掉重复的 plot 列

# 按行合并（rbind）
more_trees <- data.frame(
  plot = 6:7,
  species = c("马尾松", "杉木"),
  height = c(13.2, 9.1)
)
rbind(tree_height, more_trees)

# 数据库式合并（merge）
soil_data <- data.frame(
  plot = c(1, 2, 3, 4, 5),
  ph = c(6.5, 5.8, 7.2, 6.1, 5.5)
)

merge(tree_height, soil_data, by = "plot")

11.16.1 dplyr 的 join 函数（推荐）

library(dplyr)

# 内连接：只保留两边都有的
inner_join(tree_height, soil_data, by = "plot")

# 左连接：保留左边所有行
left_join(tree_height, soil_data, by = "plot")

# 完全连接：保留所有行
full_join(tree_height, soil_data, by = "plot")

11.17 管道操作符 %>%

管道操作符（pipe operator）是 tidyverse 的核心特性，让代码更易读、更符合人类思维习惯。

11.17.1 为什么需要管道操作符？

传统 R 代码从内向外读，不符合人类思维：

# 传统写法：从内向外读，难以理解
result <- round(mean(filter(iris, Species == "setosa")$Sepal.Length), 2)

# 或者需要创建多个中间变量
setosa_data <- filter(iris, Species == "setosa")
setosa_length <- setosa_data$Sepal.Length
setosa_mean <- mean(setosa_length)
result <- round(setosa_mean, 2)

管道操作符让代码从上到下、从左到右读，符合自然思维：

library(dplyr)

# 管道写法：从上到下读，清晰明了
result <- iris %>%
  filter(Species == "setosa") %>%
  pull(Sepal.Length) %>%
  mean() %>%
  round(2)

result

11.17.2 管道操作符的工作原理

%>% 将左边的结果作为右边函数的第一个参数：

# 这两种写法等价
x %>% f()        # 管道写法
f(x)             # 传统写法

# 这两种写法等价
x %>% f(y)       # 管道写法
f(x, y)          # 传统写法

# 使用 . 指定参数位置
x %>% f(y, .)    # 等价于 f(y, x)

11.17.3 实际应用示例

library(dplyr)

# 示例1：数据筛选与汇总
iris %>%
  filter(Species != "setosa") %>%           # 筛选
  mutate(Sepal.Area = Sepal.Length * Sepal.Width) %>%  # 新建列
  group_by(Species) %>%                     # 分组
  summarise(
    mean_area = mean(Sepal.Area),
    sd_area = sd(Sepal.Area),
    n = n()
  ) %>%                                     # 汇总统计
  arrange(desc(mean_area))                  # 排序

11.17.4 生态学实例：森林样地数据处理

# 创建模拟数据
set.seed(2024)
forest_data <- data.frame(
  plot_id = rep(1:5, each = 10),
  species = rep(c("马尾松", "杉木"), 25),
  dbh = rnorm(50, mean = 20, sd = 5),
  height = rnorm(50, mean = 12, sd = 3)
) %>%
  mutate(
    dbh = round(pmax(dbh, 5), 1),
    height = round(pmax(height, 3), 1)
  )

# 使用管道进行复杂数据处理
forest_summary <- forest_data %>%
  filter(dbh > 10) %>%                      # 筛选胸径 > 10cm 的树木
  mutate(
    biomass = 0.05 * dbh^2 * height,        # 计算生物量
    size_class = case_when(
      dbh < 15 ~ "小径木",
      dbh < 25 ~ "中径木",
      TRUE ~ "大径木"
    )
  ) %>%
  group_by(plot_id, species) %>%            # 按样地和物种分组
  summarise(
    n_trees = n(),
    mean_dbh = mean(dbh),
    mean_height = mean(height),
    total_biomass = sum(biomass),
    .groups = "drop"
  ) %>%
  arrange(plot_id, desc(total_biomass))     # 排序

head(forest_summary, 10)

11.17.5 原生管道 |>（R 4.1+）

R 4.1.0 引入了原生管道操作符 |>，不需要加载额外的包：

# magrittr 管道（需要加载 dplyr 或 magrittr）
iris %>%
  filter(Species == "setosa") %>%
  head(3)

# 原生管道（R 4.1+，无需加载包）
iris |>
  filter(Species == "setosa") |>
  head(3)

两者的区别：

特性	%>% (magrittr)	|> (原生)
需要加载包	是（dplyr/magrittr）	否
占位符	.	_（R 4.2+）
性能	稍慢	稍快
功能	更丰富	基础功能

建议：新项目使用 |>，旧项目继续使用 %>%。

11.17.6 管道操作的进阶案例

管道操作符的真正威力在于处理复杂的多步骤数据转换任务。当你需要对数据执行一系列操作时，管道操作符可以让代码逻辑清晰、易于理解和维护。进阶案例通常涉及多个数据框的连接、复杂的分组统计、条件筛选、以及数据重塑等操作。掌握这些进阶技巧，可以让你用简洁的代码完成原本需要几十行甚至上百行代码才能实现的功能。在生态学研究中，数据往往来自多个来源（野外调查、实验室测定、文献数据），需要整合、清洗、转换后才能分析。管道操作符提供了一种”流水线式”的数据处理思路，让每一步操作都清晰可见，便于检查和调试。以下案例展示了如何用管道操作符处理一个典型的森林样地数据集，从原始数据到最终的统计结果和可视化。

library(dplyr)
library(tidyr)
library(ggplot2)

# === 案例：森林碳储量数据处理 ===

# 步骤1：创建模拟数据
set.seed(2024)
forest_data <- tibble(
  plot_id = rep(1:5, each = 10),
  species = rep(c("马尾松", "杉木"), 25),
  dbh = rnorm(50, mean = 20, sd = 5),
  height = rnorm(50, mean = 12, sd = 3)
) %>%
  mutate(
    dbh = round(pmax(dbh, 5), 1),
    height = round(pmax(height, 3), 1)
  )

# 步骤2：复杂的管道操作
carbon_analysis <- forest_data %>%
  # 1. 筛选有效数据
  filter(dbh > 10, height > 5) %>%
  # 2. 计算生物量和碳储量
  mutate(
    biomass = 0.05 * dbh^2 * height,
    carbon = biomass * 0.5,  # 假设碳含量为 50%
    size_class = case_when(
      dbh < 15 ~ "小径木",
      dbh < 25 ~ "中径木",
      TRUE ~ "大径木"
    )
  ) %>%
  # 3. 按样地和物种分组统计
  group_by(plot_id, species) %>%
  summarise(
    n_trees = n(),
    mean_dbh = mean(dbh),
    mean_height = mean(height),
    total_carbon = sum(carbon),
    .groups = "drop"
  ) %>%
  # 4. 计算每个样地的物种多样性
  group_by(plot_id) %>%
  mutate(
    species_richness = n_distinct(species),
    carbon_proportion = total_carbon / sum(total_carbon)
  ) %>%
  ungroup() %>%
  # 5. 添加碳储量等级
  mutate(
    carbon_class = case_when(
      total_carbon < 50 ~ "低",
      total_carbon < 100 ~ "中",
      TRUE ~ "高"
    )
  ) %>%
  # 6. 排序
  arrange(plot_id, desc(total_carbon))

# 查看结果
print(carbon_analysis)

# === 进阶案例2：数据重塑与可视化 ===
summary_wide <- carbon_analysis %>%
  select(plot_id, species, total_carbon) %>%
  pivot_wider(
    names_from = species,
    values_from = total_carbon,
    values_fill = 0
  ) %>%
  mutate(total = 马尾松 + 杉木)

print(summary_wide)

# === 进阶案例3：多表连接 ===
# 假设有土壤数据
soil_data <- tibble(
  plot_id = 1:5,
  soil_pH = c(6.5, 5.8, 7.2, 6.1, 5.5),
  soil_carbon = c(18.5, 15.2, 22.3, 17.8, 14.5)
)

# 整合森林和土壤数据
integrated_data <- carbon_analysis %>%
  group_by(plot_id) %>%
  summarise(
    total_tree_carbon = sum(total_carbon),
    mean_dbh = mean(mean_dbh),
    .groups = "drop"
  ) %>%
  left_join(soil_data, by = "plot_id") %>%
  mutate(
    total_carbon = total_tree_carbon + soil_carbon,
    tree_carbon_ratio = total_tree_carbon / total_carbon
  )

print(integrated_data)

# === 进阶案例4：条件分组与嵌套操作 ===
nested_analysis <- forest_data %>%
  mutate(biomass = 0.05 * dbh^2 * height) %>%
  group_by(plot_id, species) %>%
  nest() %>%
  mutate(
    model = map(data, ~lm(biomass ~ dbh, data = .x)),
    r_squared = map_dbl(model, ~summary(.x)$r.squared)
  ) %>%
  select(plot_id, species, r_squared) %>%
  arrange(desc(r_squared))

print(nested_analysis)

运行结果说明：第一个案例展示了一个完整的数据处理流程，从筛选、计算、分组统计到分类和排序，所有操作都在一个管道中完成。pivot_wider() 将长格式数据转换为宽格式，便于比较不同物种的碳储量。left_join() 实现了森林数据和土壤数据的整合。嵌套分析使用 nest() 和 map() 对每个分组拟合线性模型，展示了管道操作符与函数式编程的结合。

生态学应用案例：在马尾松人工林碳储量研究中，这种复杂的管道操作非常常见。例如，你需要：1) 整合树木调查数据、土壤碳数据、凋落物碳数据；2) 按坡位、林龄、经营措施分组统计；3) 计算各碳库的比例和总碳储量；4) 分析碳储量与环境因子的关系。使用管道操作符，整个分析流程可以写成一个清晰的”数据流”，每一步的输入和输出都一目了然，便于检查和修改。这种编程风格不仅提高了代码的可读性，也大大降低了出错的概率。

11.18 正则表达式基础

处理复杂的文本模式：

# 示例数据
sample_ids <- c("Plot_A_01", "Plot_B_02", "Plot_A_03",
                "Plot_C_01", "Plot_B_04")

# 提取样地编号
str_extract(sample_ids, "[A-C]")

# 提取数字
str_extract(sample_ids, "\\d+")

# 替换模式
str_replace(sample_ids, "Plot_", "样地")

# 检测模式
str_detect(sample_ids, "A")

11.19 实战案例：森林样地数据分析

综合运用前面学到的知识：

# 创建模拟数据
set.seed(123)
forest_data <- data.frame(
  plot_id = rep(1:10, each = 5),
  species = rep(c("马尾松", "杉木", "桉树", "樟树", "楠木"), 10),
  height = rnorm(50, mean = 12, sd = 3),
  dbh = rnorm(50, mean = 20, sd = 5),
  crown_width = rnorm(50, mean = 4, sd = 1)
) %>%
  mutate(
    height = round(pmax(height, 5), 1),  # 确保高度 >= 5
    dbh = round(pmax(dbh, 10), 1),       # 确保胸径 >= 10
    crown_width = round(pmax(crown_width, 2), 1)
  )

# 数据探索
head(forest_data)
summary(forest_data)

# 按物种统计
species_stats <- forest_data %>%
  group_by(species) %>%
  summarise(
    n = n(),
    mean_height = mean(height),
    sd_height = sd(height),
    mean_dbh = mean(dbh),
    sd_dbh = sd(dbh)
  ) %>%
  arrange(desc(mean_height))

species_stats

# 计算生物量（简化公式）
forest_data <- forest_data %>%
  mutate(
    biomass = 0.05 * dbh^2 * height  # 简化的生物量估算公式
  )

# 按样地汇总
plot_summary <- forest_data %>%
  group_by(plot_id) %>%
  summarise(
    total_biomass = sum(biomass),
    mean_height = mean(height),
    species_richness = n_distinct(species)
  )

head(plot_summary)

# 可视化
par(mfrow = c(2, 2))  # 2x2 图形布局

# 1. 物种高度分布
boxplot(height ~ species, data = forest_data,
        main = "不同物种的高度分布",
        xlab = "物种", ylab = "高度 (m)",
        las = 2, col = "lightblue")

# 2. 高度与胸径关系
plot(forest_data$dbh, forest_data$height,
     col = as.numeric(factor(forest_data$species)),
     pch = 16,
     xlab = "胸径 (cm)", ylab = "高度 (m)",
     main = "高度-胸径关系")
legend("topleft", levels(factor(forest_data$species)),
       col = 1:5, pch = 16, cex = 0.7)

# 3. 样地生物量分布
barplot(plot_summary$total_biomass,
        names.arg = plot_summary$plot_id,
        main = "各样地总生物量",
        xlab = "样地编号", ylab = "生物量",
        col = "darkgreen")

# 4. 生物量直方图
hist(forest_data$biomass,
     main = "单株生物量分布",
     xlab = "生物量", col = "orange",
     breaks = 20)

par(mfrow = c(1, 1))  # 恢复默认布局

11.20 常见错误及解决方案

在学习 R 的过程中，以下是最常遇到的错误及其解决方法：

11.20.1 错误 1：对象未找到（Object not found）

# ❌ 错误示例
print(my_variable)  # Error: object 'my_variable' not found

# ✅ 解决方案
# 1. 检查变量名拼写
my_variable <- 10
print(my_variable)

# 2. 检查变量是否已创建
ls()  # 列出当前环境中的所有对象

# 3. 检查是否在正确的环境中
search()  # 查看搜索路径

原因： - 变量名拼写错误（R 区分大小写） - 变量未赋值 - 变量在不同的环境中

11.20.2 错误 2：下标越界（Subscript out of bounds）

下标越界错误是指试图访问向量、矩阵或数据框中不存在的元素。在 R 中，这类错误的表现形式比较特殊：对于向量，访问超出长度的索引会返回 NA 而不是报错；但对于数据框，访问不存在的列会直接报错”undefined columns selected”。理解这种差异对于调试代码非常重要。下标越界错误通常发生在以下场景：循环中使用了错误的索引范围、数据框列名拼写错误、数据维度与预期不符、动态索引计算错误等。在生态学数据分析中，这类错误尤其常见，因为数据集的大小和结构经常变化（如不同批次的调查数据、不同处理的实验数据），如果代码中硬编码了索引或列名，就容易出现越界问题。掌握如何预防和处理下标越界错误，是编写健壮代码的重要技能。

# === 向量索引越界 ===
x <- c(1, 2, 3)

# 情况1：访问超出长度的索引（返回 NA，不报错）
result1 <- x[5]
print(result1)  # 输出: NA

# 情况2：访问索引 0（返回空向量）
result2 <- x[0]
print(result2)  # 输出: numeric(0)

# 情况3：负索引（排除元素）
result3 <- x[-1]
print(result3)  # 输出: 2 3

# === 数据框索引越界 ===
df <- data.frame(a = 1:3, b = 4:6)

# 情况4：访问不存在的列（报错）
# df[, 5]  # Error: undefined columns selected

# 情况5：访问不存在的行（返回 NA）
result5 <- df[5, ]
print(result5)  # 输出: NA NA

# === 预防措施 ===

# 方法1：检查对象长度/维度
cat("向量长度:", length(x), "\n")
cat("数据框维度:", nrow(df), "行,", ncol(df), "列\n")

# 方法2：使用安全的索引方式
safe_index <- function(vec, idx) {
  if (idx > 0 && idx <= length(vec)) {
    return(vec[idx])
  } else {
    warning("索引超出范围")
    return(NA)
  }
}

result <- safe_index(x, 5)
print(result)

# 方法3：在循环前检查长度
for (i in 1:length(x)) {  # 正确：使用 length(x)
  print(x[i])
}

# 方法4：使用 dplyr 的安全函数
library(dplyr)
df %>% select(a, b)  # 明确指定列名，避免索引错误

运行结果说明：x[5] 返回 NA 而不是报错，这是 R 的特性，但可能导致隐蔽的错误。x[0] 返回空向量 numeric(0)，这在某些情况下会导致意外结果。df[5, ] 会返回一行全是 NA 的数据框。df[, 5] 会直接报错，因为数据框只有 2 列。使用 length() 和 dim() 可以在访问前检查对象大小，避免越界。

生态学应用案例：在马尾松林土壤数据分析中，下标越界错误的典型场景包括：

# 场景1：批量处理样地数据时索引错误
set.seed(2024)
soil_data <- data.frame(
  plot_id = 1:20,
  pH = rnorm(20, mean = 6.5, sd = 0.8),
  SOC = rnorm(20, mean = 18, sd = 5)
)

# ❌ 错误：硬编码循环次数
# for (i in 1:25) {  # 只有 20 个样地，但循环 25 次
#   cat("样地", i, "pH =", soil_data$pH[i], "\n")
# }

# ✅ 正确：使用 nrow() 动态获取行数
for (i in 1:nrow(soil_data)) {
  cat("样地", i, "pH =", soil_data$pH[i], "\n")
}

# 场景2：访问不存在的列
# ❌ 错误：列名拼写错误
# mean_TN <- mean(soil_data$TN)  # Error: 列 TN 不存在

# ✅ 正确：先检查列是否存在
if ("TN" %in% names(soil_data)) {
  mean_TN <- mean(soil_data$TN)
} else {
  cat("警告: 列 TN 不存在\n")
  mean_TN <- NA
}

# 场景3：子集化后索引失效
alkaline_plots <- soil_data[soil_data$pH > 7, ]
cat("碱性样地数量:", nrow(alkaline_plots), "\n")

# ✅ 正确：使用实际行数
if (nrow(alkaline_plots) > 0) {
  for (i in 1:nrow(alkaline_plots)) {
    print(alkaline_plots[i, ])
  }
}

这些例子展示了如何在实际分析中避免下标越界错误。关键原则是：1) 永远不要硬编码索引或长度；2) 使用 length()、nrow()、ncol() 动态获取维度；3) 在访问列之前检查列名是否存在；4) 使用 dplyr 等现代工具，用列名而不是索引访问数据。

11.20.3 错误 3：参数长度不匹配

# ❌ 错误示例
x <- c(1, 2, 3)
y <- c(10, 20)
# x + y  # Warning: longer object length is not a multiple of shorter object length

# ✅ 解决方案
# 1. 确保向量长度一致
x <- c(1, 2, 3)
y <- c(10, 20, 30)
x + y

# 2. 使用循环补齐（R 会自动循环，但要小心）
x <- c(1, 2, 3, 4)
y <- c(10, 20)  # 会循环为 c(10, 20, 10, 20)
x + y  # 结果：11, 22, 13, 24

# 3. 明确处理长度不一致的情况
if (length(x) != length(y)) {
  warning("向量长度不一致")
  # 截断到较短的长度
  min_len <- min(length(x), length(y))
  x <- x[1:min_len]
  y <- y[1:min_len]
}

11.20.4 错误 4：因子与字符混淆

# ❌ 错误示例
df <- data.frame(
  species = c("马尾松", "杉木", "桉树"),
  stringsAsFactors = TRUE  # 旧版 R 默认行为
)
class(df$species)  # "factor"
# df$species[1] <- "樟树"  # 可能产生 NA（如果 "樟树" 不在原有水平中）

# ✅ 解决方案
# 1. 使用 stringsAsFactors = FALSE（R 4.0+ 默认）
df <- data.frame(
  species = c("马尾松", "杉木", "桉树"),
  stringsAsFactors = FALSE
)
class(df$species)  # "character"
df$species[1] <- "樟树"  # 正常工作

# 2. 或者使用 tibble（推荐）
library(tibble)
df <- tibble(
  species = c("马尾松", "杉木", "桉树")
)
class(df$species)  # "character"

# 3. 转换因子为字符
df$species <- as.character(df$species)

11.20.5 错误 5：包未安装或未加载

# ❌ 错误示例
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) + geom_point()
# Error: could not find function "ggplot"

# ✅ 解决方案
# 1. 检查包是否已安装
if (!require(ggplot2)) {
  install.packages("ggplot2")
  library(ggplot2)
}

# 2. 加载包
library(ggplot2)
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) + geom_point()

# 3. 使用 :: 直接调用（不需要加载包）
ggplot2::ggplot(iris, ggplot2::aes(x = Sepal.Length, y = Sepal.Width)) + 
  ggplot2::geom_point()

11.20.6 错误 6：数据类型不匹配

# ❌ 错误示例
x <- "10"
y <- 5
# x + y  # Error: non-numeric argument to binary operator

# ✅ 解决方案
# 1. 显式类型转换
x <- as.numeric("10")
y <- 5
x + y

# 2. 检查数据类型
class(x)
is.numeric(x)

# 3. 批量转换数据框中的列
df <- data.frame(
  a = c("1", "2", "3"),
  b = c("4", "5", "6")
)
df[] <- lapply(df, as.numeric)  # 转换所有列为数值型
str(df)

11.20.7 错误 7：缺失值处理不当

# ❌ 错误示例
x <- c(1, 2, NA, 4, 5)
mean(x)  # 返回 NA

# ✅ 解决方案
# 1. 使用 na.rm = TRUE
mean(x, na.rm = TRUE)

# 2. 提前移除 NA
x_clean <- x[!is.na(x)]
mean(x_clean)

# 3. 使用 na.omit()
x_clean <- na.omit(x)
mean(x_clean)

# 4. 检查 NA 的位置
which(is.na(x))

# 5. 统计 NA 的数量
sum(is.na(x))

11.20.8 错误 8：文件路径问题

# ❌ 错误示例（Windows 路径）
# data <- read.csv("C:\Users\Data\file.csv")  # 错误：\ 是转义字符

# ✅ 解决方案
# 1. 使用正斜杠
data <- read.csv("C:/Users/Data/file.csv")

# 2. 使用双反斜杠
data <- read.csv("C:\\Users\\Data\\file.csv")

# 3. 使用相对路径（推荐）
data <- read.csv("data/file.csv")

# 4. 检查当前工作目录
getwd()

# 5. 设置工作目录
setwd("C:/Users/YourName/Project")

# 6. 使用 here 包（最推荐）
library(here)
data <- read.csv(here("data", "file.csv"))

11.20.9 生态学实例：常见数据处理错误

# 创建包含常见问题的数据
set.seed(2024)
problematic_data <- data.frame(
  plot_id = c(1, 2, 3, 4, 5),
  species = c("马尾松", "杉木", NA, "桉树", "樟树"),
  dbh = c("20.5", "15.8", "25.3", "invalid", "18.7"),  # 字符型
  height = c(12.5, NA, 15.7, 10.2, 9.8),
  count = c(10, 5, 3, -2, 8)  # 包含负值
)

# 问题诊断与修复
diagnose_and_fix <- function(data) {
  cat("=== 数据诊断 ===\n")
  
  # 1. 检查数据结构
  cat("\n1. 数据结构：\n")
  str(data)
  
  # 2. 检查缺失值
  cat("\n2. 缺失值统计：\n")
  print(colSums(is.na(data)))
  
  # 3. 检查数据类型
  cat("\n3. 数据类型问题：\n")
  if (!is.numeric(data$dbh)) {
    cat("- dbh 列不是数值型\n")
  }
  
  # 4. 修复数据
  cat("\n=== 数据修复 ===\n")
  
  # 转换 dbh 为数值型（无效值变为 NA）
  data$dbh <- suppressWarnings(as.numeric(data$dbh))
  cat("- 已转换 dbh 为数值型\n")
  
  # 移除包含 NA 的行
  data_clean <- na.omit(data)
  cat("- 已移除包含 NA 的行：", nrow(data) - nrow(data_clean), "行\n")
  
  # 移除负值
  data_clean <- data_clean[data_clean$count >= 0, ]
  cat("- 已移除负值\n")
  
  cat("\n=== 修复后数据 ===\n")
  print(data_clean)
  
  return(data_clean)
}

clean_data <- diagnose_and_fix(problematic_data)

11.21 课后练习

11.21.1 基础练习（必做）

练习 1：创建和操作数据框

创建一个包含 5 种植物的数据框，包含物种名、株高、叶面积、是否为乔木四列，并进行基本统计分析。

# 你的代码：
# 1. 创建数据框
plant_data <- data.frame(
  species = c("___", "___", "___", "___", "___"),
  height = c(___, ___, ___, ___, ___),
  leaf_area = c(___, ___, ___, ___, ___),
  is_tree = c(___, ___, ___, ___, ___)
)

# 2. 查看数据结构
str(plant_data)

# 3. 计算株高的统计量
mean_height <- ___
sd_height <- ___
cv_height <- ___

cat("平均株高:", mean_height, "\n")
cat("标准差:", sd_height, "\n")
cat("变异系数:", cv_height, "%\n")

练习 2：数据筛选

筛选出株高大于平均值的植物。

# 你的代码：
tall_plants <- plant_data[___, ]
print(tall_plants)

练习 3：数据可视化

用 hist() 绘制 iris 数据集中 Sepal.Width 的分布直方图。

# 你的代码：
hist(___, 
     main = "___",
     xlab = "___",
     col = "___")

练习 4：编写函数

编写一个函数计算 Shannon 多样性指数，并用 c(10, 5, 3, 2, 1) 测试。

# 你的代码：
shannon_index <- function(abundance) {
  # 计算相对丰度
  p <- ___
  
  # 计算 Shannon 指数: H' = -Σ(pi * ln(pi))
  H <- ___
  
  return(H)
}

# 测试
test_data <- c(10, 5, 3, 2, 1)
result <- shannon_index(test_data)
cat("Shannon 多样性指数:", result, "\n")

练习 5：管道操作符

使用管道操作符分析 iris 数据：筛选 Petal.Length > 4 的记录，按物种分组计算平均花瓣宽度。

library(dplyr)

# 你的代码：
result <- iris %>%
  filter(___) %>%
  group_by(___) %>%
  summarise(___)

print(result)

11.21.2 进阶练习（选做）

练习 6：综合数据处理

创建一个森林样地数据集，计算每个样地的生物量，并按物种分组统计。

library(dplyr)

# 你的代码：
set.seed(2024)
forest_data <- data.frame(
  plot_id = rep(1:5, each = 4),
  species = rep(c("马尾松", "杉木", "桉树", "樟树"), 5),
  dbh = rnorm(20, mean = 20, sd = 5),
  height = rnorm(20, mean = 12, sd = 3)
) %>%
  mutate(
    dbh = round(pmax(dbh, 5), 1),
    height = round(pmax(height, 3), 1),
    biomass = ___  # 使用公式: 0.05 * dbh^2 * height
  )

# 按物种分组统计
species_summary <- forest_data %>%
  group_by(___) %>%
  summarise(
    n_trees = ___,
    mean_biomass = ___,
    total_biomass = ___
  )

print(species_summary)

11.21.3 提交要求

1. 将你的代码保存为 r-basics-exercises.R 文件
2. 确保代码可以运行，没有语法错误
3. 添加必要的注释说明你的思路
4. 提交到 GitHub 仓库的 homework/ 目录
5. 在 README.md 中记录完成情况和遇到的问题

11.21.4 评分标准

项目	分值	说明
代码正确性	40%	代码能正确运行，结果准确
代码规范性	30%	变量命名清晰，代码格式规范
注释完整性	20%	关键步骤有注释说明
创新性	10%	使用了课程中学到的高级技巧

11.22 扩展阅读

• R for Data Science：数据科学经典教材
• Advanced R：深入理解 R 语言
• ggplot2 官方文档：数据可视化
• dplyr 官方文档：数据处理
• vegan 包教程：群落生态学分析
• R Graphics Cookbook：图形绘制指南