数量生态学R语言应用——探索性数据分析

《赖江山老师的数量生态学——R语言应用第二版》学习笔记，理论部分主要是赖江山老师的，代码有些我自己进行了修改优化，使用tidyverse和国内的在线地图服务改写。

• 假设检验(hypothesis testing)和建模(modelling)已经成为当今生态学研究的主要手段。然而，利用可视化工具或计算综合统计量对原始数据进行探索性数据分析(Exploratory Data Analysis，EDA)仍然必不可少，也是多维数据统计分析的开始步骤，其目的在于:

• 了解数据概况;
• 对某些变量进行转化或重新编码:
• 确定下一步分析的方向。

• 本章以案例数据为例，演示R标准程序包内与探索性数据分析相关的函数，主要内容有:

• 学习或修改R的一些基本函数;
• 学习可以应用于生态学多维数据的探索性分析技术;·以Doubs数据集内的水文数据为例，演示探索性数据分析。

数据代码下载

网上公开资源：

• 网址：http://www.numericalecology.com/numecolR/

数据集介绍

Doubs数据集变量介绍：

探索Doubs数据集基本情况

当然，以下是R语言代码的中文翻译注释：

• 加载数据。假设Doubs.RData文件在
• Doubs.RData文件包含以下对象：
• spe: 物种（群落）数据框（鱼类多度）
• env: 环境数据框
• spa: 空间数据框 - 笛卡尔坐标
• fishtraits: 鱼类物种的功能特征
• latlong: 空间数据框 - 纬度和经度

代码和解释

library(tidyverse)

Doubs = load("./NEwR-2ed_code_data/NEwR2-Data/Doubs.RData")

# 使用基本R函数探索数据框 =============
spe                       # 在控制台中显示整个数据框，对于大型数据集不推荐
spe[1:5, 1:10]            # 仅显示5行和10列
head(spe)                 # 仅显示前6行
tail(spe)                 # 仅显示最后6行
nrow(spe)                 # 行数（站点）
ncol(spe)                 # 列数（物种）
dim(spe)                  # 数据框的维度（行，列）
colnames(spe)             # 列标签（描述符=物种）
rownames(spe)             # 行标签（对象=站点）
summary(spe)              # 列的描述性统计

# 整个数据集中多度值的最小值和最大值
range(spe)
# 每个物种的最小值和最大值
apply(spe, 2, range)
# 统计每个多度类别的案例数
(ab <- table(unlist(spe)))
# 创建一个带有标题的图形窗口
dev.new(title = "多度类别的分布", 
  noRStudioGD = TRUE
)

# 多度等级数据条形图
p1 = ggplot(ab, aes(x = AbundanceClass, y = Frequency)) +
  geom_bar(stat = "identity", fill = gray(0.5)) +
  labs(title = "多度等级的分布", x = "多度等级", y = "频率") +
  theme_minimal()
p1 
# ggsave(filename = "./Rplots/p1.jpg", p1, height = 8, width = 12, units = "cm", dpi = 300)

# 缺失值的数量
sum(spe == 0)
# 群落数据集中零的比例
sum(spe == 0) / (nrow(spe) * ncol(spe))

#绘制样方分布
p2 = ggplot(spa, aes(x = X, y = Y)) +
  geom_point() +  # 添加点
  geom_path() +  # 添加连接线,一定得是geom_path，按顺序绘制连线
  scale_x_continuous(limits = c(0, 160), breaks = seq(0, 160, 30)) +  # 设置X轴范围和间隔
  scale_y_continuous(limits = c(0, 120), breaks = seq(0, 120, 30)) +  # 设置Y轴范围和间隔
  labs(title = "站点位置", x = "x坐标（公里）", y = "y坐标（公里）") +
  theme_minimal() +  # 使用简洁主题
  theme(plot.title = element_text(hjust = 0.5)) +  # 居中标题
  geom_text(aes(label = row.names(spa)), check_overlap = TRUE, vjust = -1) +  # 添加站点标签，避免重叠
  annotate("text", x = min(spa$X), y = min(spa$Y), label = "下游", size = 5, color = "red", hjust = 0) +  # 起点标注“下游”
  annotate("text", x = max(spa$X), y = max(spa$Y), label = "上游", size = 5, color = "red", hjust = 1)   # 终点标注“上游”
p2

apply()函数详解

