赖江山老师《数量生态学------R语言的应用》第二版阅读笔记1
本次是第一章绪论和第二章绘图部分
绪论
生态学（Ecology）是研究有机体及其周围环境相互关系的科学。
读赖江山翻译的《数量生态学------R语言的应用》第二版，英文版见Numerical Ecology with R（网址提供书中用到的数据以及演示代码）链接
数量生态学网站.png
本书用到的主要的数据集：
rm(list=ls())
#设置当前工作目录
setwd("E:/Rstudio_working/lai jian shan/DATA/Data")
#导入物种多度数据
spe<-read.csv("DoubsSpe.csv",row.names = 1)
View(spe) #27种鱼类在每个样方的多度。
#导入环境数据
env<-read.csv("DoubsEnv.csv",row.names = 1)
View(env) #包括11个与河流水文地形和水体化合物相关的环境变量
#导入空间坐标数据
spa<-read.csv("DoubsSpa.csv",row.names = 1)
View(spa)#样方的地理坐标（笛卡尔坐标系，x和y）
第二章 探索性数据分析(第一部分)
数据探索
数据提取
#设置工作目录
setwd("E:/Rstudio_working/lai jian shan/DATA/Data")
#加载包、数据和函数
library(vegan)
library(leaflet)
library(RgoogleMaps)
library(googleVis)
library(labdsv)
#加载后面要用到的函数
source("panelutils.R")
# 导入Doubs 数据
load("Doubs.RData")
# 这个文件包含如下这些对象
#spe：物种(群落)数据框(鱼类多度数据)
head(spe,c(5,10))
#env:环境因子数据框
head(env)
#spa：空间位置数据框-笛卡尔坐标系
head(spa,10)
#fishtraits:鱼类功能属性
head(fishtraits,c(5,10))
#latlong：空间位置数据框-经纬度
head(latlong)
物种数据：第一次接触
首先使用基础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 = "Distribution of abundance classes",
noRStudioGD = TRUE)
#所有物种混合在一起的多度分布柱状图
barplot(ab,
las=1,
xlab = "多度等级",
ylab="频度",
col=gray(5:0/5))
# 多度数据中0值的数量
sum(spe==0)
# 多度数据中0值所占比例
sum(spe==0)/(nrow(nrow(spe)*ncol(spe)))
多度分布柱状图.png
如何解释0值(缺失)在数据框内频率这么高？
答：其实造成缺失的因素有很多，但是有两种需要我们的注意：
- 真实的环境适合这个物种生存，只是我们采样的时候没采到（比如人家冬眠了，出去玩了，尚未迁徙到这里）。
- 真实的环境不适合这个物种生存，在这里生存就会被淘汰。 所以对于零值我们要小心处理，关键还是理解数据的生物学意义。
物种数据：进一步分析------生成样方位置分布图
#样方位置地图
#生成空的绘图底板(横纵坐标轴比例1:1(参数asp)，带标题)
dev.new(title = "Site Locations", width = 9, noRStudioGD = TRUE)
# 从spa数据框获取地理坐标X和Y
plot(spa,
asp=1,
type="n",
main="样方位置",
xlab="X轴坐标(km)",
ylab="Y轴坐标(km)")
#加一条连接各个样方点的蓝色线
lines(spa,col="light blue")
# 添加每个样方的编号
text(spa,row.names(spa),cex=0.8,col="red")
#添加文本
text(68, 20, "上游", cex = 1.2, col = "red")
text(15, 35, "下游", cex = 1.2, col = "red")
样方位置分布图.png
30个样方沿着Doubs河的空间分布。绘制这幅图用到的plot函数是R的基础绘图函数。可以?plot()查看其帮助文档，asp是用来调整绘图版面的长宽比列的。
1.当数据覆盖足够大的区域时，可以投影采样点到googleVis地图上
#这里由于浏览器问题所以这个步奏没办法运行
#将样方点映射到谷歌地图
#以googleVis 程序包默认的方法组织数据
#从浏览器中生成出图
nom <- latlong$Site
latlong2 <- paste(latlong$LatitudeN, latlong$LongitudeE, sep = ":")
df <- data.frame(latlong2, nom, stringsAsFactors = FALSE)
mymap1 <- gvisMap(df,
locationvar = "latlong2",
tipvar = "nom",
options = list(showTip = TRUE)
)
plot(mymap1)
2.NOT in the BOOOK-使用Teaflet 使用传单动态地图,投影到开放街道地图背景上的站点f传单使用查看器窗格I显示其输出
longitude <- latlong$LongitudeE
latitude <- latlong$LatitudeN
site <- as.character(latlong$Site)
background <- addTiles(leaflet())
mymap1 <-
addMarkers(
background,
lng = longitude,
lat = latitude,
label = site,
labelOptions = labelOptions(noHide = TRUE, textOnly = TRUE)
)
mymap1
效果图.png
把物种数据映射到采样点之上，可以看到物种是怎样随着河流变化的。
#某些鱼类的分布地图
#将绘图窗口分割为4个区域，两行两列
#ex.axis表示修改坐标轴刻度字体大小，cex.lab表示修改坐标轴名称字体大小，cex.main表示修改标题字体大小
par(mfrow=c(2,2))
plot(spa,
asp=1,
cex.axis=0.8,
col="brown",
cex=spe$Satr,
main="褐鳟",
xlab="x坐标(km)",
ylab="y坐标(km)")
lines(spa, col="light blue")
plot(spa,
asp=1,
cex.axis=0.5,
col="brown",
cex=spe$Thth,
main="茴鱼",
xlab="x坐标 (km)",
ylab="y坐标 (km)")
lines(spa, col="light blue")
plot(spa,
asp=1,
cex.axis=0.5,
col="brown",
cex=spe$Baba,
main="鲃鱼",
xlab="x坐标 (km)",
ylab="y坐标 (km)")
lines(spa, col="light blue")
plot(spa,
asp=1,
cex.axis=0.5,
col="brown",
cex=spe$Abbr,
main="欧鳊",
xlab="x坐标 (km)",
ylab="y坐标 (km)")
lines(spa, col="light blue")
某些鱼类的分布地图.png
从图上我们清楚地看到，褐鳟、茴鱼、鲃鱼、欧鳊的多度是沿着Doubs河从上游到下游分布的，也就明白为什么Verneaux 选择这4种鱼类作为不同区域的生态指示种。注意之前提到的零值问题，这里是同一条河流不会有迁移的障碍，这几种鱼的生活史也较接近不存在有冬眠的不一致的情况。
注意:plot()函数的参数cex的用法，它的作用是定义数据点标识的大小。
每个物种在多少样方中出现？可以看看物种的相对频度。
#比较物种频度
#计算每个物种出现的样方数
#按列进行计数，因此apply()第二个参数MARGIN应该设定为2
#spe>0表示先将spe内的数值转化为逻辑向量T/F，然后对逻辑值T进行列的汇总
spe.pres <- apply(spe>0,2,sum)
#按照升序的方式重新排列结果
sort(spe.pres)
#计算频度百分比
spe.relf <- 100*spe.pres/nrow(spe)
spe.pres
#设置排列结果为1位小数
round(sort(spe.relf),1)
#绘制柱状图
#将绘图窗口垂直一分为二
par(mfrow = c(1,2))
hist(spe.pres,
main = "物种出现数",
right = FALSE,
las = 1,
xlab = "出现数",
ylab = "物种数量",
breaks = seq(0, 30, by = 5),
col = "bisque"
)
hist(spe.relf,
main = "Species Relative Frequencies",
right = FALSE,
las = 1,
xlab = "出现率 (%)",
ylab = "物种数量",
breaks = seq(0, 100, by = 10),
col = "bisque"
)
柱状图.png
注意：每个样方内存在多少物种（物种的丰度）。思考频度与丰度的差异。
#样方比较：物种丰富度
#计算每个样方内物种数
#以行汇总，apply()函数第二个参数MARGIN应该设定为1
sit.pres <- apply(spe >0, 1, sum)
#按照升序的方式重新排列结果
sort(sit.pres)
#将绘图窗口垂直一分为二
par(mfrow=c(1,2))
#绘制样方沿着河流的分布位置和所含物种丰富度
#type="s",表示绘制数值直接的阶梯图
plot(sit.pres,type = "s",las=1,col="gray",
main = "物种丰富度-上下游的梯度",
xlab="样方沿着河流的编号",ylab="物种丰富度")
text(sit.pres,row.names(spe),cex = 0.8,col = "red")
#使用地理坐标绘制气泡地图
plot(spa,asp=1,main="物种丰富度地图",pch=21,col="white",
bg="brown",cex=5*sit.pres/max(sit.pres),
xlab="x坐标(km)",ylab="y坐标(km)")
lines(spa,col="light blue")
物种丰富度-上下游的梯度.png
你能否辨析沿着河流哪里是物种丰富度的热点地区？
用vegan包中的diversity()函数计算生物多样性指数
#计算生物多样性指数
# *****************
# 载入所需要的vegan程序包
library(vegan) # 如果未载入，需要执行这一步
#访问diversity() 帮助界面
?diversity
N0 <- rowSums(spe > 0) #物种丰富度
H <- diversity(spe) # Shannon熵指数
N1 <- exp(H) # Shannon 多样性指数
N2 <- diversity(spe, "inv") # Simpson多样性指数
J <- H/log(N0) # Pielou 均匀度
E1 <- N1/N0 # Shannon均匀度 (Hill比率)
E2 <- N2/N0 # Simpson均匀度 (Hill比率)
div <- data.frame(N0, H, N1, N2, E1, E2, J)
div
环境数据
使用summary()函数了解这些环境数据数值分布和空间分布与物种数据有哪些不同。
绘制部分环境数据的地图，首先是气泡图
#部分环境变量的气泡地图
#cex参数定义气泡的大小
par(mfrow=c(2,2))
plot(spa,asp=1,main="海拔",pch=21,col="white",bg="red",
cex=5*env$ele/max(env$ele),xlab="x",ylab="y")
lines(spa,col="light blue")
plot(spa,asp=1,main="流量",pch=21,col="white",bg="blue",
cex=5*env$dis/max(env$dis),xlab="x",ylab="y")
lines(spa,col="light blue")
plot(spa,asp=1,main="氧含量",pch=21,col="white",bg="green3",
cex=5*env$oxy/max(env$oxy),xlab="x",ylab="y")
lines(spa,col="light blue")
plot(spa,asp=1,main="硝酸盐浓度",pch=21,col="white",bg="brown",
cex=5*env$nit/max(env$nit),xlab="x",ylab="y")
lines(spa,col="light blue")
气泡图.png
哪幅图最能展示上下游的梯度？如何解释其他环境变量的空间分布格局？
海拔反映了环境变量的梯度。流量和海拔都好解释，想一下含氧量和硝酸盐浓度只靠这几张图能解释吗？在（150,200）处的硝酸盐浓度很高而氧含量很低，为什么？
线条图
par(mfrow=c(2,2))
plot(env$dfs, env$ele,
type = "l",
xlab = "离源头距离 (km)",
ylab = "海拔 (m)",
col = "red", main = "海拔"
)
plot(env$dfs, env$dis,
type = "l",
xlab = "离源头距离 (km)",
ylab = "流量 (m3/s)",
col = "blue",
main = "流量"
)
plot(env$dfs, env$oxy,
type = "l",
xlab = "离源头距离 (km)",
ylab = "氧含量 (mg/L)",
col = "green3",
main = "氧含量"
)
plot(env$dfs, env$nit,
type = "l",
xlab = "离源头距离 (km)",
ylab = "硝酸盐浓度 (mg/L)",
col = "brown",
main = "硝酸盐浓度"
)
线条图.png
如果需要了解任意两个环境变量之间的关系，可以使用功能强大的矩阵散点图绘图函数pairs()
另外，也可以使用panelutils.R函数为每一个双变量散点图添加LOWESS平滑线，并且=绘制每个变量频度分布表
#所有变量对之间的二维散点图
#带频度分布的柱状图和光滑拟和曲线的双变量散点图
pairs(env,panel = panel.smooth,diag.panel = panel.hist,
main="双变量散点图(带频度分布图和平滑曲线)")
注意:每个散点图展示对角线上两个变量之间的关系。散点图的横坐标对应上方或下方的变量，纵坐标对应左侧或右侧的变量
好了绪论和第二章的一部分就到这里，下一次给大家带来第二章的数据探索，这次的部分画图数据处理是经过数据处理的，详细请看第一次推送！！！
感谢你的阅读，有问题可以联系小蘑菇。
