5 ~2 v5 {. N7 i% R# L2 w5 L! ~! Z! T d! q T1 o$ A, O' e! H. e图1|原始otu表,spe.csv。前两列为分类信息。! J$ m2 ]2 b: e
' K4 B6 X( L5 f0 J- q, S K
- N9 T3 R2 N: F9 G4 r5 Y+ Z
. l0 a4 N- n( U. s8 W8 l图2|环境因子数据,env.csv。 & I: ?7 _ d" y0 [1 X( k5 k2 m 6 s4 j* g, V: T C1 b y6 X二、决策树回归模型 8 T% K0 z* H8 ]2 v2 T当因变量为定量变量时,决策树进行回归分析。分析时先基于自变量划分预测空间Rj,此时残差平方和(residual sum of squares,RSS)值最小,然后对于落入某一预测空间的样本,做相同的预测。! L" C! v* P7 v( f
4 }" X8 R8 y! |3 b
, u# ^6 o8 g& j) k* b + A- W5 ^/ r7 c) F- e+ }- `- H但是考虑到所有自变量的预测空间构建基本很难实现,因此常使用递归二进制拆分(recursive binary splitting)。每一步的拆分都使RSS减少量最大。 1 E* @" ?, R, l Y 7 @9 [7 J- Z7 o如拆分点(cutpoint)为s,则拆分空间R1和R2为: ! w/ g5 V" p3 w$ _( d ! S4 s$ u: f0 ?0 f* ^. K, {9 `, t' ]7 q: j7 l( x
3 W } y3 r& p0 R Hj和s的选择基于使RSS最小化:; w: o- O" a: [- s$ s
8 ^6 | d0 z$ [" I
6 K6 g; C' H7 i5 ~4 j$ M
; D) e7 n/ K, j& `' f( Y
后面重复此过程,寻找最佳预测空间和拆分点,从而使每个结果区域的RSS值最小,直到达到终止拆分标准,比如每个终端节点包含的样本数都不高于设定的阈值。3 W7 H1 r/ b- S
" ^/ ?- C, o! S2.1 构建回归决策树1 h5 U- @1 w9 O! r* F
使用微生物数据与环境因子数据进行决策树回归分析。为了更好的评估分类树的分类性能,不能只计算训练误差,需要估计测试误差。将数据分为训练集和测试集数据,训练集数据用于构建模型,测试集数据用于模型评估。 G U# O/ F) H6 a0 m: H$ [ , R) L0 K$ ~# }# W }# 2.1.1 将数据集分为train和test集,用train结果预测test的因变量值。 , }. `! A, x- [4 m: Y4 z. olibrary(splitstackshape) - s# f' b+ {6 T: p2 }+ t' ispe = data.frame(ID = rownames(spe),spe)# stratified提取后,样本名会消失,先提取样本名,重新构建数据框。 * {, Q+ t( S6 U$ d# z6 G, `( b q5 p9 N, t: J# E- p# ^
## train data sets,每个分类提取相同数目的样本用作训练集 8 A {4 r7 X7 F3 \/ U$ w, oset.seed(12345)7 M8 S5 s- m$ T9 ]* K2 h9 }
train.spe = stratified(spe, group=c("grazing"),size=10,replace=FALSE) 2 e6 F$ u+ r6 ^. {$ @table(train.spe$grazing) # 每个分类提取的样本数一致。0 T- F! |# H/ k/ W2 i( Y; G9 L
8 e5 t0 ?8 \9 c7 L0 n
train.env = env[rownames(env) %in% train.spe$ID,] 3 g% U. G! a- I1 Ltable(train.env$grazing) # 每个分类提取的样本数一致。; B5 ?* H+ Q0 Z7 P& X, A9 Q$ {* l
* r+ v; ]% y+ \# R1 u# D## test data sets' k1 h9 n- |5 o7 ]) C
test.spe = spe[!spe$ID %in% train.spe$ID,]# ]5 J& h4 ]( C* [/ ?/ w$ s6 L+ V4 V l
table(test.spe$grazing)2 \; _* \& c! z$ B& {
" Y! I; [( @- ^7 O q/ Z6 i P! gtest.env = env[!rownames(env) %in% rownames(train.env),] 3 o( I, O2 K% O2 C$ G2 L1 q% |table(test.env$grazing)! M& C8 E0 F5 I% n
. i) L% k+ L2 k' G
#install.packages("tree")6 j# w/ I7 m1 z+ d! D7 |
library(tree)! v8 M- ]$ p) Y6 k$ p0 v* a5 b! G
$ q' L- _2 \9 A* `3 {0 i# 2.1.2 构建回归决策树: M# R+ t. G2 o5 N( e
reg.tre = tree::tree(train.env$env1 ~.,data=train.spe[,-c(1:3)]) 3 ^4 x4 u f' v% t5 }8 R1 o% P) `2 f/ greg.tre. M1 k5 Z) G% a9 \: V( f5 t
" o% O1 p- G6 t3 t# 2.1.3 输出结果简介, _) h' W# G0 m8 \. L
## 输出表格的行为节点名(整数值表示),包含9列数据。2 C6 @7 j) v" R" h( w
reg.tre$frame 6 q* N: A; g- F H! u6 E p
## 列包括var:用于拆分节点的变量及终端节点(<leaf>);0 ^( a: _# o4 H7 w4 A5 W
reg.tre$frame$var- P" \; O% A6 e6 l
## n:每个节点的样本数量; & V1 c1 z0 c! e2 A, N1 J {- a9 wreg.tre$frame$n" Y5 }6 x2 @8 j8 s& P- i
## dev:每个节点的偏差 6 }) \7 m2 B* B8 Yreg.tre$frame$dev 5 L2 `+ R- m- O5 x/ O2 @+ ?4 z1 n% q! W## yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平; 5 m2 Y' n6 g, f' X5 J5 `% g#mean(train.env[reg.tre$where == 4,3]) # 第四个节点包含样本的因变量均值。6 e' c: W" A7 |( l7 U" @; e
reg.tre$frame$yval 2 u* O% i+ t3 W6 ?; \' c. }) |: z8 j
## split: 节点拆分,2列分别是属于左侧或右侧的标签;& v3 C/ K5 s' T5 l( b6 l
reg.tre$frame$splits5 g. B2 S4 R! H( `
## yprob:回归树,此为NULL;分类树则为因变量各水平的拟合比率,此数据有5个处理,所以有5列。% q5 o- `/ U/ ]( ?- N! g/ P' L! v
reg.tre$frame$yprob 4 D( D R3 M9 M1 q' T! R 9 v% \ i+ H1 \$ s( R## output,需要输出行名,则设置row.names=TRUE。/ q- u3 A! u2 U
write.table(reg.tre$frame,"reg_tre_res.txt",sep="\t",quote = FALSE,row.names = FALSE) 7 T0 ^; e. V1 M+ q1 r; s# Q9 C, |9 {! M4 f2 C, x
## 每个样本所属节点 / o* _% ^1 X& f% Kreg.tre$where' a7 d: x7 O" ]7 {& a" _! p
## formul形式7 L* K8 I8 l0 u/ i% z: K6 b# a% I5 `
reg.tre$terms # w3 T. q2 s7 b( e- m/ @4 d: A% k
## 自变量数据,x=FALSE则不会返回此数据8 x* O C$ [+ h) @
reg.tre$x) P' b. K& f! H
## 因变量,y=FALSE则不会返回此数据9 E6 L+ S4 G: e
reg.tre$y! Z* ^/ u' P# J$ @3 |% S. d
## 样本权重,未设置则均为1,权重值可以为分数形式。! J) ?3 G/ [3 D' K# g
reg.tre$weights3 i# [" r% j1 d$ ^, b) T9 O* J1 N
: a/ i' F( p7 O" J8 I! r# }& C6 ~## 结果描述统计; x5 ?; L6 o3 l! m
reg.tre.res = summary(reg.tre) 7 @1 N9 a! k4 }% W0 B, b) E5 O5 C4 jreg.tre.res 1 o- v% A# p: S# B* s' F5 Hreg.tre.res$used # 用于构建回归决策树的自变量/ R7 ~: g% l3 `/ {9 h# }
reg.tre.res$dev # 偏差,决策树的残差平方和。 % ?& n3 s4 F. Oreg.tre.res$df # 训练样本数减去终端节点数+ ^5 w6 i9 D8 `3 X8 F2 }) L
reg.tre.res$residuals# 每个训练样本因变量的残差 + D- P$ j7 R2 S" o1 f% _ + f0 T. U. m2 u* r## 简单绘图 3 W3 h& f2 s. pplot(reg.tre) 6 E5 \. m b# Z; Ktext(reg.tre,pretty = 0)3 E9 U2 z) U# j
; R) l+ T; Q+ U& \" ^# 2.1.4 预测测试集数据 1 r K' v3 i1 C! |1 Dreg.pred = predict(reg.tre,newdata = test.spe[,-c(1:3)]) 3 l: f6 g9 D; Q* W& F( N* k1 T2 |reg.pred - r! w6 _; G7 V" j6 d {$ [## 预测结果与原始结果绘图2 q2 W/ ?4 a/ g5 }$ N- C
plot(reg.pred,test.env$env1) ; s! P$ b9 X7 a$ V* g6 A: G4 {- Gabline(0,1) 4 U& ?. T2 k" p% i2 V+ J9 l+ d# y7 M- Q
## 计算残差平方和(MSE)和标准化均方误差(NMSE) 0 x) `# R" W/ YMSE0 = mean((reg.pred-test.env$env1)^2)% G+ p" r* N9 v7 b, s% V4 R
MSE0: A) C9 t8 ]/ O
NMSE0 = mean((test.env$env1-reg.pred)^2)/mean((test.env$env1-mean(test.env$env1))^2)% ~2 T: w' d3 s1 N5 x, b
NMSE0$ |6 `% L- M) R3 Z* a, T
" {. Y2 Z. k/ \+ R
2 I6 v# k# i e: k
7 z$ P4 W. F+ j; T, x& Y图3|回归树构建结果,reg.tre。每个节点以整数标注,tree()默认树最大生长数值为31。因子变量的分类水平不能超过32。- z2 W. K5 m$ }5 E
' C3 u- D1 `: l. ~, s x* {) s 2 m7 s# h3 i% E1 Z2 j2 [& Y: W' o8 Y
图4|回归树输出结果,reg_tre_res.txt。var:用于拆分节点的变量及终端节点(<leaf>);n:每个节点的样本数量;dev:每个节点的偏差;yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平;split: 节点拆分,2列分别是属于左侧或右侧的标签。 5 k; a) O7 [$ a$ M) c 7 [$ p: g' C* h9 F4 I+ ?2 }3 h7 H4 R& l7 v+ r( V
' @' H5 H* ^% |图5|回归树输出结果描述统计,reg.tre.res。包括终端节点数、拆分使用变量和残差平方和均值等信息。 ; N. l7 [1 P& o9 R1 K- V- ` ~. X6 ?
9 n' x t- H" r, [ G& ?7 f k
_# p3 d3 S& A' C o! M
图6|简单回归树绘图。回归决策树的每个节点上的数值是该节点处因变量的均值。& F. T2 y8 x, s5 y