apply() 函数是 R 语言中的一个非常有用的函数，它允许你对数组或矩阵的边缘（即行或列）应用一个函数，并返回结果。这个函数特别适用于当你需要对数据的子集执行操作时，比如计算行或列的统计量。下面是 apply() 函数的基本用法和参数解释：

函数描述

apply(X, MARGIN, FUN, …)

• X：一个数组，包括矩阵。
• MARGIN：一个向量，指定函数将被应用的边缘。例如，对于矩阵，1 表示行，2 表示列，c(1, 2) 表示行和列。如果 X 有命名的维度名称，可以是选择维度名称的字符向量。
• FUN：要应用的函数。如果函数是如 +, %*% 等，函数名必须使用反引号或引号括起来。
• …：传递给 FUN 的可选参数。

返回值

• 如果 FUN 的每次调用返回长度为 n 的向量，那么 apply 返回一个维度为 c(n, dim(X)[MARGIN]) 的数组，如果 n > 1。如果 n 等于 1，apply 返回一个向量（如果 MARGIN 长度为 1）或者维度为 dim(X)[MARGIN] 的数组（否则）。
• 如果 n 是 0，结果长度为 0，但不一定是“正确”的维度。
• 如果对 FUN 的调用返回不同长度的向量，apply 返回一个长度为 prod(dim(X)[MARGIN]) 的列表，如果 MARGIN 的长度大于一，则 dim 设置为 MARGIN。

在所有情况下，结果在设置维度之前会被 as.vector 强制转换为基本向量类型之一，因此（例如）因子结果将被强制转换为字符数组。

详细说明

• 如果 X 不是数组而是具有非空 dim 值的对象（例如数据框），apply 会尝试通过 as.matrix（如果是二维的，例如数据框）或 as.array 将其强制转换为数组。
• FUN 是通过 match.fun 调用找到的，通常是函数、符号（例如，反引号名称）或指定要搜索的函数的字符字符串。
• … 中的参数不能与其他参数同名，需要小心避免与 MARGIN 或 FUN 部分匹配。在通用代码中，如果 … 被传递，最好命名前三个参数：这既避免了与 MARGIN 或 FUN 的部分匹配，也确保如果传递了名为 X、MARGIN 或 FUN 的参数，会给出合理的错误消息。

示例

在上面的例子中，使用apply()函数计算了每个物种的多度最大值最小值。

# 每个物种的最小值和最大值
apply(spe, 2, range)

先看一下物种数据的组织情况，每个物种一个字段，按列分布，因此，要想计算每个物种多度的最大最小值，需要按列计算，因此apply()函数中间的参数为2。

在交互式地图上显示Doubs鱼类采样点

• Doubs鱼类采样点位置信息存储在latlong数据框中，如下图所示，包含三个字段：
• Site站点号
• LatitudeN 北纬
• LongitudeE 东经

练习代码使用的是谷歌地图服务，但是国内不能使用，因此我改成了中科星图星图地球数据云图源，图源介绍可以查阅：QGIS又多了一个在线底图新选择，星图地球数据云测评，超好用！真香！

#使用中科星图数据云
library(leaflet)
geovis_token = "6f6a548b30063f2aeeb149a5589f1362076fe636a113c690fd3864b14fe01d3e"  #从星图地球数据云网站获取
#调用影像地图和标注的服务
geovis_satellite_url = paste0("https://tiles1.geovisearth.com/base/v1/img/{z}/{x}/{y}?format=webp&tmsIds=w&token=", geovis_token)
geovis_label_url = paste0("https://tiles1.geovisearth.com/base/v1/cia/{z}/{x}/{y}?format=webp&tmsIds=w&token=", geovis_token)
geovismap = leaflet() %>% 
  addTiles(urlTemplate = geovis_satellite_url) %>% 
  addTiles(urlTemplate = geovis_label_url)
#把采样点标注到交互式地图上
longitude <- latlong$LongitudeE
latitude <- latlong$LatitudeN
site <- as.character(latlong$Site)
mymap1 <-
  addMarkers(
    geovismap,
    lng = longitude,
    lat = latitude,
    label = site,
    labelOptions = labelOptions(noHide = TRUE, textOnly = FALSE)
  )
mymap1

鱼类多度沿Doubs河的分布

