《数量生态学：R语言的应用》第二版第三章-关联测度与矩阵------Q模式

第三章-关联测度与矩阵------Q模式

生态学涉及多元统计方法，特别是排序和聚类，都是明确或不明确地基于所有可能对象或者变量之间的比较。这些比较通常采用关联测度（association
meansures）(常称为系数或者指数)的形式，不管是样方还是变量之间的比较都是基于他们组成的矩阵，因此选择合适的关联测度非常重要。

在任何分析之前，需要问下面这些问题：

比较对象（Q模式）,比较变量（R模式）

你是要处理物种数据（通常是非对称数据）还是其他类型的变量（对称数据）

你的数据是二元（二元系数）还是定量（数量系数），或这两类变量的混合，还是其他类型的数据（序数，特殊的系数）？
同时我们应该要知道，在Q模式中，关联测度是对象之间的相异或相似，例如欧氏距离、Jaccard相似系数。在R模式下，关联测度是变量之间的依赖性测度，例如协方差或相关系数。

注意：
1.测度(measure)、指数(index)、系数(coefficient)都表达相同的意思，目的都是量化对象或变量对之间的关系。
2.比较对象对的分析称为Q模式，关联测度是对象之间的相异或相似（如欧式距离、Jaccard相似系数）；比较变量对的称为R模式，关联测度是变量之间的依赖性（如协方差，相关系数）。

Q模式下对称或非对称的系数：双零问题

在两个对象中同一值为零，在这两个对象中可能蕴含的意义不同，但零值增加了对象的相似性。就物种数据而言，两个样方中都没有一个物种可能有不同的解释：不适合生存或者还没有迁徙到此地？因此物种存在的信息比物种缺失的信息更有意义。依据双零问题也可以区分两类关联测度：视双零为相似的依据（同其他值）的为对称系数，反之为非对称系数。在大部分情况下，应该优先选择非对称系数，除非可以确定引起双缺失的原因相同，例如在已知物种组成群落或生态同质区域内的控制实验。

这里可能比较难理解，简单来说就是两个样方都出现0值，但是造成0值的原因可能不一样，所以需要优先考虑非对称系数(除非可以确定引起0值的原因相同)

定性或定量数据的关联测度

变量可以是定性变量（名义的或分类的，二元的或多级的），也可以是半定量变量(序数的)或定量变量(离散的或连续的)。所以类型的变量均存在关联系数，其中大部分可以归为两类：二元变量的关联系数(二元系数，指被分析的变量是0-1的二元数据，并非关联测度数值为0-1的数据)和定量变量的关联系数(以下简称为数量系数)

image

Q模式：计算对象之间的相异矩阵

在Q模式分析中，我们需要用到6个程序包：stats(安装基础程序时已经载入)、vegan、ade4、adespatial、cluster和FD等等。

在R中，所有的相似测度方阵可以转化为相异测度方阵，距离方阵(R里面属于"dist"类对象)对角线的值(每个对象与自身的距离)均为0

Q模式：定量的物种数据

定量的物种数据通常需要使用非对称的距离测度。在物种数据分析方面，常用的系数有Bray-Curtis相异系数、弦（chord）距离、Hellinger距离和卡方距离。

Bray-Curtis 相异度(百分数差异)：Bray-Curtis相异度是生态学中用来衡量不同样地物种组成差异的测度，可以计算生物样本中不同物种组成的数量特征。包括：多度，盖度，重要值等。优点：Bray-curtis计算时，不仅考虑样本中物种的有无，而且还考虑不同物种的相对丰度。一般常规的是使用原始数据，然而，生物多样性中的数据较为复杂，数值之间存在很大的极差，为了消除这一影响，采用的方法是进行取对数。

弦（chord）距离：先将样方向量的范数标准化后(平方和为1)再计算样方对的欧氏距离；可以用vegan包内decostand()函数处理（选择参数normalize）。也可以用adespatial包内的dist.ldc()函数通过选择参数"chord"来获得弦距离。

Hellinger距离：除以样方多度总和再取平均值后计算欧氏距离，是平方根转化后的弦转化。可以用vegan包内decostand()函数处理（选择参数hellinger）。也可以用adespatial包内的dist.ldc()函数通过选择参数"hellinger"来获得hellinger距离。

