数量生态学||多元回归树
2018-08-06 本文已影响4人
周运来就是我
多元回归树（Multivariate Regression Trees，MRT）是单元回归树的拓展，是一种对一系列连续型变量递归划分成多个类群的聚类方法，是在决策树(decision-trees)基础上发展起来的一种较新的分类技术。其目的是在定量或分类解释变量的控制下递归划分一个定量变量。这样的分类方式有时也称为约束聚类或者导向聚类。同一般回归模型一样，MRT也需要因变量(响应变量，群落学中一般是物种数据)和自变量(预测变量，一般为环境因子数据)。不同的是, MRT不需要在响应变量和预测变量之间建立参数估计的回归关系而是以预测变量为分类节点，利用二歧式的分割法(binary split), 将由响应变量定义的空间(样方)逐步划分为尽可能同质的类别,每一次划分都由某一预测变量(环境因子)的一个最佳划分值来完成, 将响应变量定义的空间(样方)分成两个部分(也叫节点, node)，最佳划分原则是让两个节点内部的差异尽量小，节点间的差异尽量大。
MRT是一种强大而可靠的分类方法，即使被划分的变量含有缺失值，或者响应变量与解释变量是非线性关系，或解释变量之间存在高阶相互关系，经过交叉验证等一系列筛选过程，多元回归树都能够发挥很好的预测作用。
MRT的计算有两个一起运行的程序组成：A数据约束划分；B分组结果交叉验证。
mvpart程序包中MRT的运算流程
- 将数据随机划分成k组（默认为10组，xval＝10）。
- 从k组中随机选取一组作为“验证组”（testing set），剩余k-1组（训练组，training set）重现混合后通过约束分析，按照组内平方和最小的原则，建立回归树。
- 将以上过程重复k-1次，即依次剔除一组数据。
- 共产生k个回归树，对于每个回归树的不同分类方案，将验证组（一组数据）内的对象分配到分组结果中。计算每个回归树分类方案的CVRE。
- 对回归树进行剪枝：可以保留CVRE最小的分类方案。也可以根据“1SE”准则，保留CVRE在最小的CVRE加1个标准差范围内最小的分类方案。
- 为了获得上面运行过程的误差估计值，需要重复多次（100次或500次）将对象随机分配为k组。
- 置换检验保留具有最小CVRE误差（或者是1SE值最小，此法最常用）的回归树。
# 加载所需的程序包
library(ade4)
library(vegan) # 应该先加载ade4后加载vegan以避免冲突
library(gclus)
library(cluster)
library(RColorBrewer)
library(labdsv)
library(mvpart)
library(MVPARTwrap) # MVPARTwrap这个程序包必须从本地zip文件安装
# 导入CSV格式的数据
rm(list = ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA")
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,]
# 物种多度数据：先计算样方之间的弦距离矩阵，然后进行单连接聚合聚类
# **************************************************************
spe.norm <- decostand(spe, "normalize")
spe.ch <- vegdist(spe.norm, "euc")
spe.ch.single <- hclust(spe.ch, method="single")
mvpart程序包从CRAN上下架了，感谢赖江山老师提供的Windows版安装包。我试了一下在3.4.0上安装成功了，本文也是基于这个版本来做的。
library(mvpart) # 载入mvpart程序包
spe.ch.mvpart <- mvpart(data.matrix(spe.norm) ~ ., env, margin=0.08, cp=0,
xv="pick", xval=nrow(spe), xvmult=100, which=4)
X-Val rep : 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
Minimum tree sizes
tabmins
2 3 4 6 7 8 9 10 11
5 9 8 1 21 17 10 23 6
参数xv="pick"允许通过人机交互的方式选择所需要的分类分组。如果自动选择具有最小CVRE值的回归树，可以设定xv="min"。
下面是尝试的选择分为4组回归树结果。
回归树特征描写：
- 每个节点通过一个解释变量划分。例如第一个节点被海拔划分为含有13和16个样方的两组。分割点（海拔340m）在原始数据中的比值，它是分割点左右两个样方的平均值。如果其他变量也能产生相同的分组，也可以选择其他变量进行分组。
- 每片叶子（终端）都显示所含样方的数量和RE，并显示该组内物种多度分布情况的条形图（物种排列顺序按照响应变量矩阵内的列变量顺序排序）
- 对于一个新的观测样方，可以根据其环境变量相应的值和回归树，将新样方直接分配到某一组。样方不断分组的过程中一个环境变量也可以多次使用。
> summary(spe.ch.mvpart)
Call:
mvpart(form = data.matrix(spe.norm) ~ ., data = env, xv = "pick",
xval = nrow(spe), xvmult = 100, margin = 0.08, which = 4,
cp = 0)
n= 29
CP nsplit rel error xerror xstd
1 0.44820889 0 1.0000000 1.0758371 0.05827462
2 0.09506051 1 0.5517911 0.7332388 0.11169179
3 0.08626143 2 0.4567306 0.7165753 0.10479606
4 0.04889908 3 0.3704692 0.7021106 0.10280365
Node number 1: 29 observations, complexity param=0.4482089
Means=0.05571,0.2515,0.27,0.2729,0.05579,0.06116,0.04772,0.07049,0.1295,0.1958,0.1094,0.06915,0.1578,0.1098,0.09459,0.07749,0.07079,0.05359,0.057,0.1183,0.05701,0.03674,0.09898,0.1923,0.07156,0.1933,0.06187, Summed MSE=0.5256098
left son=2 (16 obs) right son=3 (13 obs)
Primary splits:
alt < 340 to the right, improve=0.4482089, (0 missing)
das < 204.75 to the left, improve=0.4482089, (0 missing)
deb < 24.65 to the left, improve=0.4482089, (0 missing)
amm < 0.06 to the left, improve=0.3825997, (0 missing)
oxy < 9.55 to the right, improve=0.3619088, (0 missing)
Node number 2: 16 observations, complexity param=0.08626143
Means=0.09453,0.4461,0.441,0.4128,0.09132,0.09554,0.005803,0.04497,0.1255,0.17,0.04618,0.02865,0.07062,0.06948,0.06199,0.005803,0.01124,0.01141,0.01124,0.06578,0,0,0,0.08812,0,0.01161,0.005803, Summed MSE=0.3337866
left son=4 (13 obs) right son=5 (3 obs)
Primary splits:
amm < 0.11 to the left, improve=0.2462006, (0 missing)
dbo < 4.05 to the left, improve=0.2423770, (0 missing)
oxy < 10.25 to the right, improve=0.2423770, (0 missing)
pho < 0.15 to the left, improve=0.1819782, (0 missing)
nit < 0.51 to the left, improve=0.1785430, (0 missing)
Node number 3: 13 observations, complexity param=0.09506051
Means=0.007934,0.01205,0.05962,0.1008,0.01205,0.01885,0.09931,0.1019,0.1345,0.2275,0.1871,0.119,0.2651,0.1593,0.1347,0.1657,0.1441,0.1055,0.1133,0.1829,0.1272,0.08196,0.2208,0.3205,0.1596,0.4168,0.1309, Summed MSE=0.2361686
left son=6 (3 obs) right son=7 (10 obs)
Primary splits:
amm < 0.45 to the right, improve=0.4719501, (0 missing)
dbo < 10.6 to the right, improve=0.4719501, (0 missing)
pho < 1.07 to the right, improve=0.3942148, (0 missing)
oxy < 6.75 to the left, improve=0.3942148, (0 missing)
das < 236.15 to the right, improve=0.2239291, (0 missing)
Node number 4: 13 observations
Means=0.1092,0.5207,0.4824,0.4334,0.1053,0.1033,0,0.02678,0.09794,0.1276,0.03541,0.006695,0.06573,0.05029,0.03392,0,0.006695,0,0.006695,0.04208,0,0,0,0.01639,0,0,0, Summed MSE=0.2366276
Node number 5: 3 observations
Means=0.03095,0.1228,0.2613,0.3233,0.03095,0.0619,0.03095,0.1238,0.245,0.3539,0.09285,0.1238,0.09181,0.1527,0.1836,0.03095,0.03095,0.06086,0.03095,0.1685,0,0,0,0.3989,0,0.0619,0.03095, Summed MSE=0.3165241
Node number 6: 3 observations
Means=0,0,0,0,0,0,0.05206,0,0.07647,0.3167,0,0,0.205,0.07647,0,0,0.05206,0.07647,0,0,0,0,0.1806,0.3167,0.05206,0.7619,0, Summed MSE=0.1186808
Node number 7: 10 observations
Means=0.01031,0.01567,0.07751,0.131,0.01567,0.02451,0.1135,0.1325,0.1519,0.2007,0.2432,0.1547,0.2831,0.1842,0.1751,0.2154,0.1717,0.1142,0.1473,0.2378,0.1653,0.1065,0.2329,0.3217,0.1919,0.3133,0.1701, Summed MSE=0.1265172
> printcp(spe.ch.mvpart)
mvpart(form = data.matrix(spe.norm) ~ ., data = env, xv = "pick",
xval = nrow(spe), xvmult = 100, margin = 0.08, which = 4,
cp = 0)
Variables actually used in tree construction:
[1] alt amm
Root node error: 15.243/29 = 0.52561
n= 29
CP nsplit rel error xerror xstd
1 0.448209 0 1.00000 1.07584 0.058275
2 0.095061 1 0.55179 0.73324 0.111692
3 0.086261 2 0.45673 0.71658 0.104796
4 0.048899 3 0.37047 0.70211 0.102804
# MRT的残差
par(mfrow=c(1,2))
hist(residuals(spe.ch.mvpart), col="grey")
plot(predict(spe.ch.mvpart), residuals(spe.ch.mvpart)[,2],
main="Residuals vs Predicted")
abline(h=0, lty=3, col="grey")
> # 组的组成
> spe.ch.mvpart$where
[1] 3 3 3 3 4 3 3 4 3 3 3 3 3 3 3 4 7 7 7 7 7 6 6 6 7 7 7 7 7
> spe.ch.mvpart$where
[1] 3 3 3 3 4 3 3 4 3 3 3 3 3 3 3 4 7 7 7 7 7 6 6 6 7 7 7 7 7
> # 识别组的名称
> (groups.mrt <- levels(as.factor(spe.ch.mvpart$where)))
[1] "3" "4" "6" "7"
> # 第一片叶子的鱼类物种组成
> spe.norm[which(spe.ch.mvpart$where==groups.mrt[1]),]
CHA TRU VAI LOC OMB BLA HOT TOX
1 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0 0.0000000
2 0.0000000 0.7071068 0.5656854 0.4242641 0.0000000 0.0000000 0 0.0000000
3 0.0000000 0.5735393 0.5735393 0.5735393 0.0000000 0.0000000 0 0.0000000
4 0.0000000 0.4558423 0.5698029 0.5698029 0.0000000 0.0000000 0 0.0000000
6 0.0000000 0.3779645 0.5039526 0.6299408 0.0000000 0.0000000 0 0.0000000
7 0.0000000 0.6063391 0.4850713 0.6063391 0.0000000 0.0000000 0 0.0000000
10 0.0000000 0.1543033 0.6172134 0.6172134 0.0000000 0.0000000 0 0.0000000
11 0.1856953 0.5570860 0.7427814 0.1856953 0.1856953 0.0000000 0 0.0000000
12 0.2461830 0.6154575 0.4923660 0.4923660 0.2461830 0.0000000 0 0.0000000
13 0.2373563 0.5933908 0.5933908 0.2373563 0.3560345 0.2373563 0 0.0000000
14 0.2941742 0.4902903 0.4902903 0.3922323 0.3922323 0.2941742 0 0.0000000
15 0.2822163 0.3762883 0.3762883 0.4703604 0.1881442 0.3762883 0 0.0000000
16 0.1740777 0.2611165 0.2611165 0.4351941 0.0000000 0.4351941 0 0.3481553
VAN CHE BAR SPI GOU BRO
1 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
2 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
3 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.11470787
4 0.0000000 0.11396058 0.00000000 0.00000000 0.11396058 0.22792115
6 0.1259882 0.25197632 0.00000000 0.00000000 0.12598816 0.12598816
7 0.1212678 0.12126781 0.00000000 0.00000000 0.00000000 0.00000000
10 0.3086067 0.30860670 0.00000000 0.00000000 0.15430335 0.00000000
11 0.0000000 0.18569534 0.00000000 0.00000000 0.00000000 0.00000000
12 0.0000000 0.12309149 0.00000000 0.00000000 0.00000000 0.00000000
13 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
14 0.0000000 0.09805807 0.09805807 0.00000000 0.09805807 0.09805807
15 0.2822163 0.28221626 0.18814417 0.00000000 0.18814417 0.00000000
16 0.4351941 0.17407766 0.17407766 0.08703883 0.17407766 0.08703883
PER BOU PSO ROT CAR TAN BCO PCH GRE
1 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
2 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
3 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
4 0.22792115 0 0.00000000 0 0.00000000 0.11396058 0 0 0
6 0.12598816 0 0.00000000 0 0.00000000 0.25197632 0 0 0
7 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
10 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
11 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
12 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
13 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
14 0.00000000 0 0.00000000 0 0.00000000 0.00000000 0 0 0
15 0.00000000 0 0.00000000 0 0.00000000 0.09407209 0 0 0
16 0.08703883 0 0.08703883 0 0.08703883 0.08703883 0 0 0
GAR BBO ABL ANG
1 0.00000000 0 0 0
2 0.00000000 0 0 0
3 0.00000000 0 0 0
4 0.00000000 0 0 0
6 0.12598816 0 0 0
7 0.00000000 0 0 0
10 0.00000000 0 0 0
11 0.00000000 0 0 0
12 0.00000000 0 0 0
13 0.00000000 0 0 0
14 0.00000000 0 0 0
15 0.00000000 0 0 0
16 0.08703883 0 0 0
> # 第一片叶子的环境变量组成
> env[which(spe.ch.mvpart$where==groups.mrt[1]),]
das alt pen deb pH dur pho nit amm oxy dbo
1 0.3 934 48.0 0.84 7.9 45 0.01 0.20 0.00 12.2 2.7
2 2.2 932 3.0 1.00 8.0 40 0.02 0.20 0.10 10.3 1.9
3 10.2 914 3.7 1.80 8.3 52 0.05 0.22 0.05 10.5 3.5
4 18.5 854 3.2 2.53 8.0 72 0.10 0.21 0.00 11.0 1.3
6 32.4 846 3.2 2.86 7.9 60 0.20 0.15 0.00 10.2 5.3
7 36.8 841 6.6 4.00 8.1 88 0.07 0.15 0.00 11.1 2.2
10 99.0 617 9.9 10.00 7.7 82 0.06 0.75 0.01 10.0 4.3
11 123.4 483 4.1 19.90 8.1 96 0.30 1.60 0.00 11.5 2.7
12 132.4 477 1.6 20.00 7.9 86 0.04 0.50 0.00 12.2 3.0
13 143.6 450 2.1 21.10 8.1 98 0.06 0.52 0.00 12.4 2.4
14 152.2 434 1.2 21.20 8.3 98 0.27 1.23 0.00 12.3 3.8
15 164.5 415 0.5 23.00 8.6 86 0.40 1.00 0.00 11.7 2.1
16 185.9 375 2.0 16.10 8.0 88 0.20 2.00 0.05 10.3 2.7
> # 叶子的鱼类物种组成表格和饼图
> leaf.sum <- matrix(0, length(groups.mrt), ncol(spe))
> colnames(leaf.sum) <- colnames(spe)
> for(i in 1:length(groups.mrt)){
+ leaf.sum[i,] <- apply(spe.norm[which(spe.ch.mvpart$where==groups.mrt[i]),],
+ 2, sum)
+ }
> leaf.sum
CHA TRU VAI LOC OMB BLA HOT
[1,] 1.41970277 6.7687248 6.2714982 5.634304 1.36828926 1.3430130 0.00000000
[2,] 0.09284767 0.3682695 0.7839270 0.969990 0.09284767 0.1856953 0.09284767
[3,] 0.00000000 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.15617376
[4,] 0.10314212 0.1566709 0.7750802 1.310359 0.15667090 0.2450592 1.13488718
TOX VAN CHE BAR SPI GOU BRO
[1,] 0.3481553 1.2732731 1.6589502 0.4602799 0.08703883 0.8545320 0.6537141
[2,] 0.3713907 0.7349785 1.0616448 0.2785430 0.37139068 0.2754219 0.4579960
[3,] 0.0000000 0.2294157 0.9500115 0.0000000 0.00000000 0.6150052 0.2294157
[4,] 1.3247927 1.5186125 2.0069433 2.4324424 1.54705268 2.8307202 1.8419593
PER BOU PSO ROT CAR TAN BCO
[1,] 0.4409481 0.00000000 0.08703883 0.0000000 0.08703883 0.5470478 0.00000
[2,] 0.5508437 0.09284767 0.09284767 0.1825742 0.09284767 0.5053840 0.00000
[3,] 0.0000000 0.00000000 0.15617376 0.2294157 0.00000000 0.0000000 0.00000
[4,] 1.7514544 2.15448060 1.71685385 1.1420904 1.47306052 2.3783208 1.65332
PCH GRE GAR BBO ABL ANG
[1,] 0.000000 0.0000000 0.2130270 0.0000000 0.0000000 0.00000000
[2,] 0.000000 0.0000000 1.1968310 0.0000000 0.1856953 0.09284767
[3,] 0.000000 0.5417633 0.9500115 0.1561738 2.2856126 0.00000000
[4,] 1.065468 2.3285769 3.2165264 1.9190327 3.1330002 1.70137667
> par(mfrow=c(2,2))
> for(i in 1:length(groups.mrt)){
+ pie(which(leaf.sum[i,]>0), radius=1, main=c("leaf #", groups.mrt[i]))
+ }
> # 从mvpart（）函数获得的结果对象中提取MRT结果
> # 必须加载MVPARTwrap和rdaTest程序包
> spe.ch.mvpart.wrap <- MRT(spe.ch.mvpart, percent=10, species=colnames(spe))
> summary(spe.ch.mvpart.wrap)
Portion (%) of deviance explained by species for every particular node
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--- Node 1 ---
Complexity(R2) 44.82089
alt>=340 alt< 340
~ Discriminant species :
TRU VAI LOC ABL
% of expl. deviance 19.77533931 15.26700612 10.2177021 17.23790891
Mean on the left 0.44606214 0.44096408 0.4127684 0.01160596
Mean on the right 0.01205161 0.05962155 0.1007969 0.41681637
~ INDVAL species for this node: : left is 1, right is 2
cluster indicator_value probability
TRU 1 0.9128 0.001
VAI 1 0.8258 0.001
LOC 1 0.7535 0.001
CHA 1 0.4036 0.022
ABL 2 0.9729 0.001
GRE 2 0.9231 0.001
HOT 2 0.7994 0.001
PSO 2 0.7849 0.001
GAR 2 0.7844 0.001
BBO 2 0.7692 0.001
BOU 2 0.7432 0.001
ANG 2 0.7366 0.001
GOU 2 0.7289 0.002
CAR 2 0.6998 0.001
ROT 2 0.6942 0.002
BCO 2 0.6923 0.001
SPI 2 0.6200 0.002
BAR 2 0.6170 0.002
BRO 2 0.5892 0.021
TAN 2 0.5658 0.018
PCH 2 0.5385 0.002
PER 2 0.5268 0.023
TOX 2 0.4803 0.014
Sum of probabilities = 0.786
Sum of Indicator Values = 17.76
Sum of Significant Indicator Values = 16.16
Number of Significant Indicators = 23
Significant Indicator Distribution
1 2
4 19
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--- Node 2 ---
Complexity(R2) 8.626143
amm< 0.11 amm>=0.11
~ Discriminant species :
TRU GAR
% of expl. deviance 29.3525828 27.13056017
Mean on the left 0.5206711 0.01638669
Mean on the right 0.1227565 0.39894367
eliminating null columns
~ INDVAL species for this node: : left is 1, right is 2
cluster indicator_value probability
TRU 1 0.8092 0.005
GAR 2 0.9605 0.001
TAN 2 0.8001 0.014
Sum of probabilities = 6.17
Sum of Indicator Values = 10.26
Sum of Significant Indicator Values = 2.57
Number of Significant Indicators = 3
Significant Indicator Distribution
1 2
1 2
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--- Node 3 ---
Complexity(R2) 9.506051
amm>=0.45 amm< 0.45
~ Discriminant species :
ABL
% of expl. deviance 32.0463373
Mean on the left 0.7618709
Mean on the right 0.3133000
~ INDVAL species for this node: : left is 1, right is 2
cluster indicator_value probability
ABL 1 0.7086 0.007
BAR 2 1.0000 0.004
SPI 2 1.0000 0.002
PER 2 1.0000 0.004
BOU 2 1.0000 0.002
CAR 2 1.0000 0.005
TAN 2 1.0000 0.004
ANG 2 1.0000 0.004
LOC 2 0.9000 0.019
TOX 2 0.9000 0.016
BCO 2 0.9000 0.015
PSO 2 0.7673 0.030
BBO 2 0.7080 0.037
BRO 2 0.7066 0.032
HOT 2 0.6855 0.032
VAN 2 0.6651 0.035
Sum of probabilities = 6.279
Sum of Indicator Values = 18.74
Sum of Significant Indicator Values = 13.94
Number of Significant Indicators = 16
Significant Indicator Distribution
1 2
1 15
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--- Final partition ---
Sites in the leaf # 3
[1] "1" "2" "3" "4" "6" "7" "10" "11" "12" "13" "14" "15" "16"
Sites in the leaf # 4
[1] "5" "9" "17"
Sites in the leaf # 6
[1] "23" "24" "25"
Sites in the leaf # 7
[1] "18" "19" "20" "21" "22" "26" "27" "28" "29" "30"
~ INDVAL species of the final partition:
cluster indicator_value probability
TRU 1 0.7900 0.001
VAI 1 0.5422 0.001
LOC 1 0.4506 0.035
ABL 3 0.6700 0.001
BCO 4 0.9000 0.001
BOU 4 0.8744 0.001
ANG 4 0.8461 0.001
CAR 4 0.7965 0.001
BBO 4 0.7080 0.001
PCH 4 0.7000 0.009
PSO 4 0.6568 0.001
BAR 4 0.6548 0.001
HOT 4 0.5776 0.020
GRE 4 0.5632 0.027
SPI 4 0.5424 0.030
TAN 4 0.5304 0.006
GOU 4 0.4385 0.042
Sum of probabilities = 2.848
Sum of Indicator Values = 14.7
Sum of Significant Indicator Values = 11.24
Number of Significant Indicators = 17
Significant Indicator Distribution
1 3 4
3 1 13
~ Corresponding MRT leaf number for Indval cluster number:
MRT leaf # 3 is Indval cluster no. 1
MRT leaf # 4 is Indval cluster no. 2
MRT leaf # 6 is Indval cluster no. 3
MRT leaf # 7 is Indval cluster no. 4
# MRT作为空间和时间系列约束聚类方法
spe.ch.seq <- mvpart(as.matrix(spe) ~ das, env, cp=0, xv="pick", margin=0.08,
xval=nrow(spe), xvmult=100, which=4)
summary(spe.ch.seq)
# 组的组成（终端节点的标识）
(gr <- spe.ch.seq$where)
# 重新编排聚类簇的编号
aa <- 1
gr2 <- rep(1,length(gr))
for (i in 2:length(gr)) {
if (gr[i]!=gr[i-1]) aa <- aa+1
gr2[i] <- aa
}
# 在Doubs河地图上标注样方所属的聚类簇
plot(spa, asp=1, type="n", main="MRT 分类组",
xlab="x坐标 (km)", ylab="y坐标 (km)")
lines(spa, col="light blue")
text(70, 10, "上游", cex=1.2, col="red")
text(15, 115, "下游", cex=1.2, col="red")
k <- length(levels(factor(gr2)))
for (i in 1:k) {
points(spa[gr2==i,1], spa[gr2==i,2], pch=i+20, cex=3, col=i+1, bg=i+1)
}
text(spa, row.names(spa), cex=0.8, col="white", font=2)
legend("bottomright", paste("组",1:k), pch=(1:k)+20, col=2:(k+1),
pt.bg=2:(k+1), pt.cex=2, bty="n")
