数量生态学笔记|数量生态学笔记||Q模式 Q模式|

文章图片
好了，回到我们《数量生态学》课堂上来。
数量生态学笔记||绪论
数量生态学笔记||数据探索
数量生态学笔记||物种数据转化
在进行了简单的数据探索之后，在对数据进行了符合模型的计算之后，是不是马上可以绘制Nature Science上面的高颜值图表了呢？在这之前我们还有一段路线要走，要记住我们关注的始终是生态学问题而不是别的什么？
生态学涉及许多多元统计方法，特别是排序和聚类都是明确或不明确地基于所有可能对象或者变量之间的比较。这些比较通常采用关联测度（association meansures）(常称为系数或者指数)的形式，不管是样方还是变量之间的比较都是基于他们组成的矩因此选择合适的关联测度非常重要。
在任何分析之前，跟踪你所需的关联测度，需要问下面这些问题：

你是要比较对象（Q模式）还是比较变量（R模式）
你是要处理物种数据（通常是非对称数据）还是其他类型的变量（对称数据）
你的数据是二元（二元系数）还是定量（数量系数），或是这两类变量的混合，还是其他类型的数据（序数，特殊的系数）？

关联测度的主要类别关联测度不管是距离、指数、系数目的都是量化对象对或变量对之间的关系。但是比较对象对和变量对往往需要不同的测度，按照比较对象的不同我们吧关联测度分为两类：比较对象对的分析称为Q模式，关联测度是对象对之间的距离（如欧式距离）；比较变量对的称为R模式，关联测度是变量之间的依赖性（如协方差，相关系数）。
Q模式下对称或非对称的系数：双零问题在两个对象中同一值为零，在这两个对象中可能蕴含的意义不同，但零值增加了对象的相似性。就物种数据而言，两个样方中都没有一个物种可能有不同的解释：不适合生存还是为迁徙到此地？因此物种存在的信息比物种缺失的信息更有意义。依据双零问题也可以区分两类关联测度：视零为一样的（同其他值）的为对称系数，反之为非对称系数。
Q模式：定量的物种数据（非对称的距离测度）

Bray-Curtis 相异度：多度绝对数量差

文章图片
弦（chord）距离：范数标准化在计算欧氏距离
Hellinger距离：除以样方多度总和再取平均值后计算欧氏距离

在R中我们有具体的函数来实现：

# 加载所需程序包 library(ade4) library(vegan)# 应该先加载ade4再加载vegan，以避免一些冲突 library(gclus) library(cluster) library(FD) # 导入CSV格式的数据setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA") spe <- read.csv("DoubsSpe.csv", row.names=1) env <- read.csv("DoubsEnv.csv", row.names=1) spa <- read.csv("DoubsSpa.csv", row.names=1) # 剔除无物种数据的样方8 spe <- spe[-8,] env <- env[-8,] spa <- spa[-8,] # 定量（半定量）数据的相异和距离测度 # ********************************** # 原始物种数据的Bray-Curtis相异矩阵 spe.db <- vegdist(spe)# Bray-Curtis相异系数（默认） head(spe.db) # 对数转化后物种数据的Bray-Curtis相异矩阵 spe.dbln <- vegdist(log1p(spe)) head(spe.dbln) # 弦距离矩阵 spe.norm <- decostand(spe, "nor") spe.dc <- dist(spe.norm) head(spe.dc) # Hellinger距离矩阵 spe.hel <- decostand(spe, "hel") spe.dh <- dist(spe.hel) head(spe.dh) #查看decostand（）函数和vegdist（）函数的帮助文件，寻找计算卡方距离#矩阵的参数设定，卡方距离对当前这些数据也适用。

Q模式：二元（有无）物种数据：不考虑物种丰度

Jaccard相异矩阵
S?rensen相异矩阵
Ochiai相异矩阵