log-chord距离是先将数据对数化再求弦距离。首先对原始数据加1后取自然对数（ln(y+1)）,然后对数化后的数据进行弦转化后求欧氏距离。log-chord距离可以用adespatial包内的dist.ldc()函数通过选择参数"log.chord"来获得。

在R中实现：

setwd("E:/Rstudio_working/lai jian shan/DATA/Data")
#加载包、函数和数据
library(ade4)
library(adespatial)
library(vegan)#可以先加载ade4再加载vegan，不然可能会有一些问题
library(gclus)
library(cluster)
library(FD)

#加载本章后面部分将要用到的函数
source("coldiss.R")
source("panelutils.R")

#导入Doubs数据
load("Doubs.Rdata")
#剔除无物种数据的样方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]

###Q模式相异矩阵
##  Q模式定量(半定量)数据的相异和距离测度
#原始物种数据的百分数差异(也称为Bray-Curtis)相异矩阵
#".db"意味着在Bray-Curtis距离

spe.db <- vegdist(spe)  #  Bray-Curtis相异系数（默认）
spe.db

# 对数转化后物种数据的Bray-Curtis相异矩阵
spe.dbln <- vegdist(log1p(spe)) #log1p(spe)=log(spe+1),对0值加1，避免有缺失值

# 弦距离矩阵
spe.dc <- dist.ldc(spe,"chord")
#在vegan包里面分两步计算
# 弦距离矩阵
spe.norm <- decostand(spe, "nor")
spe.dc <- dist(spe.norm)

# Hellinger距离矩阵
spe.dh <- dist.ldc(spe)
#dist.ldc函数中Hellinger是默认的距离模式
#在vegan包分两步计算
spe.hel <- decostand(spe, "hel")
spe.dh <- dist(spe.hel)

#对数转化后弦距离矩阵
spe.logchord <- dist.ldc(spe,"log.chord")
#在vegan里面分三步计算
   spe.ln <- log1p(spe)
   spe.ln.norm <- decostand(spe.ln, "nor")
   spe.logchord <- dist(spe.ln.norm)

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

当可用的仅仅是二元(有-无)物种数据，或多度的数据不适用，或包含不确定的定量数据时，可使用有-无(0-1)数据进行分析。

Jaccard相异矩阵：
对于每个样方对，Jaccard系数是两个样方共有物种数目(交集)除以两个样方所含有的全部物种数(并集)。
Sφrensen相异矩阵：Sφrensen相似系数给予两个样方共有物种数双倍的权重，其反数(1-S8)等价于基于物种有-无数据计算的百分数差异(又称Bray-Curtis)相异系数
Ochiai相异矩阵：基于物种有-无的数据算出弦距离、Hellinger距离和对数弦距离

# Q模式---二元数据的相异测度

# 使用vegdist（）函数计算Jaccard相异矩阵
spe.dj <- vegdist(spe, "jac", binary=TRUE)
#head(spe.dj) #head函数查看数据我都报错，所以用#注释了,可以用spe.dj 直接输出
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.ldc（）函数计算Sorensen相异矩阵
spe.ds <- dist.ldc(spe, "sorensen")
#head(spe.ds)

# 使用vegdist（）函数计算Sorensen相异矩阵
spe.ds2 <- vegdist(spe,method="bray",binary=TRUE)

# 使用dist.binary（）函数计算Sorensen相异矩阵
spe.ds3 <- dist.binary(spe,method=5)
#head(spe.ds)
#head(spe.ds2)
#head(sqrt(spe.ds2))
#head(spe.ds3)

# Ochiai相异矩阵
spe.och <- dist.ldc(spe,"ochiai") #或
spe.och <- dist.binary(spe, method=7)
#head(spe.och)
#这里显示Jaccard和Sorensen距离矩阵有两种数值。

#查看vegdist()、dist.binary()和dist()三个函数的帮助文件，了解这些函数能计算哪些系数。

注意：所有的二元距离函数在计算系数时，均会自动对数据进行二元转化
因此这里的数据不需要二元转化（decostand(，"pa")）。函数dist.binary（）会自动对数据进行二元转化，但函数vegist（）需要设定参数binary=TRUE。
在dist.ldc()中默认binary=FALSE，除非选择的距离参数是Jaccard和Sφrensen和Ochiai指数。

