6 g" S5 N# Y' A7 M, P. U7 u- i / V1 U$ y6 B/ _9 M7 e* w7 M# f" O- C7 ]( R
图2|环境因子数据,env.csv。# ~* f8 u9 M& O( G9 u! U, X' X- E
3 p3 d% m, t+ X9 B8 t" S' g, J4 h二、决策树回归模型 & \" A+ }5 K2 I: D" A) Y+ J" C当因变量为定量变量时,决策树进行回归分析。分析时先基于自变量划分预测空间Rj,此时残差平方和(residual sum of squares,RSS)值最小,然后对于落入某一预测空间的样本,做相同的预测。5 Q5 _6 u* x( R7 d
1 a! f/ n; h2 d" I' G
_- B# e5 l6 ?& |' M6 X
2 ]+ S" q. h. f5 \7 G$ g
但是考虑到所有自变量的预测空间构建基本很难实现,因此常使用递归二进制拆分(recursive binary splitting)。每一步的拆分都使RSS减少量最大。& p0 V7 ]) A+ I
9 N+ ^# c. x" s如拆分点(cutpoint)为s,则拆分空间R1和R2为:: u6 o" p' m- @ f& A- c
' k2 y- }' @) w m7 m1 D: _ ( D m' T$ `3 C. H6 G 7 G3 `2 I/ `% d8 e9 p/ x/ Pj和s的选择基于使RSS最小化:- y+ ~6 @% L8 T0 h: p+ j0 m& P
d4 e5 D8 X6 _9 N
: f" Y1 w0 ?0 e v6 l# U 2 p" V2 ~" g( W& ^% U后面重复此过程,寻找最佳预测空间和拆分点,从而使每个结果区域的RSS值最小,直到达到终止拆分标准,比如每个终端节点包含的样本数都不高于设定的阈值。 ' O& e& S8 |9 n- F1 Z2 l * [& Q6 ~" t7 ?: ~# D2.1 构建回归决策树 6 \* |( v4 a4 Q使用微生物数据与环境因子数据进行决策树回归分析。为了更好的评估分类树的分类性能,不能只计算训练误差,需要估计测试误差。将数据分为训练集和测试集数据,训练集数据用于构建模型,测试集数据用于模型评估。; F7 h! m' a) o
5 p9 S" t E( m6 T! j; A5 J0 R
# 2.1.1 将数据集分为train和test集,用train结果预测test的因变量值。, z0 U7 M) K3 j4 u0 n( O9 l% H
library(splitstackshape)7 b! A A% K0 C) C' L' \
spe = data.frame(ID = rownames(spe),spe)# stratified提取后,样本名会消失,先提取样本名,重新构建数据框。 : d; B: K. f1 L! B1 h! C# t 0 D/ l! r7 E, X' j" z n% r## train data sets,每个分类提取相同数目的样本用作训练集 ) Q4 ]' S M7 T# U5 vset.seed(12345) 4 T" z) e$ I5 j1 ttrain.spe = stratified(spe, group=c("grazing"),size=10,replace=FALSE) j# h @0 u2 l4 h: h! Ctable(train.spe$grazing) # 每个分类提取的样本数一致。 1 h$ L6 }* M2 ^ / f2 v2 K# d8 X" F- @6 b5 B% Qtrain.env = env[rownames(env) %in% train.spe$ID,]7 k# D3 B# R; K
table(train.env$grazing) # 每个分类提取的样本数一致。% C) t) e* G8 p: L8 R4 w5 Q
9 n/ x: g w F& l4 P## test data sets 5 `; e* Q# c& M. Q5 ]! ~test.spe = spe[!spe$ID %in% train.spe$ID,]- B# a" m4 r1 G# R
table(test.spe$grazing)) m% Y9 L2 ]' [$ N2 k& F3 k5 i4 u
1 ^( K1 P% v! V }
test.env = env[!rownames(env) %in% rownames(train.env),]" D& j/ N- S3 }5 {7 D/ F: ]+ o
table(test.env$grazing) 4 l2 j; l. ~0 W) B# D8 l& i - N- y1 u- n: w: U4 w& a) X& {' U#install.packages("tree") 8 Q4 v' B) A# jlibrary(tree) " u: ?3 b$ N5 T/ C( C& T1 u, f; }2 ]" b+ q; g' r8 s
# 2.1.2 构建回归决策树 ' ]) Z0 M( n2 R0 yreg.tre = tree::tree(train.env$env1 ~.,data=train.spe[,-c(1:3)])4 ] X! i- O+ e- J$ I. A+ M: C+ f2 [
reg.tre " U. p! D- u8 X2 G( r+ D7 y5 d# _7 F; Z) ^4 }( u
# 2.1.3 输出结果简介* s3 h5 Q* f, ^' C! d) y
## 输出表格的行为节点名(整数值表示),包含9列数据。, k3 F- z0 i' Q
reg.tre$frame 9 A6 `- z) A% @## 列包括var:用于拆分节点的变量及终端节点(<leaf>);+ o5 ]% b6 R- X8 J: V8 I
reg.tre$frame$var # M6 W3 e+ {- p, P* w$ `## n:每个节点的样本数量; + Z7 R' {' b8 e, Zreg.tre$frame$n # h+ d3 A+ P7 l' p8 ?## dev:每个节点的偏差 7 v6 A/ ~* m4 `( u. R ireg.tre$frame$dev * S. f" T0 ~) l R2 c! T/ c## yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平; ) u2 o% B( [1 d+ Q1 j! L/ K#mean(train.env[reg.tre$where == 4,3]) # 第四个节点包含样本的因变量均值。/ l# I5 v9 b% L7 h6 X" E& b
reg.tre$frame$yval ) g. d6 ^9 A* m. g* W ( `( B3 g3 E- J) u( T## split: 节点拆分,2列分别是属于左侧或右侧的标签; ' X$ m; z( w7 @* Z ^# Y! zreg.tre$frame$splits 5 n: i' B6 q- X9 S2 R! k1 b, \% ]## yprob:回归树,此为NULL;分类树则为因变量各水平的拟合比率,此数据有5个处理,所以有5列。8 u% S1 D0 g& p+ D# r4 H1 M" c- F
reg.tre$frame$yprob 1 O% l) l* O+ Q- x* t3 z- n7 ^5 M5 e: J
## output,需要输出行名,则设置row.names=TRUE。 8 C! F5 Z( M* x) s) Z" s1 E3 Xwrite.table(reg.tre$frame,"reg_tre_res.txt",sep="\t",quote = FALSE,row.names = FALSE), `% j9 J9 j1 J# d
' ]4 \ u( x9 }( |1 O## 每个样本所属节点 4 h0 }) X( \. v8 p, treg.tre$where ]2 U# T& J8 ~! p* ~7 S
## formul形式8 X, m$ f/ |; \& E3 }$ |' p# ?! R. S# t
reg.tre$terms . F$ F9 d0 P6 Z% ?+ i7 W## 自变量数据,x=FALSE则不会返回此数据 5 N! V' ^) N/ s [reg.tre$x* u1 \: X2 w& H' j
## 因变量,y=FALSE则不会返回此数据, i. ? G# C. o$ q. y
reg.tre$y K" V) p. Z$ Z8 I0 `1 W## 样本权重,未设置则均为1,权重值可以为分数形式。0 }5 n) R1 ]/ b& Z @4 ?
reg.tre$weights) ~5 _$ t2 V9 k# I K- Y
" b+ Q* Z w% o; a1 `2 M
## 结果描述统计" L; }/ y) r# \3 z) u& y- s: T
reg.tre.res = summary(reg.tre) & s/ [5 L8 _; h& R1 d. w; wreg.tre.res |4 a4 V. N- l. ?9 V" mreg.tre.res$used # 用于构建回归决策树的自变量9 L7 y9 K' c' P( n
reg.tre.res$dev # 偏差,决策树的残差平方和。2 [$ Z: s# N& {5 o, a
reg.tre.res$df # 训练样本数减去终端节点数4 s4 U; X. N/ q$ G- U( P8 y
reg.tre.res$residuals# 每个训练样本因变量的残差3 L1 c8 Q' T( f
7 C. w D5 @/ ~% V* w2 H) [1 z) E## 简单绘图, _1 {/ u1 I& S& T9 N, n6 l4 x
plot(reg.tre)4 s$ a8 L7 g7 V+ `& Q; D
text(reg.tre,pretty = 0)2 r3 _! N3 F! L3 C& U- t! A
- B9 V# r# P, ~& |# 2.1.4 预测测试集数据. n+ s {1 b8 y' Y
reg.pred = predict(reg.tre,newdata = test.spe[,-c(1:3)])" g+ P+ N2 a; ]) J" C$ c* j& f2 U
reg.pred 2 s! E# b2 F* U1 @5 n9 q4 ^; B; \8 j
## 预测结果与原始结果绘图3 |' w* ]9 @5 m6 X D$ z
plot(reg.pred,test.env$env1) 3 s$ t [8 y: ]" j& Aabline(0,1) . G& K" B: r+ w/ ~4 ^6 Z. Z8 ^5 x9 C! |& Q' F$ {( E
## 计算残差平方和(MSE)和标准化均方误差(NMSE). e+ w3 s2 a' R
MSE0 = mean((reg.pred-test.env$env1)^2) 0 d& s9 G+ Y) d p/ o6 H9 W7 T1 CMSE0, Q/ z1 N+ V" ^" b$ A
NMSE0 = mean((test.env$env1-reg.pred)^2)/mean((test.env$env1-mean(test.env$env1))^2)) r0 Q' s m5 ]" |- U
NMSE01 s6 O, m+ X* z/ a" A! E3 Q
" u/ h, M4 J/ S$ E M; ~ R
# k( S; g2 \1 \
3 g' H9 @+ f/ M6 ^3 q: S$ X
图3|回归树构建结果,reg.tre。每个节点以整数标注,tree()默认树最大生长数值为31。因子变量的分类水平不能超过32。 ; }' e* z. l3 d+ y* A& l; x9 t" i# u w
6 B ]- M7 m/ D O$ k* R + U, b/ i, O" U3 s5 ^0 }图4|回归树输出结果,reg_tre_res.txt。var:用于拆分节点的变量及终端节点(<leaf>);n:每个节点的样本数量;dev:每个节点的偏差;yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平;split: 节点拆分,2列分别是属于左侧或右侧的标签。9 U( |* ^9 e, X' `! W b2 H