作者（Jari Okanen）面向熟悉 R 操作的读者。想要运行文中实例、得到对应的结果，有时要先加载一些东西（原文 session info 部分，在6.3小节末）。

chapter2.1

library(permute)
library(vegan)
data(varespec)
vare.dis <- vegdist(varespec)
vare.mds0 <- monoMDS(vare.dis)

stressplot(vare.mds0, vare.dis)

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

vare.mds <- metaMDS(varespec, trace = FALSE)
vare.mds

plot(vare.mds, type = "t")

 

chapter 2.2

library(permute)
library(vegan)
data(varechem)
rankindex(scale(varechem), varespec, c("euc","man","bray","jac","kul"))

dis <- vegdist(decostand(varespec, "norm"), "euclid")

dis <- vegdist(decostand(varespec, "hell"), "euclidean")

d <- vegdist(varespec, "bray", binary = TRUE)
d <- designdist(varespec, "(A+B-2*J)/(A+B)")
d <- designdist(varespec, "(b+c)/(2*a+b+c)", abcd=TRUE)

 

chapter2.3

tmp <- wisconsin(sqrt(varespec))
dis <- vegdist(tmp)
vare.mds0 <- monoMDS(dis)
pro <- procrustes(vare.mds, vare.mds0)
pro

plot(pro)

plot(pro, kind = 2)

 

chapter 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

plot(vare.pca, scaling = 3)

dim(varespec)

vare.ca <- cca(varespec)
vare.ca
plot(vare.ca)

chisq.test(varespec/sum(varespec))

plot(vare.ca, scaling = 1)

 

chapter 2.5

vare.dca <- decorana(varespec)
vare.dca

plot(vare.dca, display="sites")

 

chapter 2.6

data(BCI)
mod <- decorana(BCI)
plot(mod)

names(BCI)[1:5]
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="+")

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

sel[1:14]

 

chapter 3.1

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

 

chapter 3.2

ef <- envfit(vare.mds ~ Al + Ca, varechem)

plot(vare.mds, display = "sites")

plot(ef)

tmp <- with(varechem, ordisurf(vare.mds, Al, add = TRUE))

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

 

chapter 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))
with(dune.env, ordispider(dune.ca, Management, col = "blue", label= TRUE))
with(dune.env, ordihull(dune.ca, Management, col="blue", lty=2))

 

chapter 4.1

library(permute)
library(vegan)
data(varechem)

vare.cca <- cca(varespec ~ Al + P + K, varechem)
vare.cca

plot(vare.cca)

library(scatterplot3d)
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)

 

chapter 4.2

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)

 

chapter 4.3

mod1 <- cca(varespec ~ ., varechem)
mod1

plot(procrustes(cca(varespec), mod1))

mod0 <- cca(varespec ~ 1, varechem)
mod <- step(mod0, scope = formula(mod1), test = "perm")
mod

modb <- step(mod1, scope = list(lower = formula(mod0), upper = formula(mod1)), trace = 0)
modb

modb$anova

vif.cca(mod1)

vif.cca(mod)

 

chapter 4.4

spenvcor(mod)

dune.cca <- cca(dune ~ Management, dune.env)
plot(dune.cca, display = c("lc", "wa"), type = "p")
ordispider(dune.cca, col="blue")

chapter 4.5

pred <- calibrate(mod)
head(pred)

with(varechem, plot(Al, pred[,"Al"] - Al, ylab="Prediction Error"))
abline(h=0, col="grey")

library(mgcv)
plot(mod, display = c("bp", "wa", "lc"))
ef <- with(varechem, ordisurf(mod, Al, display = "lc", add = TRUE))

chapter 4.6

dune.cca <- cca(dune ~ Management + Condition(Moisture), dune.env)
plot(dune.cca)
dune.cca

anova(dune.cca, perm.max = 2000)
with(dune.env, anova(dune.cca, strata = Moisture))

 

chapter 5.1

library(permute)
library(vegan)
betad <- betadiver(dune, "z")
adonis(betad ~ Management, dune.env, perm=200)

adonis(betad ~ A1*Management, dune.env, perm = 200)

chapter 5.2

mod <- with(dune.env, betadisper(betad, Management))
mod

plot(mod)
boxplot(mod)

anova(mod)
permutest(mod)

chapter 5.3

pc <- prcomp(varechem, scale = TRUE)
pc<- scores(pc, display = "sites", choices = 1:4)
edis <- vegdist(pc, method = "euclid")
vare.dis <- vegdist(wisconsin(sqrt(varespec)))
mantel(vare.dis, edis)

plot(vare.dis, edis)

chapter 5.4

pc <- scores(pc, choices = 1:2)
pro <- protest(vare.mds, pc)
plot(pro)
pro

 

chapter 6.1

library(permute)
library(vegan)
data（dune）

dis <- vegdist(dune)
clus <- hclust(dis, "single")
plot(clus)

cluc <- hclust(dis, "complete")
plot(cluc)

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))

chapter 6.2

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

boxplot(A1 ~ grp, data=dune.env, notch = TRUE)

ord <- cca(dune)
plot(ord, display = "sites")
ordihull(ord, grp, lty = 2, col = "red")

plot(ord, display="sites")
ordicluster(ord, cluc, col="blue")

mst <- spantree(dis, toolong = 1)
plot(mst, ord=ord, pch=21, col = "red", bg = "yellow", type = "t")

chapter 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)

vegemite(dune, use = oden, zero = "-")