PCNM（邻体矩阵主坐标分析（principalcoordinates of neighbour matrices）是什么东西呢？

  生态与环境科学很多数据取自地球的表面，因此取样点之间往往具有空间位置关系。研究取样点之间的空间关系属于空间生态学的范畴。如何研究取样点之间的空间关系（或具有什么样的空间结构）？要研究这个问题最基本条件就是得知道取样点位置坐标（如果取样点是经纬度，需要转为笛卡尔坐标系，在R里面的SoDA程序包中geoXY()函数干这个事情）。有了取样点的位置坐标，我们就可以想尽办法通过这个xy坐标来模拟各种可能存在的空间结构。如何获取或是模拟各种各样可能的空间结构，这就是通过主坐标的方法来做这个事情。首先算取样点之间的距离矩阵，然后削减（truncate）距离矩阵规模，只保留一定规模的邻体之间的距离（为什么这么干，就是扭曲一下距离矩阵，可以模拟更多的结构，大家可以不必了解这么多原理）后进行分析主坐标分析，那么得到第1到第n个主坐标都可以说是模拟潜在的空间结构。这些空间结构可以用n个变量构成的取样点-主坐标矩阵来表示，这个就是PCNM变量矩阵。这个PCNM变量矩阵可以作为普通的环境因子矩阵一样，作为解释变量来解释原始的响应变量（可以单变量，也可以是多变量）。我们可以想象一下，如果我们通过各种不同的方法模拟出来的空间结构能够很好解释响应变量（即回归的R2比较大），我们有理由相信响应变量之间存在我们所预先模拟的空间结构。那么为什么响应变量会存在这些模拟的空间结构，有可能是刚好某些环境因子具有这种空间结构，导致了这种空间结构，这些环境因子可能已经被测量，也可能还没有被测量，所有获得PCNM变量后，接下来的工作就是将PCNM作为一组解释变量，跟已经测得的环境变量一起，对响应变量进行解释（也就是变差分解），这样可以看到底是不是已经测量的环境因子是否有类似空间结构来导致响应变量具有类似空间结构，还是未测量的环境变量引起的，当然，也可能说响应变量本身具有的此类空间结构，而不是由环境因子引起部分，这部分统统都是PCNM单独解释的部分。

# 载入本章所用的程序包

library(ape)

library(spdep)

library(vegan)

library(ade4)

# 以下几个程序包可以从https://r-forge.r-project.org/R/?group_id=195下载本地安装

library(packfor)  

library(spacemakeR)

library(AEM)    

library(PCNM)  

# 导入数据

mite <-read.csv("mite.csv",header=T, row.names=1)

mite.env <- read.csv("mite_env.csv",header=T, row.names=1)

mite.xy <- read.csv("mite_xy.csv",header=T, row.names=1)

mite.h <- decostand(mite, "hellinger")

mite.xy.c <-scale(mite.xy, center=TRUE, scale=FALSE) #不对坐标进行标准化

# 1b. ...或自动构建PCNM变量

# library(PCNM) # 如果还未加载PCNM程序包

xy.d1 <- dist(mite.xy)

mite.PCNM.auto <-PCNM(xy.d1)

summary(mite.PCNM.auto)

# PCNM变量Moran指数（由第一距离等级0到削减阈值）；也见PCNM（）函数

# 产生的图（此处无显示图）

# Moran指数的期望值（代表无空间相关）

mite.PCNM.auto$expected_Moran

mite.PCNM.auto$Moran_I

# 正空间相关的特征函数

(select <-which(mite.PCNM.auto$Moran_I$Positive == TRUE))

length(select)  # I > E(I)条件下PCNM变量的数量

mite.PCNM.pos <-as.data.frame(mite.PCNM.auto$vectors)[,select]

# 2.运行基于去趋势甲螨数据的全模型PCNM分析

#------------------------------------------------------------

mite.PCNM.rda <-rda(mite.h.det, mite.PCNM.pos)

anova.cca(mite.PCNM.rda)

# 3.如果分析为显著，计算校正R2

(mite.R2a <-RsquareAdj(mite.PCNM.rda)$adj.r.squared)

# 甲螨-环境-PCNM变差分解

# *****************************

# 1. 环境变量检验和前向选择

# 将环境变量3-5重新编码成二元变量

substrate <-model.matrix(~mite.env[,3])[,-1]

shrubs <-model.matrix(~mite.env[,4])[,-1]

topo <-model.matrix(~mite.env[,5])[,-1]

mite.env2 <-cbind(mite.env[,1:2], substrate, shrubs, topo)

# 环境变量的前向选择

mite.env.rda <-rda(mite.h, mite.env2)

(mite.env.R2a <-RsquareAdj(mite.env.rda)$adj.r.squared)

mite.env.fwd <-forward.sel(mite.h, mite.env2, adjR2thresh=mite.env.R2a,nperm=9999)

env.sign <-sort(mite.env.fwd$order)

env.red <-mite.env2[,c(env.sign)]

colnames(env.red)

# 3. PCNM变量的前向选择

# 运行未去趋势甲螨数据的全模型PCNM分析

mite.undet.PCNM.rda <-rda(mite.h, mite.PCNM.pos)

anova.cca(mite.undet.PCNM.rda)

# 如果分析表明显著，计算校正R2和运行PCNM变量前向选择

(mite.undet.PCNM.R2a <-RsquareAdj(mite.undet.PCNM.rda)$adj.r.squared)

(mite.undet.PCNM.fwd <-forward.sel(mite.h, as.matrix(mite.PCNM.pos),

                                  adjR2thresh=mite.undet.PCNM.R2a))

# 根据R2a准则，如果保留12个PCNM变量，获得的校正R2已经稍大于全模

# 型的校正R2。但这个“稍微超过”也是可行，并不一定很严格。

(nb.sig.PCNM <-nrow(mite.undet.PCNM.fwd)) # 显著的PCNM变量的数量

# 按顺序排列显著的PCNM变量

(PCNM.sign <-sort(mite.undet.PCNM.fwd$order))

# 赋予所有显著PCNM变量一个新的对象

PCNM.red <-mite.PCNM.pos[,c(PCNM.sign)]

# 5.甲螨-环境-趋势-PCNM变差分解

(mite.varpart <-varpart(mite.h, env.red, PCNM.red ))

par(mfrow=c(1,2))

showvarparts(2)

plot(mite.varpart,digits=2)

# 检验单独解释部分[a],  [c] 和 [d]

#********************************

# [a]部分，环境变量单独解释部分

anova.cca(rda(mite.h,env.red, PCNM.red ))

# [b]部分，趋势单独解释部分

anova.cca(rda(mite.h,PCNM.red , env.red))

#仅有环境变量和宽尺度空间变量单独解释部分显著。

PCNM数据和程序包.zip