关联矩阵一般作为中间实体，很少用于直接研究。然而，如果对象不多，直接展示关联矩阵也很有用，能够将数据的主要特征可视化。

建议使用coldiss()函数可视化相异矩阵。coldiss()函数会使用一个能重新排列矩阵的函数order.single()(属于gclus包)，该函数可以根据对象之间的距离沿着对角线重新将对象排位。但是必须先安装gclus包。

## 图解关联矩阵
#使用coldiss()函数生成热图或格状图
 
#coldiss()函数使用说明：
#coldiss(D=dist.object,nc=4,byrank=TRUE,diag=FALSE)
#D应该是相异矩阵
#如果D为最大值大于1的距离矩阵，此时D会除以max(D)
#nc颜色种类数量
#byrank=TRUE   等大小分级，即每个颜色所包含的值的数量一样多
#byrank=FALSE  等区间分级，即每个颜色所包含的值的区间一样长
#如果diag=TRUE 表示样方号放置在矩阵对角线上

## 比较从物种数据获得差异矩阵和距离矩阵

#等区间分级的4种颜色(方便比较)

#图解基于原始鱼类多度数据的百分数差异(也称为Bray-Curtis)相异矩阵
#矩阵
coldiss(spe.db,byrank = FALSE,diag = TRUE)

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

基于原始数据的百分数相异矩阵热图

基于对数转化后数据的百分数相异矩阵热图

在未转化的相异矩阵中，数量多的物种之间的多度差异与数量少的物种之间的多度差异有同等权重。

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

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

#对数转化后的弦距离矩阵

coldiss(spe.logchord,byrank = FALSE,diag = TRUE)

弦距离矩阵

Hellinger距离矩阵

对数转化后的弦距离矩阵

#Jaccard 相异矩阵
coldiss(spe.dj, byrank = FALSE, diag = TRUE)

我们可以比较当前的Jaccard距离热图和前面的其他距离矩阵热图。Jaccard图是基于二元数据计算。我们可以仔细观察Jaccard图和前面的数量系数图之间的差异是否比数量系数图之间的差异大呢？

这些案例均是处理物种数据，Doubs样带具有强烈的生态梯度特征(例如氧含量和硝酸盐浓度)。Doubs样带的环境背景很清楚，可以假设在特定的某一段河流，物种的缺失可能是某种相同的原因造成的，因此可以计算对称系数的关联矩阵。

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

Q模式：定量数据(除物种多度数据外的数据)

对双零有明确解释的定量数据，欧氏距离是对称距离测度的最佳选择。注意欧氏距离的值没有上限，但受变量纲量影响较大，所以前面我们的数据转化就派上用场了。

此处用标准化后的环境因子变量（env）计算样方的欧氏距离。先剔除dfs变量(离源头距离)，因为它属于空间变量而非环境因子变量。同样使用coldiss()函数可视化距离矩阵。

#剔除env数据框内dfs变量
env2 <- env[, -1]  #删除第一列的数据，然后赋值给新变量

#由标准化后的env2数据框计算的欧氏距离矩阵
env.de <- dist(scale(env2))
coldiss(env.de, nc = 16, diag = TRUE)

标准化后的env2数据框计算的欧氏距离矩阵

注意：可以利用scale()函数对环境变量进行快速标准化

相异矩阵的热图很合适快速比较，例如，可以同时绘制基于物种多度和基于环境因子的Hellinger距离图，为了便于比较，两个图均选择等数量的分级

#物种数据的Hellinger距离矩阵
#(nc=16等数量的分级)
coldiss(spe.dh,nc=16,diag = TRUE)

物种数据的Hellinger距离矩阵

欧氏距离理所应当然可以用于计算基于地理坐标变量的地理距离矩阵。
地理坐标可以是一维或者二维的直角坐标系(笛卡儿坐标)，其单位也可以多种多样(例如cm、m、km属于相同投影带的UTM坐标)。如果是球体系统坐标(经纬度)，在计算欧氏距离之前必须先转化。SoDA程序包内geoXY()函数可以完成球坐标系统的专业。需要注意的是，标准化数据会改变两个维度的比率，因此一般地理坐标(x-y)不应该标准化(如果需要可以标准化)

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

