生命科学-简书专题统计分析与数据挖掘R语言与统计分析

《数量生态学------R语言的应用》第二版阅读笔记1

2021-04-19  本文已影响0人  fafu生信小蘑菇

赖江山老师《数量生态学------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. 真实的环境不适合这个物种生存,在这里生存就会被淘汰。 所以对于零值我们要小心处理,关键还是理解数据的生物学意义。

物种数据:进一步分析------生成样方位置分布图

#样方位置地图
#生成空的绘图底板(横纵坐标轴比例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="双变量散点图(带频度分布图和平滑曲线)")

注意:每个散点图展示对角线上两个变量之间的关系。散点图的横坐标对应上方或下方的变量,纵坐标对应左侧或右侧的变量

好了绪论和第二章的一部分就到这里,下一次给大家带来第二章的数据探索,这次的部分画图数据处理是经过数据处理的,详细请看第一次推送!!!

感谢你的阅读,有问题可以联系小蘑菇。

上一篇下一篇

猜你喜欢

热点阅读