vegan

---

title: “vegan study”
author: “hardy huang”
date: “2021/3/31”
output: html_document

knitr::opts_chunk$set(echo = TRUE)

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

Including Plots

You can also embed plots, for example:

plot(pressure)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

#

#

library(vegan)
library(MASS)
data("varespec")
vare.dis<-vegdist(varespec)
vare.mds0<-isoMDS(vare.dis)

# stressplot()评估拟合结果
stressplot(vare.mds0,vare.dis)

ordiplot(vare.mds0,type = "t")

#metaMDS迭代算法
vare.mds<-metaMDS(varespec,trace=FALSE)
vare.mds
plot(vare.mds,type = "t")

#2.3比较排序图：普氏旋转

tmp<-wisconsin(sqrt(varespec))
dis<-vegdist(tmp)
vare.mds0<-isoMDS(dis,trace = 0)
pro<-procrustes(vare.mds,vare.mds0)
pro

plot(pro)
plot(pro,kind = 2)

#2.4 特征向量法

vare.pca<-rda(varespec)
vare.pca
plot(vare.pca)
sum(apply(varespec, 2, var))#总方差
biplot(vare.pca,scaling = -1)
#将所有物种标准化为单位方差
vare.pca<-rda(varespec,scale = TRUE)
vare.pca
biplot(vare.pca)
dim(varespec)

对应分析

vare.ca<-cca(varespec)
vare.ca
plot(vare.ca)
chisq.test(varespec/sum(varespec))
plot(vare.ca,scaling = 1)

#2.5去趋势对应分析（DCA）

vare.dca<-decorana(varespec)
vare.dca
plot(vare.dca,display = "sites")

#2.6排序图

data("BCI")
head("BCI")
mod<-decorana(BCI)
plot(mod)
#变量名太长使用make.cepnames(),缩写拉丁字母
shnam<-make.cepnames(names(BCI))
shnam[1:5]

pl<-plot(mod,dis="sp")
identify(pl,"sp",labels = shnam)
stems<-colSums(BCI)
plot(mod,dis="sp",type="n")
sel<-orditorp(mod,dis="sp",lab=shnam,priority = stems,pcol="gray",pch="+")

#ordilabel()函数可以在不透明的标签上绘制文本，覆盖在下面的其他标签之上


plot(mod,dis="sp",type = "n")
ordilabel(mod,dis="sp",labels = shnam,priority = stems)

#3使用环境因子被动解释排序轴

data(varechem)
ef<-envfit(vare.mds,varechem,permu=999)
ef
plot(vare.mds,display="sites")
plot(ef,p.max=0.1)

#3.2曲面拟合

ef<-envfit(vare.mds~Al+Ca,varechem)
plot(vare.mds,display = "sites")
plot(ef)

temp<-with(varechem,ordisurf(vare.mds,Al,add=TRUE))
with(varechem,ordisurf(vare.mds,Ca,add = TRUE,col="green4"))

#3.3因子变量

data(dune)
data(dune.env)
dune.ca<-cca(dune)
ef<-envfit(dune.ca,dune.env,permutations = 999)
ef
plot(dune.ca,display="sites")
plot(ef)

#

plot(dune.ca,display = "sites",type="p")
with(dune.env,ordiellipse(dune.ca,Management,kind = "se",conf =0.95,col="red"))
with(dune.env,ordispider(dune.ca,Management,col = "blue",label = TRUE))
with(dune.env,ordihull(dune.ca,Management,col = "black",lty = 2))

#

#4.1模型说明

vare.cca<-cca(varespec~Al+P+K,varechem)
vare.cca
plot(vare.cca)
library(vegan3d)
ordiplot3d(vare.cca,type="h")
dune.cca<-cca(dune~Management,dune.env)#因子变量
plot(dune.cca)
dune.cca

vare.cca<-cca(dune~Moisture,dune.env)
plot(vare.cca)

#

anova(vare.cca)
mod<-cca(varespec~Al+P+K,varechem)
anova(mod,by="term",step=200)
anova(mod,by="margin",perm=500)
anova(mod,by="axis",perm=1000)

#

mod1<-cca(varespec~.,varechem)
mod1
plot(procrustes(cca(varespec),mod1))
# 减少约束变量，找到重要的环境变量
mod0<-cca(varespec~1,varechem)
mod<-step(mod0,scope = formula(mod1),test="perm")
mod

#

spenvcor(mod)
# WA坐标和相应的LC坐标相结合
dune.cca<-cca(dune~Management,dune.env)
plot(dune.cca,display = c("lc","wa"),type = "p")
ordispider(dune.cca,col = "blue")

#

# vegan双序图
# 当scaling=2时，双序图箭头和样方的关系最准，
# 当scaling=1时，他们和物种的关系最好。
pred<-calibrate(mod)
head(pred)
with(varechem,plot(Al,pred[,"Al"]-Al,ylab = "prediction error"))#预测的残差图
abline(h=0,col="grey")
plot(mod,display = c("bp","wa","lc"))
ef<-with(varechem,ordisurf(mod,Al,display = "lc",add = TRUE))#Al的投影
ef<-with(varechem,ordisurf(mod,P,display = "lc",add = TRUE))#P的投影

#

dune.cca<-cca(dune~Management+Condition(Moisture),dune.env)
plot(dune.cca)
dune.cca
anova(dune.cca,perm.max=2000)
anova(cca(dune~Management,dune.env))
with(dune.env,anova(dune.cca,strata = Moisture))

#
#

betad<-betadiver(dune,"z")
adonis(betad~Management,dune.env,perm=200)
adonis(betad~A1*Management,dune.env,perm=200)

#

mod<-with(dune.env,betadisper(betad,Management))
mod
plot(mod)
boxplot(mod)
anova(mod)
permutest(mod)

#

pc<-prcomp(varechem,scale=TRUE )
pc<-scores(pc,display="site",choices=1:4)
edis<-vegdist(pc,method = "euclid")
vare.dis<-vegdist(wisconsin(sqrt(varespec)))
mantel(vare.dis,edis)
# 绘制相异矩阵散点图
plot(vare.dis,edis)

#

# 现在比较metaMDS排序结果和环境变量PCA分析的前两轴
pc<-scores(pc,choices = 1:2)
pro<-protest(vare.mds,pc)
plot(pro)
pro

#6分类
#6.1聚类分析

dis<-vegdist(dune)
clus<-hclust(dis,"single")#单连接聚类
plot(clus)

cluc<-hclust(dis,"complete")#完全连接聚类
plot(cluc)


clua<-hclust(dis,"average")#平均距离聚类
plot(clua)

range(dis)


cor(dis,cophenetic(clus))
cor(dis,cophenetic(cluc))
cor(dis,cophenetic(clua))

#6.2类的显示和解释

plot(cluc)
rect.hclust(cluc,4)
grp<-cutree(cluc,4)

boxplot(A1~grp,data = dune.env,notch=TRUE)
#聚类结果的排序
ord<-cca(dune)
plot(ord,display="site")
ordihull(ord,grp,lty=2,col = "red")

#排序中添加聚类树
plot(ord,display = "site")
ordicluster(ord,clus,col = "blue")
mst<-spantree(dis,toolong = 1)
plot(mst,ord=ord,pch=21,col="red",bg="yellow",type="t")

#6.3分类群落表

wa<-scores(ord,display = "sites",choices = 1)
den<-as.dendrogram(clua)
oden<-reorder(den,wa,mean)
op<-par(mfrow=c(2,1),mar=c(3,5,1,2)+.1)
plot(den)
plot(oden)
par(op)