网站开发用什么技术可行性,平面设计找素材的网站,百度关键词搜索量查询,桂林做网站的公司上一节数量生态学笔记||冗余分析(RDA)概述中#xff0c;我们回顾了RDA的计算过程#xff0c;不管这个过程我们有没有理解透彻#xff0c;我希望你能知道的是#xff1a;RDA是响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析。本节我们就是具体来看一个RD…上一节数量生态学笔记||冗余分析(RDA)概述中我们回顾了RDA的计算过程不管这个过程我们有没有理解透彻我希望你能知道的是RDA是响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析。本节我们就是具体来看一个RDA的分析案例来看看里面的参数以及结果的解读。# CHAPTER 6 - CANONICAL ORDINATION# ********************************# 载入所需程序包library(ade4)library(vegan)library(packfor)# 可在http://r-forge.r-project.org/R/?group_id195下载但是好像在R 3.5.1上加载不了所以这篇我用R3.4来做的。packfor已经不用函数都搬到adespatial# 如果是MacOS X系统packfor程序包内forward.sel函数的运行需要加载# gfortran程序包。用户必须从cran.r-project.org网站内选择MacOS X# 然后选择tools安装gfortran程序包。rm(list ls())setwd(D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA)library(MASS)library(ellipse)library(FactoMineR)# 附加函数source(evplot.R)source(hcoplot.R)# 导入CSV数据文件spe env spa # 删除没有数据的样方8spe env spa # 提取环境变量das(离源头距离)以备用das # 从环境变量矩阵剔除das变量env # 将slope变量(pen)转化为因子(定性)变量pen2 pen2[env$pen quantile(env$pen)[4]] steeppen2[env$pen quantile(env$pen)[3]] moderatepen2[env$pen quantile(env$pen)[2]] lowpen2 table(pen2)# 生成一个含定性坡度变量的环境变量数据框env2env2 env2$pen # 将所有解释变量分为两个解释变量子集# 地形变量(上下游梯度)子集envtopo names(envtopo)#水体化学属性变量子集envchem names(envchem)# 物种数据Hellinger转化spe.hel 使用vegan包运行RDAvegan包运行RDA有两种不同的模式。第一种是简单模式直接输入用逗号隔开的数据矩阵对象到rda()函数式中为响应变量矩阵为解释变量矩阵为偏RDA分析需要的协变量矩阵。此公式有一个缺点不能有因子变量(定性变量)。如果有因子变量建议使用第二种模式式中为响应变量矩阵。解释变量矩阵包括定量变量(var1)、因子变量(factorA)以及变量2和变量3的交互作用项协变量(var4)被放到Condition()里。所用的数据都放在XWdata的数据框里。这个公式与lm()函数以及其他回归函数一样左边是响应变量右边是解释变量。# 基于Hellinger 转化的鱼类数据RDA解释变量为对象env2包括的环境变量# 关注省略模式的公式spe.rda summary(spe.rda) # 2型标尺(默认)#这里使用一些默认的选项即 scaleFALSE(基于协方差矩阵的RDA)和#scaling2RDA结果的摘录RDA formula :Call:rda(formula spe.hel ~ alt pen deb pH dur pho nit amm oxy dbo, data env2)方差分解(Partitioning of variance)总方差被划分为约束和非约束两部分。约束部分表示响应变量矩阵的总方差能被解释变量解释的部分如果用比例来表示其值相当于多元回归的。在RDA中这个解释比例值也称作双多元冗余统计。然而类似多元回归的未校正的RDA的是有偏差的需要进一步校正。Partitioning of variance:Inertia ProportionTotal 0.5025 1.0000Constrained 0.3654 0.7271Unconstrained 0.1371 0.2729特征根以及对方差的贡献率(Eigenvalues, and their contribution to the variance ):当前这个RDA分析产生了12个典范轴(特征根用RDA1 至RDA12表示)和16个非约束轴(特征根用PC1至PC16表示)。输出结果不仅包含每轴特征根同时也给出累积方差解释率(约束轴)或承载轴(非约束轴)最终的累计值必定是1.12 个典范轴累积解释率也代表响应变量总方差能够被解释变量解释的部分。两个特征根的重要区别典范特征根RDAx是响应变量总方差能够被解释变量解释的部分而残差特征根RCx响应变量总方差能够被残差轴解释的部分与RDA无关。Eigenvalues, and their contribution to the varianceImportance of components:RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11 RDA12 PC1 PC2 PC3Eigenvalue 0.2281 0.0537 0.03212 0.02321 0.008707 0.007218 0.004862 0.002919 0.002141 0.001160 0.0009134 0.0003406 0.04580 0.02814 0.01529Proportion Explained 0.4539 0.1069 0.06393 0.04618 0.017328 0.014363 0.009676 0.005809 0.004260 0.002308 0.0018176 0.0006778 0.09115 0.05600 0.03042Cumulative Proportion 0.4539 0.5607 0.62467 0.67085 0.688176 0.702539 0.712215 0.718025 0.722284 0.724592 0.7264100 0.7270878 0.81824 0.87424 0.90466PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16Eigenvalue 0.01399 0.009841 0.007676 0.004206 0.003308 0.002761 0.002016 0.001752 0.0009851 0.0005921 0.0004674 0.0002127 0.0001004Proportion Explained 0.02784 0.019584 0.015276 0.008371 0.006583 0.005495 0.004013 0.003486 0.0019604 0.0011783 0.0009301 0.0004233 0.0001998Cumulative Proportion 0.93250 0.952084 0.967360 0.975731 0.982314 0.987809 0.991822 0.995308 0.9972684 0.9984468 0.9993768 0.9998002 1.0000000累积约束特征根(Accumulated constrained eigenvalues)表示在本轴以及前面所有轴的典范轴所能解释的方差占全部解释方差的比例累积。Accumulated constrained eigenvaluesImportance of components:RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11 RDA12Eigenvalue 0.2281 0.0537 0.03212 0.02321 0.008707 0.007218 0.004862 0.002919 0.002141 0.001160 0.0009134 0.0003406Proportion Explained 0.6242 0.1470 0.08792 0.06352 0.023832 0.019755 0.013308 0.007990 0.005859 0.003174 0.0024999 0.0009322Cumulative Proportion 0.6242 0.7712 0.85913 0.92265 0.946483 0.966237 0.979545 0.987535 0.993394 0.996568 0.9990678 1.0000000物种得分(Species scores)双序图和三序图内代表响应变量的箭头的顶点坐标。与PCA相同坐标依赖标尺Scaling的选择。Scaling 2 for species and site scores* Species are scaled proportional to eigenvalues* Sites are unscaled: weighted dispersion equal on all dimensions* General scaling constant of scores: 1.93676Species scoresRDA1 RDA2 RDA3 RDA4 RDA5 RDA6CHA 0.13383 0.11623 -0.238180 0.018611 0.043221 -0.029737TRU 0.64238 0.06648 0.123713 0.181572 -0.009691 0.029793VAI 0.47475 0.07015 -0.010218 -0.115369 -0.045317 -0.030033LOC 0.36260 0.06972 0.041240 -0.190586 -0.046881 0.006448OMB 0.13079 0.10709 -0.239224 0.043603 0.065881 0.003458BLA 0.06587 0.12475 -0.216900 -0.004157 0.021793 -0.004195HOT -0.17417 0.06778 -0.008426 -0.016419 -0.079730 0.044706TOX -0.12683 0.16052 -0.035733 -0.016087 -0.089768 -0.001880VAN -0.07963 0.04200 0.007636 -0.059179 -0.033596 -0.121440CHE -0.10903 -0.17552 -0.090099 -0.168373 0.019444 0.008745BAR -0.18528 0.21154 -0.073087 -0.006879 -0.012995 0.040484SPI -0.16064 0.15513 -0.014309 -0.002488 -0.060810 0.011045GOU -0.20537 0.02484 -0.007973 -0.017742 -0.049137 -0.096231BRO -0.10734 0.02848 0.090055 0.012324 0.075184 -0.057088PER -0.09164 0.10506 0.070393 -0.057443 0.013870 -0.009906BOU -0.20907 0.16002 0.025500 0.012078 -0.011477 0.022035PSO -0.22799 0.11121 0.018800 -0.009474 -0.027431 0.024517ROT -0.16098 0.01348 0.041628 0.032398 0.054117 -0.094582CAR -0.17384 0.14901 0.022262 0.009534 0.004991 -0.005396TAN -0.14025 0.10632 0.078290 -0.122627 0.054162 0.031256BCO -0.18594 0.12222 0.053881 0.026170 0.044015 0.014577PCH -0.14630 0.08894 0.061880 0.034763 0.083530 0.004396GRE -0.30881 0.01606 0.039366 0.029254 -0.011141 -0.052412GAR -0.31982 -0.16601 -0.018225 -0.115454 0.054341 0.064772BBO -0.23897 0.09090 0.051627 0.010224 0.007004 0.036497ABL -0.43215 -0.22639 -0.108190 0.138807 -0.083920 0.008460ANG -0.19442 0.14149 0.033659 0.017387 0.008110 0.017638样方得分(Site scores (weighted sums of species scores))物种得分的加权和使用响应变量矩阵计算获得的样方坐标。Site scores (weighted sums of species scores)RDA1 RDA2 RDA3 RDA4 RDA5 RDA61 0.40151 -0.154306 0.55539 1.600773 0.191866 0.9168932 0.53523 -0.025084 0.43389 0.294615 -0.518456 0.4588603 0.49430 -0.014605 0.49409 0.169038 -0.246166 0.1634214 0.33452 0.001173 0.51626 -0.321009 0.088716 -0.2198375 0.02794 -0.194357 0.44612 -0.559210 0.853768 -1.1156546 0.24422 -0.130778 0.41372 -0.696264 0.181514 -0.2734737 0.46590 -0.125982 0.31674 -0.137834 -0.548635 -0.0617039 0.03662 -0.605060 -0.07022 -1.260916 0.669108 1.16498610 0.31381 -0.198654 0.10764 -0.635139 -0.741448 -0.99023611 0.48116 -0.039598 -0.37851 0.181924 0.221494 0.25451112 0.49162 0.014263 -0.37983 0.163103 0.223730 0.32467213 0.49848 0.212367 -0.67408 0.518823 0.400091 0.22162214 0.38202 0.229538 -0.75771 0.223651 0.515712 -0.13974015 0.28739 0.218713 -0.71887 -0.210821 0.176392 -0.55318516 0.09129 0.400192 -0.34443 -0.376097 -0.366868 -0.57523017 -0.05306 0.423994 -0.41009 -0.188492 -0.726152 0.15187618 -0.14185 0.385926 -0.36814 -0.217143 -0.644298 -0.00105219 -0.28204 0.275528 -0.01877 -0.371457 -0.691725 -0.06223020 -0.39683 0.209468 0.11547 -0.177972 -0.387121 0.04869021 -0.42851 0.278256 0.22010 -0.005993 -0.027083 -0.04220922 -0.46553 0.251819 0.22784 0.040192 0.152965 0.03218523 -0.28123 -1.145599 -0.50543 0.300015 -0.004403 1.15720624 -0.40893 -0.752909 -0.26785 0.428851 -0.645606 0.64308425 -0.35205 -0.770380 -0.12186 0.459170 0.078892 -1.72597326 -0.46923 0.093958 0.23058 0.107865 0.310432 0.13255627 -0.47021 0.213521 0.24887 0.084219 0.331685 0.12543928 -0.47266 0.233922 0.27053 0.105776 0.381436 0.09371929 -0.37457 0.393260 0.10423 0.202115 0.282621 0.02183430 -0.48932 0.321417 0.31431 0.278218 0.487541 -0.151031样方约束——解释变量的线性组合(Site constraints (linear combinations of constraining variables))使用解释变量矩阵计算获得的样方坐标是拟合的(fitted)样方坐标。Site constraints (linear combinations of constraining variables)RDA1 RDA2 RDA3 RDA4 RDA5 RDA61 0.55135 0.002395 0.47774 0.626878 -0.210700 0.315112 0.29737 0.105715 0.64862 0.261161 -0.057741 0.093223 0.36834 -0.185376 0.59788 0.324212 -0.002385 0.310904 0.44348 -0.066086 0.33260 -0.344463 -0.279591 -0.370795 0.27004 -0.366721 0.17992 -0.453274 0.716614 -0.065456 0.37148 -0.255624 0.40847 0.217259 0.023374 0.343607 0.53874 -0.179999 0.06845 -0.424896 0.024884 -0.334549 -0.04438 -0.362632 0.12371 -1.180662 0.348188 0.2635210 0.16289 -0.154212 0.22252 -0.241425 -0.573048 -0.0286711 0.29912 0.176150 -0.08233 0.003924 0.164806 -0.4460312 0.36843 0.197492 -0.41095 0.300566 -0.053922 -0.0113913 0.42626 0.190107 -0.59764 0.100988 0.118714 -0.2102214 0.34373 0.134362 -0.80378 0.063879 0.665797 0.4812615 0.21385 0.237182 -0.56341 -0.001099 -0.028564 0.0165516 0.03056 0.352192 -0.12110 -0.202316 0.058413 -0.4354217 -0.10499 0.178587 -0.26925 0.046988 -0.608314 0.2123718 -0.11204 0.221631 -0.24024 -0.302957 -0.251346 -0.0144819 -0.05479 0.311860 -0.30701 0.010366 -0.481829 -0.1285520 -0.25684 0.303770 -0.06768 -0.036587 -0.562578 0.1369821 -0.39177 0.196355 0.01877 -0.281086 -0.383524 0.3931022 -0.21361 0.180414 0.01066 0.074301 -0.036849 -0.0242923 -0.21654 -1.016853 -0.57298 0.548175 0.182594 0.5144324 -0.52578 -0.645438 -0.11182 -0.240149 -0.512492 0.3242025 -0.38886 -0.867381 -0.08079 0.482839 -0.106743 -1.2130526 -0.48440 0.031510 0.14065 -0.114545 0.425712 -0.1798927 -0.61221 0.138191 0.32316 -0.015795 0.232397 0.3828828 -0.46921 0.459843 0.22002 0.078870 0.278747 -0.3650429 -0.38450 0.344487 0.20622 0.353008 0.504693 0.1774730 -0.42572 0.338080 0.24960 0.345839 0.404692 -0.13777解释变量双序图得分(Biplot scores for constraining variables)排序图内解释解释变量箭头的坐标按照下面的过程获得运行解释变量与拟合的样方坐标之间的相关分析然后将所有相关系数转化为双序图内坐标。所有的变量包括个水平的因子口可以有自己的坐标对因子变量在排序轴的坐标用各个因子的形心表示更合适。Biplot scores for constraining variablesRDA1 RDA2 RDA3 RDA4 RDA5 RDA6alt 0.8239 -0.203257 0.46604 -0.16936 0.003229 0.10407penmoderate -0.3592 -0.008729 -0.21727 -0.18278 0.157934 0.50094pensteep 0.2791 0.156027 -0.06876 0.01927 0.176390 -0.15469penvery_steep 0.6129 -0.148496 0.45379 0.03618 -0.191046 -0.04715deb -0.7770 0.254952 -0.17470 0.30995 0.227580 -0.11938pH 0.1023 0.178431 -0.30131 0.03959 0.298584 0.04848dur -0.5722 0.044963 -0.56040 -0.14813 0.275617 -0.24527pho -0.4930 -0.650488 -0.19868 0.29286 0.055893 -0.39345nit -0.7755 -0.203836 -0.23285 0.19744 -0.078849 -0.35073amm -0.4744 -0.687577 -0.16648 0.28470 -0.051768 -0.33852oxy 0.7632 0.575528 -0.16425 0.08026 -0.136143 0.13748dbo -0.5171 -0.791805 -0.15652 0.22064 0.075568 -0.09105因子解释变量形心(Centroids for factor constraints)因子变量各个水平形心点的坐标即每个水平所用标识为一的样方的形心。Centroids for factor constraintsRDA1 RDA2 RDA3 RDA4 RDA5 RDA6penlow -0.2800 0.005530 -0.09025 0.07614 -0.07860 -0.18390penmoderate -0.2093 -0.005086 -0.12660 -0.10650 0.09203 0.29189pensteep 0.1965 0.109867 -0.04842 0.01357 0.12420 -0.10892penvery_steep 0.3908 -0.094679 0.28933 0.02307 -0.12181 -0.03006在rda()函数中大家感兴趣的典范特征系数(即每个解释变量与每个典范轴之间的回归系数)可用coef()函数获得#如何从rda()输出结果中获得典范系数coef(spe.rda)alt 0.0004482548 7.559499e-05 0.0005205825 0.0003883845 0.001857206 -6.313946e-05 -0.001355362 0.001120849 -0.0002530083 0.001189659penmoderate -0.0123961693 -1.660194e-02 0.0160069104 -0.0278562534 0.276128846 1.310695e-01 -0.022918427 0.018830063 -0.3113354204 -0.278006278pensteep 0.0478390238 4.893302e-02 0.1022577908 0.1347997771 0.393812929 -1.795824e-01 0.046319741 0.123642821 0.0963820046 -0.447099975penvery_steep 0.0180005587 -5.691933e-02 0.2322637617 0.1002359565 0.036814635 -1.742060e-01 0.517299284 0.067564014 -0.2262450630 -0.590022907deb -0.0014063766 4.440084e-03 0.0089298486 0.0164715901 0.013318449 2.705757e-03 -0.002419805 0.010394632 -0.0006430624 0.004427139pH 0.0188227657 -3.163592e-02 -0.0482021704 0.1141647498 0.412886847 1.091403e-01 0.139806409 -0.436295510 -0.0215841003 -0.904063764dur 0.0025580344 -1.955496e-03 -0.0065901935 -0.0093556696 0.005228707 -6.098098e-03 0.002195518 0.010536248 -0.0006844877 0.003353110pho 0.1031541920 4.583584e-02 -0.1000153096 -0.1050243435 0.422991234 -3.694132e-01 0.035874664 -0.701043138 -0.2865315085 0.245917738nit -0.0123824749 1.041485e-01 0.0625396719 0.0774218297 0.234401221 -3.541252e-02 -0.240428544 0.128403162 -0.0686364968 0.113090041amm -0.1084411839 -4.407786e-01 0.0057247742 0.0538542716 -1.812468883 4.798631e-03 0.281937862 1.068013480 0.3084215704 -1.217501183oxy 0.0686692124 1.980446e-02 -0.0894153251 0.1200884061 0.032052566 3.880445e-02 -0.058026043 0.061374900 -0.0196444146 0.089881042dbo 0.0108322463 -2.696114e-02 -0.0253225230 0.0745175780 0.067082880 9.276548e-02 -0.019719504 0.047865971 0.0359365102 0.065416035RDA11 RDA12alt 0.0006826822 0.0009471677penmoderate 0.0398080898 -0.2974896027pensteep 0.2445035939 -0.3475535793penvery_steep 0.2457103975 -0.1848717482deb -0.0022565029 0.0064858596pH 0.0696045266 0.5756301035dur 0.0007758175 0.0062181193pho -0.0015544897 -0.6309008260nit 0.3983543655 0.0942246573amm -1.5964107514 0.8979015239oxy 0.0627415675 0.0258937429dbo 0.1113928484 0.0403158432提取。解读和绘制vegan包输出的RDA结果校正# 提取校正R2# **********# 从rda的结果中提取未校正R2(R2 [1] 0.7270878# 从rda的结果中提取校正R2(R2adj [1] 0.5224036#可以看出校正R2总是小于R2。校正R2作为被解释方差比例的无偏估计后#面的变差分解部分所用的也是校正R2。# RDA三序图现在绘制RDA的排序图。如果一张排序图中有三个实体样方、响应变量、解释变量这种排序图称为三序图(triplot)为了区分响应变量和解释变量定量解释变量用箭头表示响应变量用不带箭头的线表示。# RDA三序图# **********# 1型标尺距离三序图plot(spe.rda, scaling1, mainRDA三序图spe.helenv2 - 1型标尺- 加权和样方坐标)#此排序图同时显示所有的元素样方、物种、定量解释变量(用箭头表示)#和因子变量的形心。为了与定量解释变量区分物种用不带箭头的线表示。spe.sc arrows(0, 0, spe.sc[, 1], spe.sc[, 2], length0, lty1, colred)plot(spe.rda, mainRDA三序图spe.helenv2 - 2型标尺- 加权和样方坐标)spe2.sc arrows(0, 0, spe2.sc[, 1], spe2.sc[, 2], length0, lty1, colred)# 样方坐标是环境因子线性组合# 1型标尺plot(spe.rda, scaling1, displayc(sp, lc, cn),mainRDA三序图spe.helenv2 - 1型标尺- 拟合的样方坐标)arrows(0, 0, spe.sc[, 1], spe.sc[, 2], length0, lty1, colred)# 2型标尺plot(spe.rda, displayc(sp, lc, cn),mainRDA三序图spe.helenv2 - 2型标尺- 拟合的样方坐标)arrows(0, 0, spe2.sc[,1], spe2.sc[,2], length0, lty1, colred)RDA 结果的置换检验# RDA所有轴置换检验anova.cca(spe.rda, step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(formula spe.hel ~ alt pen deb pH dur pho nit amm oxy dbo, data env2)Df Variance F Pr(F)Model 12 0.36537 3.5522 0.001 ***Residual 16 0.13714---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1# 每个典范轴逐一检验anova.cca(spe.rda, byaxis, step1000)Permutation test for rda under reduced modelForward tests for axesPermutation: freeNumber of permutations: 999Model: rda(formula spe.hel ~ alt pen deb pH dur pho nit amm oxy dbo, data env2)Df Variance F Pr(F)RDA1 1 0.228081 26.6098 0.001 ***RDA2 1 0.053696 6.2646 0.003 **RDA3 1 0.032124 3.7478 0.361RDA4 1 0.023207 2.7075 0.763RDA5 1 0.008707 1.0159 1.000RDA6 1 0.007218 0.8421 1.000RDA7 1 0.004862 0.5673 1.000RDA8 1 0.002919 0.3406 1.000RDA9 1 0.002141 0.2497 1.000RDA10 1 0.001160 0.1353 1.000RDA11 1 0.000913 0.1066 1.000RDA12 1 0.000341 0.0397 1.000Residual 16 0.137141---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1#每个典范轴的检验只能输入由公式模式获得的rda结果。有多少个轴结果是#显著的# 使用Kaiser-Guttman准则确定残差轴spe.rda$CA$eig[spe.rda$CA$eig mean(spe.rda$CA$eig)]PC1 PC2 PC3 PC4 PC50.045802781 0.028143080 0.015288209 0.013987518 0.009841239#很明显还有一部分有意思的变差尚未被目前所用的这套环境变量解释。偏RDA分析偏RDA固定地形变量影响后水体化学属性的效应# 简单模式X和W可以是分离的定量变量表格spechem.physio spechem.physioCall: rda(X spe.hel, Y envchem, Z envtopo)Inertia Proportion RankTotal 0.5025 1.0000Conditional 0.2087 0.4153 3Constrained 0.1602 0.3189 7Unconstrained 0.1336 0.2659 18Inertia is varianceEigenvalues for constrained axes:RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA70.09137 0.04591 0.00928 0.00625 0.00387 0.00214 0.00142Eigenvalues for unconstrained axes:PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC80.04642 0.02072 0.01746 0.01325 0.00974 0.00588 0.00512 0.00400(Showed only 8 of all 18 unconstrained eigenvalues)summary(spechem.physio)# 公式模式X和W必须在同一数据框内class(env)spechem.physio2 Condition(alt pen deb), dataenv)spechem.physio2Call: rda(formula spe.hel ~ pH dur pho nit amm oxy dbo Condition(alt pen deb), data env)Inertia Proportion RankTotal 0.5025 1.0000Conditional 0.2087 0.4153 3Constrained 0.1602 0.3189 7Unconstrained 0.1336 0.2659 18Inertia is varianceEigenvalues for constrained axes:RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA70.09137 0.04591 0.00928 0.00625 0.00387 0.00214 0.00142Eigenvalues for unconstrained axes:PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC80.04642 0.02072 0.01746 0.01325 0.00974 0.00588 0.00512 0.00400(Showed only 8 of all 18 unconstrained eigenvalues)#上面两个分析的结果完全相同。偏RDA检验(使用公式模式获得的RDA结果以便检验每个轴)anova.cca(spechem.physio2, step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(formula spe.hel ~ pH dur pho nit amm oxy dbo Condition(alt pen deb), data env)Df Variance F Pr(F)Model 7 0.16024 3.0842 0.001 ***Residual 18 0.13360---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1anova.cca(spechem.physio2, step1000, byaxis)Permutation test for rda under reduced modelForward tests for axesPermutation: freeNumber of permutations: 999Model: rda(formula spe.hel ~ pH dur pho nit amm oxy dbo Condition(alt pen deb), data env)Df Variance F Pr(F)RDA1 1 0.091373 12.3108 0.001 ***RDA2 1 0.045907 6.1851 0.010 **RDA3 1 0.009277 1.2499 0.964RDA4 1 0.006251 0.8422 0.993RDA5 1 0.003866 0.5209 0.996RDA6 1 0.002142 0.2886 1.000RDA7 1 0.001425 0.1920 0.997Residual 18 0.133599---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1# 偏RDA三序图(使用拟合值的样方坐标)# 1型标尺plot(spechem.physio, scaling1, displayc(sp, lc, cn),mainRDA三序图spe.helchem ︳Tope- 1型标尺- 拟合的样方坐标)spe3.sc arrows(0, 0, spe3.sc[, 1], spe3.sc[, 2], length0, lty1, colred)# 2型标尺plot(spechem.physio, displayc(sp, lc, cn),mainRDA三序图spe.helchem ︳Tope- 2型标尺- 拟合的样方坐标)spe4.sc arrows(0, 0, spe4.sc[,1], spe4.sc[,2], length0, lty1, colred)解释变量向前选择每个变量的共线性程度可以用变量的方差膨胀因子(variance inflation factor,VIF)度量VIF是衡量一个变量的回归系数的方差由共线性引起的膨胀比例。如果VIF值超过20表示共线性很严重。实际上VIF超过10则可能会有共线性问题需要处理。# 两个RDA结果的变量方差膨胀因子(VIF)# ********************************************# 本章第一个RDA结果包括所有环境因子变量vif.cca(spe.rda)alt penmoderate pensteep penvery_steep deb pH dur pho nit20.397021 2.085921 2.987679 4.594983 6.684187 1.363575 3.049937 30.614913 18.953092amm oxy dbo40.114746 6.854703 17.980889vif.cca(spechem.physio) # 偏RDAalt pen deb pH dur pho nit amm oxy dbo16.188416 1.873974 6.711460 1.205235 3.268620 25.363359 16.080319 30.685907 6.904214 17.782727# 使用双终止准则(Blanchet等2008a)前向选择解释变量# 1.包含所有解释变量的RDA全模型spe.rda.all Call: rda(formula spe.hel ~ alt pen deb pH dur pho nit amm oxy dbo, data env)Inertia Proportion RankTotal 0.5025 1.0000Constrained 0.3689 0.7341 10Unconstrained 0.1336 0.2659 18Inertia is varianceEigenvalues for constrained axes:RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA100.22803 0.05442 0.03382 0.03008 0.00749 0.00566 0.00443 0.00281 0.00138 0.00079Eigenvalues for unconstrained axes:PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC80.04642 0.02072 0.01746 0.01325 0.00974 0.00588 0.00512 0.00400(Showed only 8 of all 18 unconstrained eigenvalues)# 2.全模型校正R2(R2a.all [1] 0.5864353# 3.使用packfors 包内forward.sel()函数选择变量# library(packfor) #如果尚未载入packfors包需要运行这一步forward.sel(spe.hel, env, adjR2threshR2a.all)Testing variable 1Testing variable 2Testing variable 3Testing variable 4Procedure stopped (adjR2thresh criteria) adjR2cum 0.594764 with 4 variables (superior to 0.586435)variables order R2 R2Cum AdjR2Cum F pval1 alt 1 0.32806549 0.3280655 0.3031790 13.182488 0.0012 oxy 9 0.16402853 0.4920940 0.4530243 8.396715 0.0013 dbo 10 0.09733024 0.5894243 0.5401552 5.926448 0.0014 pen 2 0.06323025 0.6526545 0.5947636 4.368924 0.007#注意正如这个函数的提示信息所示选择最后一个变量其实违背了#adjR2thresh终止准则所以pen变量最终不应该在被选变量列表中。# 使用vegan包内ordistep()函数前向选择解释变量。该函数可以对因子变量进# 行选择也可以运行解释变量的逐步选择和后向选择。step.forward directionforward, pstep 1000)Start: spe.hel ~ 1Df AIC F Pr(F) alt 1 -28.504 13.1825 0.005 ** oxy 1 -27.420 11.7086 0.005 ** deb 1 -26.872 10.9840 0.005 ** nit 1 -26.716 10.7795 0.005 ** dbo 1 -23.172 6.4340 0.005 ** dur 1 -22.499 5.6673 0.005 ** pho 1 -22.189 5.3200 0.005 ** amm 1 -22.047 5.1620 0.005 ** pen 1 -20.155 3.1305 0.005 ** pH 1 -17.489 0.4839 0.815---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Step: spe.hel ~ altDf AIC F Pr(F) oxy 1 -34.620 8.3967 0.005 ** dbo 1 -32.103 5.5373 0.005 ** amm 1 -30.777 4.1281 0.010 ** pho 1 -30.560 3.9032 0.010 ** nit 1 -29.451 2.7810 0.035 * pen 1 -29.049 2.3847 0.040 * deb 1 -27.972 1.3504 0.170 dur 1 -27.954 1.3332 0.290 pH 1 -27.426 0.8403 0.515---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Step: spe.hel ~ alt oxyDf AIC F Pr(F) dbo 1 -38.789 5.9264 0.005 ** pho 1 -37.052 4.1280 0.010 ** amm 1 -36.527 3.6055 0.015 * pen 1 -36.399 3.4797 0.015 * dur 1 -34.740 1.8964 0.095 . deb 1 -34.388 1.5714 0.125 nit 1 -33.474 0.7474 0.620 pH 1 -33.035 0.3605 0.950---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Step: spe.hel ~ alt oxy dboDf AIC F Pr(F) pen 1 -41.639 4.3689 0.015 * dur 1 -39.394 2.2555 0.025 * deb 1 -38.436 1.4019 0.190 pho 1 -37.789 0.8420 0.455 nit 1 -37.577 0.6611 0.655 amm 1 -37.583 0.6656 0.700 pH 1 -37.316 0.4399 0.910---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Step: spe.hel ~ alt oxy dbo penDf AIC F Pr(F) dur 1 -41.502 1.5255 0.130 deb 1 -41.058 1.1528 0.380 pho 1 -40.822 0.9570 0.470 amm 1 -40.587 0.7641 0.630 nit 1 -40.560 0.7418 0.635 pH 1 -40.375 0.5912 0.760 RsquareAdj(rda(spe.hel ~ alt, dataenv))$adj.r.squared[1] 0.303179 RsquareAdj(rda(spe.hel ~ altoxy, dataenv))$adj.r.squared[1] 0.4530243 RsquareAdj(rda(spe.hel ~ altoxydbo, dataenv))$adj.r.squared[1] 0.5401552 RsquareAdj(rda(spe.hel ~ altoxydbopen, dataenv))$adj.r.squared[1] 0.5947636 #简约的RDA分析 # ************** spe.rda.pars spe.rda.parsCall: rda(formula spe.hel ~ alt oxy dbo, data env)Inertia Proportion RankTotal 0.5025 1.0000Constrained 0.2962 0.5894 3Unconstrained 0.2063 0.4106 25Inertia is varianceEigenvalues for constrained axes:RDA1 RDA2 RDA30.21802 0.05088 0.02729Eigenvalues for unconstrained axes:PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC80.05835 0.04357 0.02568 0.01740 0.01451 0.01230 0.00787 0.00646(Showed only 8 of all 25 unconstrained eigenvalues) anova.cca(spe.rda.pars, step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(formula spe.hel ~ alt oxy dbo, data env)Df Variance F Pr(F)Model 3 0.29619 11.963 0.001 ***Residual 25 0.20632---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 anova.cca(spe.rda.pars, step1000, byaxis)Permutation test for rda under reduced modelForward tests for axesPermutation: freeNumber of permutations: 999Model: rda(formula spe.hel ~ alt oxy dbo, data env)Df Variance F Pr(F)RDA1 1 0.218022 26.4181 0.001 ***RDA2 1 0.050879 6.1651 0.001 ***RDA3 1 0.027291 3.3069 0.002 **Residual 25 0.206319---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 vif.cca(spe.rda.pars)alt oxy dbo1.223386 3.544201 3.402515 (R2a.pars [1] 0.5401552# 1型标尺plot(spe.rda.pars, scaling1, displayc(sp, lc, cn),mainRDA三序图spe.helaltoxydbo - 1型标尺 - 拟合的样方坐标)spe4.sc scores(spe.rda.pars, choices1:2, scaling1, displaysp)arrows(0, 0, spe4.sc[, 1], spe4.sc[, 2], length0, lty1, colred)# 2型标尺plot(spe.rda.pars, displayc(sp, lc, cn),mainRDA三序图spe.helaltoxydbo - 2型标尺 - 拟合的样方坐标)spe5.sc scores(spe.rda.pars, choices1:2, displaysp)arrows(0, 0, spe5.sc[,1], spe5.sc[,2], length0, lty1, colred)#如果第三典范轴也显著可以选择绘制轴1和轴3、轴2和轴3的三序图。变差分解# 变差分解说明图par(mfrowc(1,3))showvarparts(2) # 两组解释变量showvarparts(3) #三组解释变量showvarparts(4) #四组解释变量# 1.带所有环境变量的变差分解spe.part.all spe.part.allPartition of variance in RDACall: varpart(Y spe.hel, X envchem, envtopo)Explanatory tables:X1: envchemX2: envtopoNo. of explanatory tables: 2Total variation (SS): 14.07Variance: 0.50251No. of observations: 29Partition table:Df R.squared Adj.R.squared Testable[ab] X1 7 0.60579 0.47439 TRUE[bc] X2 3 0.41526 0.34509 TRUE[abc] X1X2 10 0.73414 0.58644 TRUEIndividual fractions[a] X1|X2 7 0.24135 TRUE[b] 0 0.23304 FALSE[c] X2|X1 3 0.11205 TRUE[d] Residuals 0.41356 FALSE---Use function ‘rda’ to test significance of fractions of interestplot(spe.part.all, digits2)#这些图内校正R2是正确的数字但是韦恩图圆圈大小相同未与R2的大小成比例。 # 2.分别对两组环境变量进行前向选择 spe.chem R2a.all.chem forward.sel(spe.hel, envchem, adjR2threshR2a.all.chem, nperm9999)Testing variable 1Testing variable 2Testing variable 3Testing variable 4Procedure stopped (adjR2thresh criteria) adjR2cum 0.487961 with 4 variables (superior to 0.474388)variables order R2 R2Cum AdjR2Cum F pval1 oxy 6 0.30247973 0.3024797 0.2766456 11.708553 0.00012 dbo 7 0.09015052 0.3926303 0.3459095 3.859122 0.00433 nit 4 0.11522115 0.5078514 0.4487936 5.852965 0.00054 amm 5 0.05325801 0.5611094 0.4879610 2.912325 0.0083 spe.topo R2a.all.topo forward.sel(spe.hel, envtopo, adjR2threshR2a.all.topo, nperm9999)Testing variable 1Testing variable 2Testing variable 3Procedure stopped (alpha criteria): pvalue for variable 3 is 0.228900 (superior to 0.050000)variables order R2 R2Cum AdjR2Cum F pval1 alt 1 0.32806549 0.3280655 0.3031790 13.182488 0.00012 pen 2 0.05645105 0.3845165 0.3371717 2.384674 0.0469 # 解释变量简约组合(基于变量选择的结果) names(env)[1] alt pen deb pH dur pho nit amm oxy dbo envchem.pars envtopo.pars # 变差分解 (spe.part Partition of variance in RDACall: varpart(Y spe.hel, X envchem.pars, envtopo.pars)Explanatory tables:X1: envchem.parsX2: envtopo.parsNo. of explanatory tables: 2Total variation (SS): 14.07Variance: 0.50251No. of observations: 29Partition table:Df R.squared Adj.R.squared Testable[ab] X1 3 0.50785 0.44879 TRUE[bc] X2 2 0.38452 0.33717 TRUE[abc] X1X2 5 0.66351 0.59036 TRUEIndividual fractions[a] X1|X2 3 0.25318 TRUE[b] 0 0.19561 FALSE[c] X2|X1 2 0.14156 TRUE[d] Residuals 0.40964 FALSE---Use function ‘rda’ to test significance of fractions of interest plot(spe.part, digits2) # 所有可测部分的置换检验 # [ab]部分的检验 anova.cca(rda(spe.hel, envchem.pars), step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(X spe.hel, Y envchem.pars)Df Variance F Pr(F)Model 3 0.25520 8.5992 0.001 ***Residual 25 0.24731---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # [bc]部分的检验 anova.cca(rda(spe.hel, envtopo.pars), step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(X spe.hel, Y envtopo.pars)Df Variance F Pr(F)Model 2 0.19322 8.1216 0.001 ***Residual 26 0.30929---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # [abc]部分的检验 env.pars anova.cca(rda(spe.hel, env.pars), step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(X spe.hel, Y env.pars)Df Variance F Pr(F)Model 5 0.33342 9.0704 0.001 ***Residual 23 0.16909---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # [a]部分的检验 anova.cca(rda(spe.hel, envchem.pars, envtopo.pars), step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(X spe.hel, Y envchem.pars, Z envtopo.pars)Df Variance F Pr(F)Model 3 0.14020 6.3565 0.001 ***Residual 23 0.16909---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # [c]部分的检验 anova.cca(rda(spe.hel, envtopo.pars, envchem.pars), step1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(X spe.hel, Y envtopo.pars, Z envchem.pars)Df Variance F Pr(F)Model 2 0.078219 5.3197 0.001 ***Residual 23 0.169091---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #各个部分置换检验有不显著的吗 # 3.无变量nit(硝酸盐浓度)的变差分解 envchem.pars2 (spe.part2 Partition of variance in RDACall: varpart(Y spe.hel, X envchem.pars2, envtopo.pars)Explanatory tables:X1: envchem.pars2X2: envtopo.parsNo. of explanatory tables: 2Total variation (SS): 14.07Variance: 0.50251No. of observations: 29Partition table:Df R.squared Adj.R.squared Testable[ab] X1 2 0.39263 0.34591 TRUE[bc] X2 2 0.38452 0.33717 TRUE[abc] X1X2 4 0.65265 0.59476 TRUEIndividual fractions[a] X1|X2 2 0.25759 TRUE[b] 0 0.08832 FALSE[c] X2|X1 2 0.24885 TRUE[d] Residuals 0.40524 FALSE---Use function ‘rda’ to test significance of fractions of interest plot(spe.part2, digits2)RDA 作为多元方差分析(MANOVA)的工具# 基于RDA的双因素MANOVA# **************************# 生成代表海拔的因子变量(3个水平每个水平含9个样方)alt.fac # 生成近似模拟pH值的因子变量pH.fac 2, 1, 2, 3, 2, 1, 1, 3, 3))# 两个因子是否平衡table(alt.fac, pH.fac)table(alt.fac, pH.fac)# 用Helmert对照法编码因子和它们的交互作用项alt.pH.helm contrastslist(alt.faccontr.helmert, pH.faccontr.helmert))alt.pH.helm(Intercept) alt.fac1 alt.fac2 pH.fac1 pH.fac2 alt.fac1:pH.fac1 alt.fac2:pH.fac1 alt.fac1:pH.fac2 alt.fac2:pH.fac21 1 -1 -1 -1 -1 1 1 1 12 1 -1 -1 1 -1 -1 -1 1 13 1 -1 -1 0 2 0 0 -2 -24 1 -1 -1 1 -1 -1 -1 1 15 1 -1 -1 0 2 0 0 -2 -26 1 -1 -1 -1 -1 1 1 1 17 1 -1 -1 0 2 0 0 -2 -28 1 -1 -1 1 -1 -1 -1 1 19 1 -1 -1 -1 -1 1 1 1 110 1 1 -1 1 -1 1 -1 -1 111 1 1 -1 -1 -1 -1 1 -1 112 1 1 -1 0 2 0 0 2 -213 1 1 -1 0 2 0 0 2 -214 1 1 -1 1 -1 1 -1 -1 115 1 1 -1 -1 -1 -1 1 -1 116 1 1 -1 -1 -1 -1 1 -1 117 1 1 -1 1 -1 1 -1 -1 118 1 1 -1 0 2 0 0 2 -219 1 0 2 1 -1 0 2 0 -220 1 0 2 -1 -1 0 -2 0 -221 1 0 2 1 -1 0 2 0 -222 1 0 2 0 2 0 0 0 423 1 0 2 1 -1 0 2 0 -224 1 0 2 -1 -1 0 -2 0 -225 1 0 2 -1 -1 0 -2 0 -226 1 0 2 0 2 0 0 0 427 1 0 2 0 2 0 0 0 4attr(,assign)[1] 0 1 1 2 2 3 3 3 3attr(,contrasts)attr(,contrasts)$alt.fac[1] contr.helmertattr(,contrasts)$pH.fac[1] contr.helmert#检查当前对照法产生的表格哪一列代表海拔因子、pH值和交互作用项# 检查Helmert对照表属性1每个变量的和为1apply(alt.pH.helm[, 2:9], 2, sum)# 检查Helmert对照表属性2变量之间不相关cor(alt.pH.helm[, 2:9])alt.fac1 alt.fac2 pH.fac1 pH.fac2 alt.fac1:pH.fac1 alt.fac2:pH.fac1 alt.fac1:pH.fac2 alt.fac2:pH.fac2alt.fac1 1 0 0 0 0 0 0 0alt.fac2 0 1 0 0 0 0 0 0pH.fac1 0 0 1 0 0 0 0 0pH.fac2 0 0 0 1 0 0 0 0alt.fac1:pH.fac1 0 0 0 0 1 0 0 0alt.fac2:pH.fac1 0 0 0 0 0 1 0 0alt.fac1:pH.fac2 0 0 0 0 0 0 1 0alt.fac2:pH.fac2 0 0 0 0 0 0 0 1# 使用函数betadisper()(vegan包)(Marti Anderson检验)验证组内协方差矩阵# 的齐性spe.hel.d1 # 海拔因子(spe.hel.alt.MHV Homogeneity of multivariate dispersionsCall: betadisper(d spe.hel.d1, group alt.fac)No. of Positive Eigenvalues: 26No. of Negative Eigenvalues: 0Average distance to median:1 2 30.5208 0.5175 0.3881Eigenvalues for PCoA axes:PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA86.5329 1.7407 1.2269 1.0591 0.6117 0.4683 0.3987 0.3207anova(spe.hel.alt.MHV)Analysis of Variance TableResponse: DistancesDf Sum Sq Mean Sq F value Pr(F)Groups 2 0.1032 0.051602 0.8763 0.4292Residuals 24 1.4133 0.058889permutest(spe.hel.alt.MHV) # 置换检验Permutation test for homogeneity of multivariate dispersionsPermutation: freeNumber of permutations: 999Response: DistancesDf Sum Sq Mean Sq F N.Perm Pr(F)Groups 2 0.1032 0.051602 0.8763 999 0.439Residuals 24 1.4133 0.058889 # pH值因子 (spe.hel.pH.MHV Homogeneity of multivariate dispersionsCall: betadisper(d spe.hel.d1, group pH.fac)No. of Positive Eigenvalues: 26No. of Negative Eigenvalues: 0Average distance to median:1 2 30.6658 0.6716 0.7019Eigenvalues for PCoA axes:PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA86.5329 1.7407 1.2269 1.0591 0.6117 0.4683 0.3987 0.3207 anova(spe.hel.pH.MHV)Analysis of Variance TableResponse: DistancesDf Sum Sq Mean Sq F value Pr(F)Groups 2 0.00676 0.0033802 0.1587 0.8542Residuals 24 0.51124 0.0213018 permutest(spe.hel.pH.MHV) #置换检验Permutation test for homogeneity of multivariate dispersionsPermutation: freeNumber of permutations: 999Response: DistancesDf Sum Sq Mean Sq F N.Perm Pr(F)Groups 2 0.00676 0.0033802 0.1587 999 0.855Residuals 24 0.51124 0.0213018 #组内协方差齐性可以继续分析。# 首先检验交互作用项。海拔因子和pH因子构成协变量矩阵interaction.rda anova(interaction.rda, step1000, perm.max1000)#交互作用是否显著显著的交互作用表示一个因子的影响依赖另一个因子#的水平这将妨害主因子变量的分析。# 检验海拔因子的效应此时pH值因子和交互作用项作为协变量矩阵factor.alt.rda anova(factor.alt.rda, step1000, perm.max1000, stratapH.fac)#海拔因子影响是否显著#检验pH值因子的效应此时海拔值因子和交互作用项作为协变量矩阵factor.pH.rda alt.pH.helm[, c(2:3, 6:9)])anova(factor.pH.rda, step1000, perm.max1000, strataalt.fac)#pH值影响是否显著# RDA和显著影响的海拔因子三序图alt.rda.out plot(alt.rda.out, scaling1, displayc(sp, wa, cn),mainMultivariate ANOVA, factor altitude - scaling 1 - wa scores)spe.manova.sc arrows(0, 0, spe.manova.sc[, 1], spe.manova.sc[, 2], length0, colred)基于距离的RDA分析 # 基于距离的RDA分析(db-RDA) # **************************** # 1.分步计算 spe.bray spe.pcoa spe.scores # 交互作用的检验。从协变量矩阵获得Helmert对照编码海拔因子和pH值因子 interact.dbrda anova(interact.dbrda, step1000, perm.max1000)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(X spe.scores[1:27, ], Y alt.pH.helm[, 6:9], Z alt.pH.helm[, 2:5])Df Variance F Pr(F)Model 4 0.021857 0.441 1Residual 18 0.223016 #交互作用是否显著如果没有可以继续检验主因子的效应(此处未显示) # 2.直接使用vegan包内capscale()函数运行。只能以模型界面运行。响应变量 #可以是原始数据矩阵。 interact.dbrda2 anova(interact.dbrda2, step1000, perm.max1000)Permutation test for capscale under reduced modelPermutation: freeNumber of permutations: 999Model: capscale(formula spe[1:27, ] ~ alt.fac * pH.fac Condition(alt.fac pH.fac), distance bray, add TRUE)Df SumOfSqs F Pr(F)Model 4 0.4667 0.4811 1Residual 18 4.3658 # 或者响应变量可以是相异矩阵。 interact.dbrda3 anova(interact.dbrda3, step1000, perm.max1000)Permutation test for capscale under reduced modelPermutation: freeNumber of permutations: 999Model: capscale(formula spe.bray ~ alt.fac * pH.fac Condition(alt.fac pH.fac), add TRUE)Df SumOfSqs F Pr(F)Model 4 0.4667 0.4811 1Residual 18 4.3658非线性关系的RDA分析 # 二阶解释变量的RDA # ******************* # 生成das和das正交二阶项(由poly()函数获得)矩阵 das.df colnames(das.df) # 验证两个变量是否显著 forward.sel(spe, das.df)Testing variable 1Testing variable 2variables order R2 R2Cum AdjR2Cum F pval1 das 1 0.44777219 0.4477722 0.4273193 21.892865 0.0012 das2 2 0.07870749 0.5264797 0.4900550 4.321662 0.005 # RDA和置换检验 spe.das.rda anova(spe.das.rda)Permutation test for rda under reduced modelPermutation: freeNumber of permutations: 999Model: rda(formula spe ~ das das2, data as.data.frame(das.df))Df Variance F Pr(F)Model 2 36.304 14.454 0.001 ***Residual 26 32.652---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # 三序图(拟合的样方坐标2型标尺) plot(spe.das.rda, scaling2, displayc(sp, lc, cn), mainRDA三序图spedas das2 - 2型标尺 - 拟合的样方坐标) spe6.sc arrows(0, 0, spe6.sc[, 1], spe6.sc[, 2], length0, lty1, colred)# 4种鱼类的分布地图# ******************par(mfrowc(2, 2))plot(spa$x, spa$y, asp1, colbrown, cexspe$TRU,xlabx (km), ylaby (km), main褐鳟)lines(spa$x, spa$y, collight blue)plot(spa$x, spa$y, asp1, colbrown, cexspe$OMB,xlabx (km), ylaby (km), main鳟鱼)lines(spa$x, spa$y, collight blue)plot(spa$x, spa$y, asp1, colbrown, cexspe$ABL,xlabx (km), ylaby (km), main欧鮊鱼)lines(spa$x, spa$y, collight blue)plot(spa$x, spa$y, asp1, colbrown, cexspe$TAN,xlabx (km), ylaby (km), main鲤鱼)lines(spa$x, spa$y, collight blue)自写代码角为了能够正确自写RDA分析代码有必要参考Legendre和Legendre(1998)第11.1节相关内容。下面是计算步骤(基于协方差矩阵的RDA)1.计算中心化的物种数据矩阵与标准化解释变量矩阵的多元线性回归获得拟合值矩阵2.计算拟合值矩阵的PCA3.计算两类样方坐标4.结果输出。下面代码解释部分使用的公式编码与Legendre和Legendre(1998)一致。下面的代码集中在RDA约束部分目的是使读者对RDA数学过程感兴趣而不是最优化程序。myRDA # 1.数据的准备# *************Y.mat Yc X.mat Xcr # 2.多元线性回归的计算# *********************# 回归系数矩阵 (eq. 11.4)B # 拟合值矩阵(eq. 11.5)Yhat # 残差矩阵Yres #维度n p m # 3. 拟合值PCA分析# ******************# 协方差矩阵 (eq. 11.7)S # 特征根分解eigenS # 多少个典范轴kc 0.00000001))# 典范轴特征根ev # 矩阵Yc(中心化)的总方差(惯量)trace sum(diag(cov(Yc)))# 正交特征向量(响应变量的贡献1型标尺)U row.names(U) # 样方坐标(vegan包内wa 坐标1型标尺eq. 11.12)F row.names(F) # 样方约束(vegan包内lc 坐标2型标尺eq. 11.13)Z row.names(Z) # 典范系数 (eq. 11.14)CC row.names(CC) # 解释变量# 物种-环境相关corXZ # 权重矩阵的诊断D # 解释变量双序图坐标coordX coordX2 row.names(coordX) row.names(coordX2) # 相对特征根平方根转化(为2型标尺)U2 row.names(U2) F2 row.names(F2) Z2 row.names(Z2) # 未校正R2R2 # 校正R2R2adj # 4.残差的PCA# *******************# 与第5章相同写自己的代码可以从这里开始...# eigenSres # evr # 5.输出Outputresult names(result) Species_sc1, wa_sc1, lc_sc1, Biplot_sc1, Species_sc2,wa_sc2, lc_sc2, Biplot_sc2)result}#将此函数应用到Doubs鱼类数据和环境数据的RDA分析doubs.myRDA summary(doubs.myRDA)Length Class ModeTotal_variance 1 -none- numericR2 1 -none- numericR2a 1 -none- numericCan_ev 10 -none- numericCan_coeff 100 -none- numericSpecies_sc1 270 -none- numericwa_sc1 290 -none- numericlc_sc1 290 -none- numericBiplot_sc1 100 -none- numericSpecies_sc2 270 -none- numericwa_sc2 290 -none- numericlc_sc2 290 -none- numericBiplot_sc2 100 -none- numeric
