( O$ ^/ p- p$ I- d2 m: s5 } 7 i% K0 G! u. R% U1 R1 U* i2 o# ~& S7 w9 Y" q# O
图2|环境因子数据,env.csv。8 p3 @, G' Z) J* @2 d
" D; r. H( D( a/ {8 \1 N
二、决策树回归模型 I: V3 U+ W: ?4 \; d! e* b: `当因变量为定量变量时,决策树进行回归分析。分析时先基于自变量划分预测空间Rj,此时残差平方和(residual sum of squares,RSS)值最小,然后对于落入某一预测空间的样本,做相同的预测。 $ V |4 }( v: q: m $ ]' w4 G3 F! F, h8 ^8 ]' [& u" f4 d; G) a) p" g
# ~5 V/ T# ]" ^( T/ K2 k
但是考虑到所有自变量的预测空间构建基本很难实现,因此常使用递归二进制拆分(recursive binary splitting)。每一步的拆分都使RSS减少量最大。 8 k$ i: }' g4 A- {" ?$ ~: Q4 l" m4 ]) R J: h; T1 x* }
如拆分点(cutpoint)为s,则拆分空间R1和R2为:5 f) p/ u& F' r* G
# H% f7 D+ ?6 y/ w0 j! |
8 \2 D# P2 u! U) C- X: J: ^9 X& H' g6 E; L. { X
j和s的选择基于使RSS最小化:: P! n7 U. T$ B8 Z# f7 s% |
* _% u2 x. t" D& ]. X( x+ p" h0 `& S" _9 s
?7 v5 [% I& k+ g2 o后面重复此过程,寻找最佳预测空间和拆分点,从而使每个结果区域的RSS值最小,直到达到终止拆分标准,比如每个终端节点包含的样本数都不高于设定的阈值。 % u5 x, S6 n( N ( x8 q9 }! F0 m1 r; I. D+ b2.1 构建回归决策树/ o- s" x2 r0 Y8 S: m! E
使用微生物数据与环境因子数据进行决策树回归分析。为了更好的评估分类树的分类性能,不能只计算训练误差,需要估计测试误差。将数据分为训练集和测试集数据,训练集数据用于构建模型,测试集数据用于模型评估。( Y# f& z9 A$ R, o( ^
, Y0 `5 a7 b4 J( c
# 2.1.1 将数据集分为train和test集,用train结果预测test的因变量值。" I, P! w7 _, x0 B
library(splitstackshape) y: D% g9 _6 n* Cspe = data.frame(ID = rownames(spe),spe)# stratified提取后,样本名会消失,先提取样本名,重新构建数据框。+ Z: d2 h; {$ |
1 M8 K& G& D3 S2 o/ [
## train data sets,每个分类提取相同数目的样本用作训练集# [* N- Y2 i$ I6 W
set.seed(12345)' U" h1 A& H" s7 k/ E# R
train.spe = stratified(spe, group=c("grazing"),size=10,replace=FALSE) 9 Z3 Y8 S4 g! e2 W1 ttable(train.spe$grazing) # 每个分类提取的样本数一致。' K H- v+ R9 Y
. A) L! f# Y* z& r9 u0 f
train.env = env[rownames(env) %in% train.spe$ID,]. P! d1 a; l( s5 u& c
table(train.env$grazing) # 每个分类提取的样本数一致。 8 _5 G: C" I$ }. l: M" `7 z6 I1 K8 {+ ^
## test data sets 6 K8 S0 R1 G, R! h6 }test.spe = spe[!spe$ID %in% train.spe$ID,]/ Y5 q4 R% `! I m
table(test.spe$grazing) * h- b6 u' t, d' B& L3 E: h7 v" w. O, G; k
test.env = env[!rownames(env) %in% rownames(train.env),]& v. N+ d8 K( W0 Y S0 y! R
table(test.env$grazing): l7 O. A8 k9 U; |6 _
* q4 s* e3 d1 T9 N3 r" _, \3 i#install.packages("tree")* P7 s, H/ L; }# R0 `* P# ]
library(tree) % @+ {7 _" y9 t% V" }- T( q! M0 @. l% i: u5 K
# 2.1.2 构建回归决策树 & ~5 E1 S6 B; Z) u9 \* Yreg.tre = tree::tree(train.env$env1 ~.,data=train.spe[,-c(1:3)]), b8 n5 X% \6 u: b, P k/ I) {3 j
reg.tre4 {+ F) j; f; e& S* ~' E
" S9 T1 B, U J: n. q' l8 Q1 @% |
# 2.1.3 输出结果简介 ) \- R& y& o z4 Q1 _## 输出表格的行为节点名(整数值表示),包含9列数据。# z5 u5 @& ]% o* w; F0 {( M
reg.tre$frame . }) m! D" d& v1 Z* `# V; x## 列包括var:用于拆分节点的变量及终端节点(<leaf>);) }5 O9 Y1 ?( W1 ]' S
reg.tre$frame$var0 d# t2 x- W0 y% f) D8 p* c. V8 @0 T
## n:每个节点的样本数量; v3 E8 L# A& v5 greg.tre$frame$n! _5 x$ `& T# ?) J, N' ^
## dev:每个节点的偏差( F4 H" [- Q/ ]% A' ^. j, @+ Z( a
reg.tre$frame$dev 5 f& d8 K6 P$ S0 l1 ]## yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平;/ q) k; d) A) \1 D8 v+ C
#mean(train.env[reg.tre$where == 4,3]) # 第四个节点包含样本的因变量均值。 / q# P* a( v! o- n. Zreg.tre$frame$yval) c( Y* u/ \' P3 g+ i: Y/ E% u
0 v N, A0 P' F5 o: x0 H
## split: 节点拆分,2列分别是属于左侧或右侧的标签;2 F5 ~" S C* n8 K1 h
reg.tre$frame$splits ( ~( T4 X y+ w5 Q## yprob:回归树,此为NULL;分类树则为因变量各水平的拟合比率,此数据有5个处理,所以有5列。 & ^2 z% y% `" \reg.tre$frame$yprob 3 S/ R1 U' [% d9 o3 o) M: g& N+ w# j" P: L2 O# V0 z+ X* ]9 T
## output,需要输出行名,则设置row.names=TRUE。 2 H) s" D% I: z4 k. R" fwrite.table(reg.tre$frame,"reg_tre_res.txt",sep="\t",quote = FALSE,row.names = FALSE)/ b- W6 K7 o' b* t* x
1 N4 z; ?! U( z( Z. P
## 每个样本所属节点 : v! y- r, w4 c7 b0 Vreg.tre$where7 b) o. Q" g2 x# K- O' j1 B" w& Q
## formul形式! z- [- k' k& g! w4 y
reg.tre$terms * A8 ?) }* I4 A" y## 自变量数据,x=FALSE则不会返回此数据 + Q3 V6 y' {1 E& V1 A2 h' i4 Nreg.tre$x. n7 [( m3 g* q0 b
## 因变量,y=FALSE则不会返回此数据% H0 Q) p# u. f$ E
reg.tre$y1 Q( r- T* h( T
## 样本权重,未设置则均为1,权重值可以为分数形式。 q. D& e! o4 L6 greg.tre$weights# v- Z5 e! ]! }
) u1 |; T0 N7 t, g% W
## 结果描述统计2 ]" W I: u. W" \
reg.tre.res = summary(reg.tre), h1 w: h6 o8 p+ {
reg.tre.res/ m) `# r3 p' o# x) X2 z# b
reg.tre.res$used # 用于构建回归决策树的自变量1 z! o9 Y" `1 h6 C# \9 }/ i" R# ^
reg.tre.res$dev # 偏差,决策树的残差平方和。 ) u% D5 P0 e, Treg.tre.res$df # 训练样本数减去终端节点数 - J" w) [! c8 ]reg.tre.res$residuals# 每个训练样本因变量的残差 # ~0 G2 X3 ]* B 4 S- j/ A) w2 N- Q! X& S## 简单绘图 0 V# f2 A% b0 G' M7 yplot(reg.tre) " H; |, C7 k) s3 m, E: Ctext(reg.tre,pretty = 0)1 g$ j5 h6 M V! p6 H4 i
+ c8 q8 `$ a. a, ?8 C6 K8 z1 O5 B$ N
# 2.1.4 预测测试集数据* _4 J$ V T R- a! o+ c
reg.pred = predict(reg.tre,newdata = test.spe[,-c(1:3)]) 7 o4 s( Y2 _; J+ w9 x9 ?reg.pred & g2 G& k$ _2 B( p, g6 {
## 预测结果与原始结果绘图# C6 F( B; e) @
plot(reg.pred,test.env$env1) 3 }% }8 E9 v: W+ ^* zabline(0,1) * O* B" @5 C. G B. T5 @5 H* b; m . S7 E- I1 m- a" m, E## 计算残差平方和(MSE)和标准化均方误差(NMSE)/ h* d$ l/ T j1 ~# F
MSE0 = mean((reg.pred-test.env$env1)^2) & @8 k# {/ K/ P0 aMSE0 8 C2 P" i- m- c! d, J' sNMSE0 = mean((test.env$env1-reg.pred)^2)/mean((test.env$env1-mean(test.env$env1))^2) / d3 ~+ j, u# g0 ANMSE02 q& p/ g- R; c, T {
8 z& x- g& ~$ D, X; O
2 Z* s- D" K7 L7 M$ a v9 F
. m- [1 j5 U( s k! b图3|回归树构建结果,reg.tre。每个节点以整数标注,tree()默认树最大生长数值为31。因子变量的分类水平不能超过32。 * t* _, z" P, R0 I. k% E, V% Z% O! S
" z) k5 m) W6 Q1 ]( B' w
& @5 l* {8 ]5 t
图4|回归树输出结果,reg_tre_res.txt。var:用于拆分节点的变量及终端节点(<leaf>);n:每个节点的样本数量;dev:每个节点的偏差;yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平;split: 节点拆分,2列分别是属于左侧或右侧的标签。6 p1 _$ t! ^# ?( m6 @) w. q4 V
; v0 d+ l) r" y! Z" E' o* }! V" ~" }( r& T# l, J
; K! L# H2 c( Q5 M
图5|回归树输出结果描述统计,reg.tre.res。包括终端节点数、拆分使用变量和残差平方和均值等信息。+ P' v7 V: |, f+ o z+ d/ ~0 @
6 A% x* T6 d1 F- P# l, R- P2 J; B3 y
4 ?# f5 u% M8 [# w5 S
/ |" t z7 w& C
图6|简单回归树绘图。回归决策树的每个节点上的数值是该节点处因变量的均值。2 |# G+ x, }. u r0 \3 q# m, |4 N
$ h$ E3 ~; R: p) @( {% [/ M* x+ [
5 V, z1 x( f! s9 W3 o, K5 h
% W) z1 t9 K$ L9 X9 E/ e. W图7|测试数据因变量实际值与模型预测值散点图并添加趋势线。7 D& P' w, d3 o% c7 R
l2 n& t# B" _+ N/ e/ d `) }; s9 d1 r- o6 D2 ]
0 ?% n j( ~" N" K$ w5 t/ s4 o图8|均方误差与标准化均方误差。评价模型预测好坏的一个准则为标准化均方误差(normalized mean squares error,NMSE)。! W' e3 L0 ^ O/ n" u5 T( t
5 E9 M/ V: _* Z% ^+ d5 o