博文

数量生态学笔记||非约束排序||PCoA

已有 6137 次阅读 2018-8-14 22:32 |系统分类:科研笔记

主坐标分析（Principal Coordinates Analysis ,PCoA, =Classical Multidimensional Scaling ，cMDScale）是一种探索和可视化数据相似或不同之处的方法。它从一个相似矩阵或相异矩阵（=距离矩阵）开始，并在一个低维度空间中为每个对象分配一个位置。

PCoA能够将样本之间的相似性距离（虚拟距离），经过投影后，在低维度空间进行距离展示，以最大限度地保留原始样本的距离关系，使相似的样本在图形中的距离更为接近，相异的样本距离更远。因此相比于PCA，PCoA以样本距离为整体考虑，更符合生态学数据特征，应用也更为广泛。

PCoA是一种多维标度方法，是一个探索性的过程方法，其目的主要是：

减少（观察）项目
如果可能，在数据中揭示现有结构
揭示相关特征
寻找尽可能低维度的空间（“最小化条件”）
空间必须满足“单调条件”
解释空间的轴，依照假设提供关于感知和评判过程的信息

# 导入本章所需的程序包 
library(ade4)
library(vegan)
library(gclus)  
library(ape)
rm(list = ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA")
# 导入CSV文件数据 
spe <- read.csv("DoubsSpe.csv", row.names=1)
env <- read.csv("DoubsEnv.csv", row.names=1)
spa <- read.csv("DoubsSpa.csv", row.names=1)
# 删除没有数据的样方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]

利用cmdscale 包和vegan包多Doubs数据进行PCoA分析

利用wascores(）函数可以以多度加权平均方式将物种被动投影到样方的PCoA的排序图中。

在非约束排序中被动加入物种，如何解释呢？也同PCA吗？以及后面的约束排序有何区别？

 基于鱼类物种数据Bray-Curtis相异矩阵的PCoA分析
# *********************************************
spe.bray <- vegdist(spe)
spe.b.pcoa <- cmdscale(spe.bray, k=(nrow(spe)-1), eig=TRUE)
# 绘制样方主坐标排序图并用加权平均方法将物种投影到样方PCoA排序图
ordiplot(scores(spe.b.pcoa)[,c(1,2)], type="t", main="PCoA分析（带物种投影）") 
abline(h=0, lty=3)
abline(v=0, lty=3)
# 添加物种
spe.wa <- wascores(spe.b.pcoa$points[,1:2], spe)
text(spe.wa, rownames(spe.wa), cex=0.7, col="red")

# 使用pcoa（）函数运行PCoA分析和物种向量投影
# *****************************************
spe.h.pcoa <- pcoa(dist(spe.h))
# 双序图
par(mfrow=c(1,2))
# 第一个双序图：被动加入Hellinger转化的物种数据
biplot.pcoa(spe.h.pcoa, spe.h, dir.axis2=-1) 
abline(h=0, lty=3)
abline(v=0, lty=3)
# 第二个双序图：被动加入Hellinger转化后标准化的物种数据
spe.std <- apply(spe.h, 2, scale)          
biplot.pcoa(spe.h.pcoa, spe.std, dir.axis2=-1) 
abline(h=0, lty=3)
abline(v=0, lty=3)
#如何比较当前PCoA结果与PCA结果？

> # 基于欧氏和非欧氏距离的PCoA结果比较
> # ***********************************
> # 基于Hellinger距离矩阵PCoA
> is.euclid(dist(spe.h))
[1] TRUE
> #如何比较当前PCoA结果与PCA结果？
> # 基于欧氏和非欧氏距离的PCoA结果比较
> # ***********************************
> # 基于Hellinger距离矩阵PCoA
> is.euclid(dist(spe.h))
[1] TRUE
> summary(spe.h.pcoa)
           Length Class      Mode     
correction   2    -none-     character
note         1    -none-     character
values       5    data.frame list     
vectors    783    -none-     numeric  
trace        1    -none-     numeric  
> spe.h.pcoa$values
    Eigenvalues Relative_eig Broken_stick Cumul_eig Cumul_br_stick
1  7.2228938501 5.133437e-01  0.144128028 0.5133437      0.1441280
2  1.7987448715 1.278400e-01  0.107090991 0.6411837      0.2512190
3  1.2970422885 9.218307e-02  0.088572472 0.7333668      0.3397915
4  1.0780684157 7.662021e-02  0.076226793 0.8099870      0.4160183
5  0.6150272794 4.371107e-02  0.066967534 0.8536980      0.4829858
6  0.4691296076 3.334186e-02  0.059560127 0.8870399      0.5425459
7  0.4122804127 2.930149e-02  0.053387287 0.9163414      0.5959332
8  0.3236126046 2.299972e-02  0.048096282 0.9393411      0.6440295
9  0.1942121662 1.380300e-02  0.043466652 0.9531441      0.6874962
10 0.1685395834 1.197840e-02  0.039351426 0.9651225      0.7268476
11 0.1235468620 8.780692e-03  0.035647722 0.9739032      0.7624953
12 0.0835046563 5.934822e-03  0.032280719 0.9798380      0.7947760
13 0.0759645928 5.398937e-03  0.029194299 0.9852370      0.8239703
14 0.0513764625 3.651415e-03  0.026345296 0.9888884      0.8503156
15 0.0407307694 2.894807e-03  0.023699794 0.9917832      0.8740154
16 0.0313000342 2.224548e-03  0.021230658 0.9940077      0.8952461
17 0.0232647295 1.653465e-03  0.018915843 0.9956612      0.9141619
18 0.0151628431 1.077650e-03  0.016737194 0.9967389      0.9308991
19 0.0133146837 9.462979e-04  0.014679581 0.9976852      0.9455787
20 0.0103052835 7.324146e-04  0.012730263 0.9984176      0.9583090
21 0.0077422970 5.502586e-04  0.010878411 0.9989678      0.9691874
22 0.0063077278 4.483013e-04  0.009114743 0.9994161      0.9783021
23 0.0040023911 2.844570e-04  0.007431241 0.9997006      0.9857334
24 0.0021331293 1.516052e-04  0.005820935 0.9998522      0.9915543
25 0.0013971526 9.929809e-05  0.004277725 0.9999515      0.9958320
26 0.0004271756 3.036012e-05  0.002796244 0.9999819      0.9986283
27 0.0002552902 1.814392e-05  0.001371742 1.0000000      1.0000000
> # 基于Bray-Curtis相异矩阵的PCoA
> is.euclid(spe.bray)
[1] FALSE
> spe.bray.pcoa <- pcoa(spe.bray)
> spe.bray.pcoa$values    # 观察第18轴及之后的特征根
     Eigenvalues  Relative_eig Rel_corr_eig Broken_stick Cum_corr_eig Cumul_br_stick
1   3.695331e+00  5.464785e-01 0.4332316505  0.144128028    0.4332317      0.1441280
2   1.098472e+00  1.624459e-01 0.1343469468  0.107090991    0.5675786      0.2512190
3   7.104740e-01  1.050674e-01 0.0896903953  0.088572472    0.6572690      0.3397915
4   4.149729e-01  6.136766e-02 0.0556797818  0.076226793    0.7129488      0.4160183
5   3.045604e-01  4.503947e-02 0.0429718837  0.066967534    0.7559207      0.4829858
6   1.917884e-01  2.836235e-02 0.0299924247  0.059560127    0.7859131      0.5425459
7   1.569703e-01  2.321333e-02 0.0259850463  0.053387287    0.8118981      0.5959332
8   1.319099e-01  1.950731e-02 0.0231007251  0.048096282    0.8349989      0.6440295
9   1.294251e-01  1.913984e-02 0.0228147364  0.043466652    0.8578136      0.6874962
10  8.667896e-02  1.281839e-02 0.0178948785  0.039351426    0.8757085      0.7268476
11  4.615780e-02  6.825978e-03 0.0132311063  0.035647722    0.8889396      0.7624953
12  3.864487e-02  5.714940e-03 0.0123664087  0.032280719    0.9013060      0.7947760
13  2.745800e-02  4.060586e-03 0.0110788585  0.029194299    0.9123848      0.8239703
14  1.306508e-02  1.932111e-03 0.0094223096  0.026345296    0.9218072      0.8503156
15  7.087896e-03  1.048183e-03 0.0087343669  0.023699794    0.9305415      0.8740154
16  4.039469e-03  5.973709e-04 0.0083835090  0.021230658    0.9389250      0.8952461
17  1.300594e-03  1.923365e-04 0.0080682790  0.018915843    0.9469933      0.9141619
18  0.000000e+00  0.000000e+00 0.0079145195  0.016737194    0.9549078      0.9308991
19 -3.534426e-05 -5.226833e-06 0.0074650365  0.014679581    0.9623729      0.9455787
20 -3.940676e-03 -5.827610e-04 0.0068877933  0.012730263    0.9692607      0.9583090
21 -8.956051e-03 -1.324452e-03 0.0062368816  0.010878411    0.9754975      0.9691874
22 -1.461149e-02 -2.160799e-03 0.0060783322  0.009114743    0.9815759      0.9783021
23 -1.598905e-02 -2.364517e-03 0.0054490158  0.007431241    0.9870249      0.9857334
24 -2.145686e-02 -3.173116e-03 0.0044461644  0.005820935    0.9914711      0.9915543
25 -3.017013e-02 -4.461666e-03 0.0039677639  0.004277725    0.9954388      0.9958320
26 -3.432671e-02 -5.076355e-03 0.0035909666  0.002796244    0.9990298      0.9986283
27 -3.760052e-02 -5.560497e-03 0.0009702192  0.001371742    1.0000000      1.0000000
28 -6.037087e-02 -8.927857e-03 0.0000000000  0.000000000    1.0000000      1.0000000
29 -6.880061e-02 -1.017448e-02 0.0000000000  0.000000000    1.0000000      1.0000000
> # 基于Bray-Curtis相异矩阵平方根的PCoA
> is.euclid(sqrt(spe.bray))
[1] TRUE
> spe.braysq.pcoa <- pcoa(sqrt(spe.bray))
> spe.braysq.pcoa$values  # 观察特征根
   Eigenvalues Relative_eig Broken_stick Cumul_eig Cumul_br_stick
1   3.21560824  0.354677497  0.140256109 0.3546775      0.1402561
2   1.03573822  0.114240607  0.104541823 0.4689181      0.2447979
3   0.80071738  0.088318107  0.086684680 0.5572362      0.3314826
4   0.54426205  0.060031412  0.074779918 0.6172676      0.4062625
5   0.44091188  0.048632019  0.065851347 0.6658996      0.4721139
6   0.38896013  0.042901808  0.058708489 0.7088015      0.5308224
7   0.33391888  0.036830827  0.052756109 0.7456323      0.5835785
8   0.29313251  0.032332143  0.047654068 0.7779644      0.6312325
9   0.25584733  0.028219635  0.043189782 0.8061841      0.6744223
10  0.23652013  0.026087869  0.039221528 0.8322719      0.7136439
11  0.19443183  0.021445584  0.035650099 0.8537175      0.7492940
12  0.18072725  0.019933986  0.032403346 0.8736515      0.7816973
13  0.15112925  0.016669364  0.029427156 0.8903209      0.8111245
14  0.12732111  0.014043356  0.026679903 0.9043642      0.8378044
15  0.10582112  0.011671935  0.024128883 0.9160361      0.8619332
16  0.09697075  0.010695751  0.021747930 0.9267319      0.8836812
17  0.08837033  0.009747135  0.019515787 0.9364790      0.9031970
18  0.07362981  0.008121274  0.017414947 0.9446003      0.9206119
19  0.07259890  0.008007567  0.015430820 0.9526079      0.9360427
20  0.06753369  0.007448880  0.013551121 0.9600568      0.9495938
21  0.06587845  0.007266309  0.011765406 0.9673231      0.9613593
22  0.05253221  0.005794236  0.010064726 0.9731173      0.9714240
23  0.05203557  0.005739458  0.008441350 0.9788568      0.9798653
24  0.04946184  0.005455578  0.006888555 0.9843123      0.9867539
25  0.04130873  0.004556300  0.005400459 0.9888686      0.9921543
26  0.03812350  0.004204973  0.003971888 0.9930736      0.9961262
27  0.03326408  0.003668986  0.002598262 0.9967426      0.9987245
28  0.02953257  0.003257404  0.001275510 1.0000000      1.0000000
> # 基于Bray-Curtis相异矩阵的PCoA（Lingoes校正） 
> spe.brayl.pcoa <- pcoa(spe.bray, correction="lingoes")
> spe.brayl.pcoa$values   # 观察特征根
     Eigenvalues    Corr_eig Rel_corr_eig Broken_stick Cum_corr_eig Cum_br_stick
1   3.695331e+00 3.764131306 0.4332316505  0.144128028    0.4332317    0.1441280
2   1.098472e+00 1.167272861 0.1343469468  0.107090991    0.5675786    0.2512190
3   7.104740e-01 0.779274609 0.0896903953  0.088572472    0.6572690    0.3397915
4   4.149729e-01 0.483773541 0.0556797818  0.076226793    0.7129488    0.4160183
5   3.045604e-01 0.373361024 0.0429718837  0.066967534    0.7559207    0.4829858
6   1.917884e-01 0.260589051 0.0299924247  0.059560127    0.7859131    0.5425459
7   1.569703e-01 0.225770961 0.0259850463  0.053387287    0.8118981    0.5959332
8   1.319099e-01 0.200710550 0.0231007251  0.048096282    0.8349989    0.6440295
9   1.294251e-01 0.198225738 0.0228147364  0.043466652    0.8578136    0.6874962
10  8.667896e-02 0.155479574 0.0178948785  0.039351426    0.8757085    0.7268476
11  4.615780e-02 0.114958410 0.0132311063  0.035647722    0.8889396    0.7624953
12  3.864487e-02 0.107445488 0.0123664087  0.032280719    0.9013060    0.7947760
13  2.745800e-02 0.096258614 0.0110788585  0.029194299    0.9123848    0.8239703
14  1.306508e-02 0.081865696 0.0094223096  0.026345296    0.9218072    0.8503156
15  7.087896e-03 0.075888508 0.0087343669  0.023699794    0.9305415    0.8740154
16  4.039469e-03 0.072840082 0.0083835090  0.021230658    0.9389250    0.8952461
17  1.300594e-03 0.070101207 0.0080682790  0.018915843    0.9469933    0.9141619
18  0.000000e+00 0.068765269 0.0079145195  0.016737194    0.9549078    0.9308991
19 -3.534426e-05 0.064859937 0.0074650365  0.014679581    0.9623729    0.9455787
20 -3.940676e-03 0.059844562 0.0068877933  0.012730263    0.9692607    0.9583090
21 -8.956051e-03 0.054189118 0.0062368816  0.010878411    0.9754975    0.9691874
22 -1.461149e-02 0.052811562 0.0060783322  0.009114743    0.9815759    0.9783021
23 -1.598905e-02 0.047343750 0.0054490158  0.007431241    0.9870249    0.9857334
24 -2.145686e-02 0.038630480 0.0044461644  0.005820935    0.9914711    0.9915543
25 -3.017013e-02 0.034473899 0.0039677639  0.004277725    0.9954388    0.9958320
26 -3.432671e-02 0.031200098 0.0035909666  0.002796244    0.9990298    0.9986283
27 -3.760052e-02 0.008429746 0.0009702192  0.001371742    1.0000000    1.0000000
28 -6.037087e-02 0.000000000 0.0000000000  0.000000000    1.0000000    1.0000000
29 -6.880061e-02 0.000000000 0.0000000000  0.000000000    1.0000000    1.0000000
> # 基于Bray-Curtis相异矩阵的PCoA（Cailliez校正）
> spe.brayc.pcoa <- pcoa(spe.bray, correction="cailliez")
> spe.brayc.pcoa$values   # 观察特征根
     Eigenvalues   Corr_eig Rel_corr_eig Broken_stick Cum_corr_eig Cum_br_stick
1   3.695331e+00 5.20461437 0.4442681027  0.144128028    0.4442681    0.1441280
2   1.098472e+00 1.60465006 0.1369736134  0.107090991    0.5812417    0.2512190
3   7.104740e-01 1.09152082 0.0931726832  0.088572472    0.6744144    0.3397915
4   4.149729e-01 0.68985417 0.0588862466  0.076226793    0.7333006    0.4160183
5   3.045604e-01 0.52129425 0.0444978992  0.066967534    0.7777985    0.4829858
6   1.917884e-01 0.38710929 0.0330438139  0.059560127    0.8108424    0.5425459
7   1.569703e-01 0.36447367 0.0311116277  0.053387287    0.8419540    0.5959332
8   1.319099e-01 0.29671255 0.0253275095  0.048096282    0.8672815    0.6440295
9   1.294251e-01 0.27559544 0.0235249448  0.043466652    0.8908064    0.6874962
10  8.667896e-02 0.21600584 0.0184383509  0.039351426    0.9092448    0.7268476
11  4.615780e-02 0.15419396 0.0131620626  0.035647722    0.9224069    0.7624953
12  3.864487e-02 0.15378333 0.0131270104  0.032280719    0.9355339    0.7947760
13  2.745800e-02 0.11812808 0.0100834637  0.029194299    0.9456173    0.8239703
14  1.306508e-02 0.08848541 0.0075531526  0.026345296    0.9531705    0.8503156
15  7.087896e-03 0.07304055 0.0062347726  0.023699794    0.9594053    0.8740154
16  4.039469e-03 0.06999353 0.0059746780  0.021230658    0.9653799    0.8952461
17  1.300594e-03 0.05712927 0.0048765792  0.018915843    0.9702565    0.9141619
18  0.000000e+00 0.05587583 0.0047695843  0.016737194    0.9750261    0.9308991
19 -3.534426e-05 0.05432215 0.0046369623  0.014679581    0.9796631    0.9455787
20 -3.940676e-03 0.04912221 0.0041930931  0.012730263    0.9838562    0.9583090
21 -8.956051e-03 0.04100207 0.0034999542  0.010878411    0.9873561    0.9691874
22 -1.461149e-02 0.03777775 0.0032247250  0.009114743    0.9905808    0.9783021
23 -1.598905e-02 0.03451234 0.0029459878  0.007431241    0.9935268    0.9857334
24 -2.145686e-02 0.02959507 0.0025262474  0.005820935    0.9960531    0.9915543
25 -3.017013e-02 0.02436729 0.0020800022  0.004277725    0.9981331    0.9958320
26 -3.432671e-02 0.01902747 0.0016241925  0.002796244    0.9997573    0.9986283
27 -3.760052e-02 0.00284371 0.0002427403  0.001371742    1.0000000    1.0000000
28 -6.037087e-02 0.00000000 0.0000000000  0.000000000    1.0000000    1.0000000
29 -6.880061e-02 0.00000000 0.0000000000  0.000000000    1.0000000    1.0000000
#如果要选择承载最大比例变差的前两轴去了解数据的结构，你会选择上面哪
#种结果呢？