前面已经绘制了河流和采样点，鱼在哪里呢？接下来以Doubs河流生态区4种指示鱼类为例，展示鱼类多度沿Doubs河的分布情况。

使用spa数据框中的空间位置，来绘制位置图，使用spe数据框中的物种多度来表示鱼类物种沿河变化的情况。

#使用spa数据框来描述鱼类多度沿河变化的情况
#每种鱼类绘制一个图，使用patchwork组合起来
p1 = ggplot(spa, aes(x = X, y = Y, color = spe$Satr)) +
  geom_point(aes(size = spe$Satr), alpha = 0.6) +
  geom_path() +
  labs(title = "Brown trout", x = "x coordinate (km)", y = "y coordinate (km)") +
  theme(legend.position = "none")

p2 = ggplot(spa, aes(x = X, y = Y, color = spe$Thth)) +
  geom_point(aes(size = spe$Thth), alpha = 0.6) +
  geom_path() +
  labs(title = "Grayling", x = "x coordinate (km)", y = "y coordinate (km)") +
  theme(legend.position = "none")

p3 = ggplot(spa, aes(x = X, y = Y, color = spe$Baba)) +
  geom_point(aes(size = spe$Baba), alpha = 0.6) +
  geom_path() +
  labs(title = "Barbel", x = "x coordinate (km)", y = "y coordinate (km)") +
  theme(legend.position = "none")

p4 = ggplot(spa, aes(x = X, y = Y, color = spe$Abbr)) +
  geom_point(aes(size = spe$Abbr), alpha = 0.6) +
  geom_path() +
  labs(title = "Common bream", x = "x coordinate (km)", y = "y coordinate (km)") +
  theme(legend.position = "none")

# 使用patchwork组合图层
final_plot <- (p1 + p2 + p3 + p4) + plot_layout(ncol = 2, nrow = 2)
final_plot

每个物种在多少个样方中出现？可以计算物种相对频度并绘制直方图。

# 计算每个物种出现的站点数
spe_pres <- spe %>% 
  summarise(across(everything(), ~sum(. > 0))) %>% 
  gather(key = "species", value = "occurrences") %>% 
  arrange(occurrences)

# 排序并计算百分比频率
spe_relf <- spe_pres %>% 
  mutate(relative_frequency = 100 * occurrences / nrow(spe))

# 输出结果四舍五入
round(spe_relf %>% arrange(occurrences) %>% pull(relative_frequency), 1)

# 绘制直方图
# 物种出现次数的直方图
p5 = ggplot(spe_pres, aes(x = occurrences)) +
  geom_histogram(binwidth = 5, color = "black", fill = "bisque") +
  labs(title = "物种出现次数", x = "出现次数", y = "物种数量") +
  theme_minimal()

# 物种相对频率的直方图
p6 = ggplot(spe_relf, aes(x = relative_frequency)) +
  geom_histogram(binwidth = 10, color = "black", fill = "bisque") +
  labs(title = "物种相对频率", x = "出现频率 (%)", y = "物种数量") +
  theme_minimal()
p56 = p5+p6
p56

接下来了解每个样方内有多少个物种存在

#计算每个样方内有多少个物种存在
# 计算每个站点的物种丰富度
sit_pres <- spe %>%
  rowSums(. > 0) %>%
  enframe(name = "site", value = "species_richness") %>% 
  mutate(site = as.integer(site))

# 绘制物种丰富度与站点位置（上游-下游梯度）的关系图
p7 = ggplot(sit_pres, aes(x = site, y = species_richness)) +
  geom_step(color = "gray") +   #绘制阶梯线图
  geom_point(color = "gray") +
  labs(title = "物种丰富度与 \n 上游-下游梯度的关系", 
       x = "站点编号", 
       y = "物种丰富度") +
  theme_minimal() +
  geom_text(aes(label = site), size = 3, color = "red", hjust = 1.5)

# 使用地理坐标绘制气泡图
p8 = spa %>%
  bind_cols(sit_pres) %>%
  ggplot(aes(x = X, y = Y, size = species_richness)) +
  geom_point(color = "brown", fill = "white", shape = 21) +
  scale_size_continuous(range = c(1, 10)) +
  labs(title = "物种丰富度地图", 
       x = "x 坐标 (km)", 
       y = "y 坐标 (km)") +
  theme_minimal() +
  geom_path(color = "light blue", linewidth = 0.8)
p78 = p7+p8
p78