1）排序圖添加環(huán)境梯度等值線
主要用到vegan包的ordisurf函數(shù)裹虫，主要思路是：
- 將樣點(diǎn)攜帶的環(huán)境因子信息構(gòu)建樣條曲線（splines）
- 樣點(diǎn)的坐標(biāo)信息和樣條曲線的坐標(biāo)信息統(tǒng)一到一個(gè)坐標(biāo)系
> ?ordisurf
Function ordisurf fits a smooth surface for given variable and plots the result on ordination diagram.
ordisurf(x, y, choices = c(1, 2), knots = 10,
family = "gaussian", col = "red", isotropic = TRUE,
thinplate = TRUE, bs = "tp", fx = FALSE, add = FALSE,
display = "sites", w = weights(x, display), main, nlevels = 10,
levels, npoints = 31, labcex = 0.6, bubble = FALSE,
cex = 1, select = TRUE, method = "REML", gamma = 1,
plot = TRUE, lwd.cl = par("lwd"), ...)
# 注意此處的x， R document 的解釋如下：
# For `ordisurf` an ordination configuration, either a matrix or a result known by `scores`
> LatNMDS <- metaMDS(LatOTU.R) # NMDS分析
... ...
... Procrustes: rmse 1.962234e-06 max resid 5.019163e-06
... Similar to previous best
*** Solution reached
# 提取作圖數(shù)據(jù)雳窟，此處筆者主要想展示緯度梯度
> LatNMDS_data <- data.frame(
+ NMDS1 = LatNMDS$points[,1],
+ NMDS2 = LatNMDS$points[,2],
+ MAT= metadata1$MAT,
+ Lat = metadata1$Latitude,
+ Type = metadata1$Type,
+ nam = rownames(LatOTU.R)
+ )
> library(metR)
> ordi <- ordisurf(LatNMDS_data[,1:2], as.numeric(LatNMDS_data$Lat), plot = FALSE, bs = 'ds')
> ordi.grid <- ordi$grid
> ordi.sp <- expand.grid(x = ordi.grid$x, y = ordi.grid$y)
> ordi.sp$z <- as.vector(ordi.grid$z)
> ordi.sp.na <- data.frame(na.omit(ordi.sp))
> ggplot(LatNMDS_data,aes(NMDS1,NMDS2))+
+ geom_point(size=4, aes(shape = Type))+
+ annotate("text",
+ x=min(LatNMDS_data$NMDS1),
+ y=min(LatNMDS_data$NMDS2),
+ hjust=0,
+ vjust=0,
+ label=paste("Stress:",round(LatNMDS$stress,4)))+
+ geom_text(label=LatNMDS_data$nam,
+ size=2,
+ hjust=-0.2,
+ vjust=-0.5)+
+ stat_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..))+
+ geom_text_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..), size = 5)+
+ scale_color_gradient(low = "firebrick3",high = "royalblue3")
結(jié)果如下圖所示：
NMDS+LatGradient
既然函數(shù)ordisurf()的輸入是排序后樣點(diǎn)的坐標(biāo)值，那么其他排序方法也可以添加環(huán)境因子梯度捣作，如PCA分析
# vegan自帶數(shù)據(jù)，意將pH梯度展示
# scores()函數(shù)提取坐標(biāo)值
> data("varespec")
> data("varechem")
> vare.pca <- prcomp(varespec)
> pdata = scores(vare.pca, choices=c(1,2)) %>% as.data.frame()
>
> ordi <- ordisurf(pdata[,1:2], as.numeric(varechem$pH), plot = FALSE, bs = 'ds')
> ordi.grid <- ordi$grid
> ordi.sp <- expand.grid(x = ordi.grid$x, y = ordi.grid$y)
> ordi.sp$z <- as.vector(ordi.grid$z)
> ordi.sp.na <- data.frame(na.omit(ordi.sp))
> ggplot(pdata,aes(PC1,PC2))+
+ geom_point(size=4)+
+ stat_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..))+
+ geom_text_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..), size = 5)+
+ scale_color_gradient(low = "firebrick3",high = "royalblue3")
結(jié)果如圖：
PCA+pHGradient
NMDS分析原理
image
image
image
image
image
image
image
image
image
image
image
image