#基于一维dfs变量(离源头距离)的欧氏距离矩阵
dfs.df <- as.data.frame(env$dfs,row.names = rownames(env))
riv.de <- dist(dfs.df)
coldiss(riv.de,nc=16,diag = TRUE)

基于二维空间坐标的欧氏距离矩阵

基于一维dfs变量(离源头距离)的欧氏距离矩阵

Q模式：二元数据（除物种有-无数据外的数据）

对于二元数据，最简单的对称相似测度是"简单匹配系数S1"，对于每组样方，S1是双1的数量加上双0的数量除以变量数。

由于当前的环境变量是特定的定量变量，不能简单转化为二元数据。此处需要自己创建虚拟数据来演示S1的计算。随机产生数据方法会经常用到，这里我们也练习如何在R内生成满足不同模拟需求的数据集。

创建随机数的方法

#基于随机数据的案例
#生成30个对象、5个二元变量的数据集，每个变量有预先设置的固定的0和1的数量
#变量1：10个1和20个0，顺序随机
var1 <- sample(c(rep(1,10),rep(0,20)))
var1
# 变量2：15个0在一起，15个1在一起
var2 <- c(rep(0,15),rep(1,15))
var2
#变量3：3个1，3个0交替出现，直到总数量达到30为止
var3 <- rep(c(1,1,1,0,0,0),5)
var3
#变量4：5个1，10个0交替出现，直到总数量达到30为止
var4 <- rep(c(rep(1,5),rep(0,10)),2)
var4
#变量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)
var5
#将变量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 模式：混合类型，包括分类（定性多级）变量

Gower相似系数当作一种对称指数；当数据框内一个变量被当做一个因子时，最简单的匹配规则被应用，即如果一个因子在两个对象中有相同的水平，表示该对象对相似指数为1，反之为0。Gower相异指数可以利用cluster程序包内daisy()函数计算。应避免使用vegdist()函数计算Gower相异系数，因此该函数只适用与定量数据和有-无数据计算，对多级变量并不适用。

只要每个变量给予合适的定义，daisy()函数就可以处理混合变量的数据。当数据中存在缺失值时，该函数会自动排除与含有缺失值样方对的计算。

FD程序包里gowdis()函数是计算Gower相似系数最完善的函数，可以计算混合变量(包括非对称的二元变量)的距离，也可以像daisy()函数一样设置变量的权重和处理缺失值。

#为计算Gower指数(S15)的虚拟数据
#随机生成30个平均值为0、标准差为1的正态分布数据
var.g1 <- rnorm(30,0,1)
var.g1
#随机生成30个从0到5均匀分布的数据
var.g2 <- runif(30,0,5)
var.g2
#生成3个水平的因子变量(每个水平10个重复)
var.g3 <- gl(3,10,labels=c("A","B","C"))
var.g3
#生成与var.g3正交的两个水平的因子变量
var.g4 <- gl(2,5,30,labels=c("D","E"))
var.g4

var.g3 和var.g4 组合在一起代表一个双因素交叉平衡设计

dat2 <- data.frame(var.g1,var.g2,var.g3,var.g4)
summary(dat2)

dat2数据集

summary 查看

# 使用daisy（）函数计算Gower相异矩阵
# 完整的Gower相异矩阵（基于4个变量）
dat2.S15 <- daisy(dat2, "gower")
range(dat2.S15)
coldiss(dat2.S15, diag=TRUE)

完整的Gower相异矩阵（基于4个变量）

# 仅使用2个正交的因子变量计算gower相异矩阵
dat2partial.S15 <- daisy(dat2[,3:4], "gower")
coldiss(dat2partial.S15, diag=TRUE)

# 在dat2partial.S15矩阵内相异系数值代表什么？
levels(factor(dat2partial.S15))
#对象对所对应的数值分享共同的因子水平2、1和无因子水平。最高相异系数
#值的对象对不分享共同的因子水平。

仅使用2个正交的因子变量计算gower相异矩阵

# 使用FD程序包内gowdis（）函数计算Gower相异矩阵
library(FD)  #如果FD还没载入
?gowdis
dat2.S15.2 <- gowdis(dat2)
range(dat2.S15.2)
coldiss(dat2.S15.2, diag=TRUE)