# 二元数据的相异测度 # ****************** #注意：所有的二元距离函数在计算系数时，均会自动对数据进行二元转化 #因此这里的数据不需要二元转化（decostand(，"pa")）。函数dist.binary（） #会自动对数据进行二元转化，但函数vegist（）需要设定参数binary=TRUE。 # 使用vegdist（）函数计算Jaccard相异矩阵 spe.dj <- vegdist(spe, "jac", binary=TRUE) head(spe.dj) head(sqrt(spe.dj)) # 使用dist（）函数计算Jaccard相异矩阵 spe.dj2 <- dist(spe, "binary") head(spe.dj2) # 使用dist.binary（）函数计算Jaccard相异矩阵 spe.dj3 <- dist.binary(spe, method=1) head(spe.dj3) # 使用dist.binary（）函数计算S?rensen相异矩阵 spe.ds <- dist.binary(spe, method=5) head(spe.ds) # 使用vegdist（）函数计算S?rensen相异矩阵 spe.ds2 <- vegdist(spe, binary=TRUE) head(spe.ds2) head(sqrt(spe.ds2)) # Ochiai相异矩阵 spe.och <- dist.binary(spe, method=7) head(spe.och) #这里显示Jaccard和S?rensen距离矩阵有两种数值。可以回看3.3节引言部 #分了解这两种数值的差异。

# 图解关联矩阵 # ************ # 运行coldiss（）函数时会自动要求加载gclus程序包，我们也可以预先加载 library(gclus) # 导入我们自己编写的函数coldiss（） source("coldiss.R")#如果函数文件没有在当前工作目录下，需要指定文件路径 # 使用coldiss（）函数产生的彩图（在一些数据分析文献中也称作热图或格状图） # ****************************************************************** #绘制两种（原始和重排）相异矩阵图并标色。屏幕输出颜色分类：红紫色= #相异系数接近0（最大相似系数）；青绿色=相异系数接近1（最小相似系数） #coldiss（）函数使用说明： # coldiss(D=dist.object, nc=4, byrank=TRUE, diag=FALSE) # D应该是一个相异矩阵 #如果D为最大值大于1的距离矩阵，此时D会除以max(D) # nc颜色种类数量 # byrank=TRUE等大小分级，即每个颜色所包含的值的数量一样多 # byrank=FALSE等区间分级，即每个颜色所包含的值的区间一样长 #如果diag=TRUE，表示样方号放置在矩阵对角线上 # 图解基于原始鱼类多度数据的Bray-Curtis相异矩阵 # 等区间分级的4种颜色（方便比较）

coldiss(spe.db, byrank=FALSE, diag=TRUE) # 图解基于对数转化数据的Bray-Curtis相异矩阵

文章图片

# 图解基于对数转化数据的Bray-Curtis相异矩阵 coldiss(spe.dbln, byrank=FALSE, diag=TRUE)

文章图片

# 弦距离矩阵 coldiss(spe.dc, byrank=FALSE, diag=TRUE)

文章图片

# Hellinger距离矩阵 coldiss(spe.dh, byrank=FALSE, diag=TRUE)

文章图片

# Jaccard距离矩阵 coldiss(spe.dj, byrank=FALSE, diag=TRUE) #请比较当前的Jaccard距离热图和前面的其他距离矩阵热图。Jaccard图是 #基于二元数据计算。这是否影响结果呢？Jaccard图和前面的数量系数图之 #间的差异是否比数量系数图之间的差异大呢？

文章图片

# 简单匹配相异系数（在ade4程序包内也称为Sokal & Michener指数) spe.s1 <- dist.binary(spe, method=2) coldiss(spe.s1^2, byrank=FALSE, diag=TRUE) #比较一下当前的对称相异矩阵与之前的Jaccard矩阵。哪个受双零问题影响 #更突出？

文章图片
Q 模式：定量数据（除物种丰度数据外的数据）对双零有明确解释的定量数据，欧氏距离是对称距离测度的最佳选择。注意欧氏距离受变量纲量影响较大，所以前面我们的数据转化就派上用场了。
此处用标准化后的环境因子变量（env）计算样方的欧氏距离。

# 剔除env数据框内das空间变量 env2 <- env[,-1] # 由标准化后的env2数据框计算的欧氏距离矩阵 env.de <- dist(scale(env2)) coldiss(env.de, diag=TRUE)

