上节课我们简单了解了《数量生态学》的基本内容,了解了书中用到的数据集。这节课的索引图见手记。
就像SPSS有一个菜单叫描述统计。本章也属于这个范畴。
首先对物种群落数据有一个简单的统计描述,也是看看我们的数据格式是否正常。
rm(list=ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\Run")
#导入物种多度数据
spe<-read.csv("../DATA/DoubsSpe.csv",row.names = 1)
str(spe)
#导入环境数据
env<-read.csv("../DATA/DoubsEnv.csv",row.names = 1)
#导入空间坐标数据
spa<-read.csv("../DATA/DoubsSpa.csv",row.names = 1)
我们对物种群落数据做一个简单的描述性统计,同时也是看看我们的数据格式是否正常。
# 基础函数
# ********
spe #在控制台显示整个数据框的内容,但对于大样本的数据框
#并不建议直接显示
spe[1:5,1:10] #只展示前5行和前10列
head(spe) #只展示前几行
nrow(spe) #提取数据框总行数
ncol(spe) #提取数据框总列数
dim(spe) #提取数据框的维度(显示数据框多少行,多少列)
colnames(spe) #提取列名,在这里是物种名
rownames(spe) #提取行名,在这里一行代表一个样方
summary(spe) #以列为单位,对列变量进行描述性统计
#比较多度的中位值和平均值。大部分是对称分布吗?
如果多度分布是对称的,中位数应该和平均值差别不大。但是显然并不是对称的。
# 多度数据总体分布情况
# *******************
# 整个多度数据值的范围
range(spe)
# 计算每种多度值的数量
ab <- table(unlist(spe))
ab
# 所有种混和在一起的多度分布柱状图
barplot(ab, las=1, xlab="多度等级", ylab="频度", col=gray(5:0/5))
# 多度数据中0值的数量
sum(spe==0)
# 多度数据中0值所占比例
sum(spe==0) / (nrow(spe)*ncol(spe))
#请观察多度频率分布柱状图,如何解读为什么0值(缺失)在数据框内频
#率这么高?
造成缺失的因素有很多,但是有两种需要我们的注意:
1.真实的环境适合这个物种生存,只是我们采样的时候没有采集到
2.真实的环境不适合这个物种生产
对于零值我们需要小心处理,关键还是理解数据的生物学意义。
数据的进一步处理(分布位置)
数据探索也是数据和模型相互磨合的过程,不仅看数据还需要了解试验设计。
样方位置地图
# **************
# 生成空的绘图窗口(横纵坐标轴比例1:1,带标题)
# 从spa数据框获取地理坐标x和y
plot(spa, asp=1, type="n", main="样方位置",
xlab="x坐标 (km)", ylab="y坐标 (km)")
# 加一条连接各个样方点的蓝色线(代表Doubs河)
lines(spa, col="light blue")
# 添加每个样方的编号
text(spa, row.names(spa), cex=0.8, col="red")
# 添加文本
text(70, 10, "上游", cex=1.2, col="red")
text(20, 120, "下游", cex=1.2, col="red")
下面我们把物种数据映射到采样点上,看看物种是怎样随着河流变化的。
某些鱼类的分布地图
******************
将绘图窗口分割为4个绘图区域,每行两个
par(mfrow=c(2,2))
plot(spa, asp=1, col="brown", cex=speOMB, main="茴鱼",
xlab="x坐标 (km)", ylab="y坐标 (km)")
lines(spa, col="light blue")
plot(spa, asp=1, col="brown", cex=speBCO, main="欧鳊",
xlab="x坐标 (km)", ylab="y坐标 (km)")
lines(spa, col="light blue")
观察所生成的4张图,你就会明白为什么Verneaux 选择这4种鱼类作为不同区域的生态指示种,看了后面将要展示的环境因子空间分布情况会更清楚。
Tip:
理论小知识
par(mfrowc(1,2))实现一页多图的功能。
其中,通过设定函数par()的各个参数来调整我们的图形
mfrow=c(2,2)就是画4张画
mfrow=c(3,5)就是画14张画
例如,par(mfrow=c(2,3))一个图显示2行,3列
另一个引起我们注意的是plot()函数参数cex的用法,它的作用是定义数据点标识的大小。
下面我们解决每个物种在多少样方中出现?我们可以看看物种的相对频度
比较物种频度
**************
计算每个物种出现的样方数
按列进行计数,因此函数apply()第二个参数MARGIN应该设定为2
spe.pres <- apply(spe > 0, 2, sum)
按照升序的方式重新排列结果
sort(spe.pres)
计算频度百分比
spe.relf <- 100*spe.pres/nrow(spe)
round(sort(spe.relf), 1) # 设置排列结果为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="物种相对频度", right=FALSE, las=1,
xlab="出现率(%)", ylab="物种数量",
breaks=seq(0, 100, by=10), col="bisque")
提示:请关注函数apply()的使用,这里是对物种数据框spe的列进行汇总运算。注意第一部分spe>0表示先将spe内数值转化为逻辑向量TRUE/FALSE,然后对逻辑值TRUE进行列的计数汇总。
现在读者已经了解每个物种存在于多少个样方内,接下来要了解每个样方内有多少个物种存在(物种丰富度)
我们可以清楚地看出沿河流物种的整体分布。
最后我们用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