基于距离的Moran特征根图(dbMEN / PCNM)
通过主坐标矩阵模拟空间结构,进而作为解释变量解释群落组成等。如果解释率很高,可结合环境因子及变差分解进一步探讨,是因为响应变量自身的空间结构或者未测量的环境变量(独自解释的部分),还是已经测量的环境因子(共享部分)具有这样的空间结构
- 需要注意的是这个分析基于经纬度(要先把GPS所得坐标转换为笛卡尔坐标),具体转换方法:
# PCNM,现在改为dbMEN-------------------------------------------------------------------------
library(SoDA)
a=read.csv("latlong.csv",row.names=1)
xy=geoXY(a$LatitudeN,a$LongitudeE)#转换为笛卡尔坐标系后用于后续分析
此处要感谢张金龙老师和赖江山老师的解惑,后面再次翻阅赖老师的数量生态学时发现其实书中明确给出了这个信息,还是看得不够认真😅
R语言中具体的实现步骤:
library(adespatial)
data("mite.xy")#加载数据
mite.hel <- decostand(mite,method = "hellinger")#物种数据转化
anova(rda(mite.hel,mite.xy))#结果显著,进行去趋势
mite.hel.det <- resid(lm(as.matrix(mite.hel)~.,data = mite.xy))#相当于进行线性回归,保留残差
mite.dbmen <- dbmem(mite.xy,silent = FALSE)
# 基于去趋势数据构建全模型
mite.dbmen.rda <- rda(mite.hel.det~.,mite.dbmen)
anova(mite.dbmen.rda)
# 若模型显著,则计算校正R2,进行前向选择
(mite.r2a <- RsquareAdj(mite.dbmen.rda)$adj.r.squared)
mite.dbmen.fwd <- forward.sel(mite.hel.det,as.matrix(mite.dbmen),adjR2thresh = mite.r2a)
nrow(mite.dbmen.fwd)#筛选所得显著的dbmen数量
(mite.sign <- sort(mite.dbmen.fwd$order))#排序
mite.dbmen.red <- mite.dbmen[,mite.sign]#提取显著dbmen
# 对筛选后的变量进行检验
mite.dbmen.rda2 <- rda(mite.hel.det~.,mite.dbmen.red)
anova(mite.dbmen.rda2)
RsquareAdj(mite.dbmen.rda2)$adj.r.squared#校正R2未明显下降
anova(mite.dbmen.rda2,by="axis")
至此分析就结束了,其中的mite.dbmen.red
就是我们需要的数据,可通过write.csv()
等导出至本地用于后续分析
一个简便的function,一步实现上述过程
source("quickMEM.R")
library(ade4)
mite.dbmen.quick <- quickMEM(mite.hel,mite.xy)
mite.dbmen.var <- mite.dbmen.quick[[1]][,mite.dbmen.quick[[3]][,1]]#提取最终所得的轴坐标,用于后续分析
如果报错?有可能是没有加载所需的程序包,例如ade4
quickMEM.R(点击即可下载,也可以去网上寻找数量生态学的代码附件)
对于嵌套数据而言:
#嵌套数据(如不同湖泊的样点同时分析)
?create.dbMEM.model#嵌套数据需使用该功能指定样本具体所属组别
create.dbMEM.model(coord=cartesian, nsites=nsites.per.group)
其中coord为笛卡尔坐标,nsites为每个组(同一组的所有样本需放在一起)的样本数量
# -------------------------------------------------------------------------
# 环境变量数据整理,将分类变量变成0,1数值矩阵(由于变差分解varpart()功能只能接受定量变量)
substrate <- model.matrix(~mite.env[,3])[,-1]
shrub <- model.matrix(~mite.env[,4])[,-1]
topo <- model.matrix(~mite.env[,5])[,-1]
mite.env2 <- cbind(mite.env[,c(1,2)],substrate,shrub,topo)
colnames(mite.env2) <- c("SubsDens","WatrCont","Sphagn2","Sphagn3","Sphagn4",
"Litter","Barepeat","Interface","Few","Many","Topo")
#这个哑变量设置过程有点没弄懂,但是不影响这个主题的分析过程,继续往下
mite.env.rda <- rda(mite.hel~.,mite.env2)
mite.env.r2a <- RsquareAdj(mite.env.rda)
mite.env.fwd <- forward.sel(mite.hel,mite.env2,adjR2thresh = mite.env.r2a)
mite.env.sign <- sort(mite.env.fwd$order)
mite.env.red <- mite.env2[,mite.env.sign]
#方法1:
res.part <- varpart(mite.hel,mite.dbmen.red,mite.env.red)
plot(res.part,Xnames=c("spatial","env"),bg=c("red", "lightgreen"))
anova(rda(mite.hel,mite.env.red,mite.dbmen.red))
anova(rda(mite.hel,mite.dbmen.red,mite.env.red))
#方法2:
library(rdacca.hp)
?rdacca.hp
iv <- list(mite.env.red,mite.dbmen.red)
(mite.var <- rdacca.hp(mite.hel,iv,method = "RDA",var.part = TRUE))
permu.hp(mite.hel,iv,method = "RDA")
plot(mite.var)
这篇推文就懒得上传图片了。。。。
参考来源:
[1].数量生态学(第二版) 赖江山译
[2].The accumulation of species and recovery of species composition along a 70 year succession in a tropical secondary forest(文章附有代码),参见作者的解读:一个群落beta多样性分析的完整R脚本