文章图片

# 物种数据的Hellinger距离矩阵（等数量的分级） coldiss(spe.dh, diag=TRUE) #因为具有相同的样方位置排序，可以比较图3.2的左边热图与前面基于环 #境变量的热图。你是否能观察到一些共同的特征？

文章图片
欧式距离理所当然可用于计算地理坐标变量的地理矩阵。

# 基于二维空间坐标的欧氏距离矩阵 spa.de <- dist(spa) coldiss(spa.de, diag=TRUE)

文章图片

# 基于一维das变量（离源头距离）的欧氏距离矩阵 das.df <- as.data.frame(env$das, row.names=rownames(env)) riv.de <- dist(das.df) coldiss(riv.de, diag=TRUE) #为什么基于x-y的欧氏距离图和基于das的欧氏距离图有这样的差异？

文章图片
Q 模式：二元数据（处物种有无数据外的数据）对于二元数据，最简单的对称相似测度是“简单匹配系数S1”，对于每组样方，S1是双1的数量加上双0的数量除以变量数。

# 生成30个对象、5个二元变量的数据集，每个变量有预先设置的固定的0和1的 # 数量 # 变量1：10个1和20个0，顺序随机 var1 <- sample(c(rep(1,10), rep(0,20))) # 变量2：15个0在一起，15个1在一起 var2 <- c(rep(0,15), rep(1,15)) # 变量3：3个1，3个0交替出现，直到总数量达到30为止 var3 <- rep(c(1,1,1,0,0,0),5) # 变量4：5个1，10个0交替出现，直到总数量达到30为止 var4 <- rep(c(rep(1,5), rep(0,10)), 2) # 变量5：前16元素是7个1和9个0的随机排列，接着是4个0和10个1 var5.1 <- sample(c(rep(1,7), rep(0,9))) var5.2 <- c(rep(0,4), rep(1,10)) var5 <- c(var5.1, var5.2) # 将变量1至变量5合成一个数据框 dat <- data.frame(var1, var2, var3, var4, var5) dim(dat) # 简单匹配系数的计算（在ade4程序包中也称为Sokal & Michener指数） dat.s1 <- dist.binary(dat, method=2) coldiss(dat.s1, diag=TRUE)

文章图片
Q 模式：混合类型，包括分类（定性多级）变量

Cower相似系数

# 为计算Gower指数（S15）的虚拟数据 # 随机生成30个平均值为0、标准差为1的正态分布数据 var.g1 <- rnorm(30,0,1) # 随机生成30个从0到5均匀分布的数据 var.g2 <- runif(30,0,5) # 生成3个水平的因子变量（每个水平10个重复） var.g3 <- gl(3,10) # 生成与var.g3正交的2个水平的因子变量 var.g4 <- gl(2,5,30) #var.g3和var.g4组合在一起代表一个双因素交叉平衡设计 dat2 <- data.frame(var.g1,var.g2,var.g3,var.g4) summary(dat2) # 使用daisy（）函数计算Gower相异矩阵 # 完整的Gower相异矩阵（基于4个变量） dat2.S15 <- daisy(dat2, "gower") range(dat2.S15) coldiss(dat2.S15, diag=TRUE)

文章图片

# 仅使用2个正交的因子变量计算gower相异矩阵 dat2partial.S15 <- daisy(dat2[,3:4], "gower") coldiss(dat2partial.S15, diag=TRUE) # 在dat2partial.S15矩阵内相异系数值代表什么？ levels(factor(dat2partial.S15)) #对象对所对应的数值分享共同的因子水平2、1和无因子水平。最高相异系数 #值的对象对不分享共同的因子水平。

文章图片

# 使用FD程序包内gowdis（）函数计算Gower相异矩阵 library(FD)#如果FD还没载入 ?gowdis dat2.S15.2 <- gowdis(dat2) range(dat2.S15.2) coldiss(dat2.S15.2, diag=TRUE) # 仅使用两个正交的因子变量计算距离矩阵 dat2partial.S15.2 <- gowdis(dat2[,3:4])