) ? l4 C3 c( j/ J2 Y# k* a # D$ a x! c- k图1|原始otu表,spe.csv。前两列为分类信息。- A2 R$ J+ V! t$ s
0 P) j! c3 ~" g
0 g. B" k$ ?: M$ j9 L
, r Z2 @" y4 _+ C q$ C
图2|环境因子数据,env.csv。 6 E+ ~$ }! E/ a4 e0 g1 |* M( u1 I: ]4 p" [ K' R @; q* n, \) g
二、决策树回归模型 3 S/ l; w. X6 x* a% d) Y+ H当因变量为定量变量时,决策树进行回归分析。分析时先基于自变量划分预测空间Rj,此时残差平方和(residual sum of squares,RSS)值最小,然后对于落入某一预测空间的样本,做相同的预测。% U* B$ D, U/ {. W2 C; h
4 Q! O* c$ z: t8 {) `7 j) @& y5 J( b
" q# }, f+ D: o* V( r3 u
但是考虑到所有自变量的预测空间构建基本很难实现,因此常使用递归二进制拆分(recursive binary splitting)。每一步的拆分都使RSS减少量最大。 , }0 I; G6 h p; ~3 P2 H5 v1 T! w" {- P& V
如拆分点(cutpoint)为s,则拆分空间R1和R2为: 8 H7 J# a+ v7 |, Y7 s6 J 1 W/ r9 c- N$ R ) `' q( H. [ Z/ P) x. ?2 v F% ?( Q8 N7 S7 x
j和s的选择基于使RSS最小化:8 B8 ^0 R i- J3 W! Z
# e$ O0 W" H6 Q5 {/ B I ( ^5 _: N! ]* u1 a# e: j0 L8 T D ! i( N" E1 |3 |1 ]% n* F后面重复此过程,寻找最佳预测空间和拆分点,从而使每个结果区域的RSS值最小,直到达到终止拆分标准,比如每个终端节点包含的样本数都不高于设定的阈值。 b: P6 C6 D0 w) Z: \- g. V
% C1 U3 I4 u9 S! ^
2.1 构建回归决策树 - Y/ _% S7 @& F! P ~) y' L使用微生物数据与环境因子数据进行决策树回归分析。为了更好的评估分类树的分类性能,不能只计算训练误差,需要估计测试误差。将数据分为训练集和测试集数据,训练集数据用于构建模型,测试集数据用于模型评估。 5 ?' Z2 |, v' o1 t! q: {" {" `- z* P$ _& x1 U
# 2.1.1 将数据集分为train和test集,用train结果预测test的因变量值。 ( `4 k. M8 U' g" W* U e$ Mlibrary(splitstackshape)0 a* g i- @" ~6 c0 v4 @" m( x
spe = data.frame(ID = rownames(spe),spe)# stratified提取后,样本名会消失,先提取样本名,重新构建数据框。" X. O5 M g3 h" u# z0 T8 H8 _
! k7 U/ F; T# l
## train data sets,每个分类提取相同数目的样本用作训练集 * A }4 N' e3 S# F$ [4 q* eset.seed(12345) ]! `' ~7 l! l# v0 ]9 o; D; i& T/ z5 Rtrain.spe = stratified(spe, group=c("grazing"),size=10,replace=FALSE)) s- z6 Y0 H* I
table(train.spe$grazing) # 每个分类提取的样本数一致。 - z2 i. b* i' E' D) N8 a* I: Q9 p; X# m% S1 s, L2 p' |
train.env = env[rownames(env) %in% train.spe$ID,]& P( @9 L5 O0 x7 Z3 b
table(train.env$grazing) # 每个分类提取的样本数一致。) h' l& S# |" }( C& r, v
) F1 w. r) O& y) L5 J
## test data sets' y' ^+ ]# i# H9 ^( T
test.spe = spe[!spe$ID %in% train.spe$ID,] ! t+ \6 a& m. ltable(test.spe$grazing) ( R e3 J, b) v! M# l" D+ y* ~- b/ C! N5 ]7 C, N
test.env = env[!rownames(env) %in% rownames(train.env),] , T2 `9 Y( ?! W4 s6 b# B0 F: ctable(test.env$grazing)6 }% w. R) c+ l* b
5 s4 d' C" R6 B" U1 J. a1 U+ `#install.packages("tree") ! v2 U Y" Q8 l5 m( hlibrary(tree) 2 M3 U% f ~- y2 L& z& l5 s1 W& b4 e$ S! o# O. K/ T; M
# 2.1.2 构建回归决策树8 e/ S. v) z' @' @$ A
reg.tre = tree::tree(train.env$env1 ~.,data=train.spe[,-c(1:3)])+ x, F& g' z8 v
reg.tre ' E+ e( _' Z8 u% R1 s1 R% t3 l! e, Y D2 [
# 2.1.3 输出结果简介 " X3 a4 |5 u) b6 T+ M## 输出表格的行为节点名(整数值表示),包含9列数据。9 d- I* X- S( X" v( h2 \
reg.tre$frame . Q; B; [0 E. M1 p/ R7 K& E## 列包括var:用于拆分节点的变量及终端节点(<leaf>);# k7 _! s' K8 G8 H6 t- }) Z' u
reg.tre$frame$var - d3 p6 c* {2 e# `/ R## n:每个节点的样本数量; 9 B) X: z0 I' F0 Z3 ]2 freg.tre$frame$n 0 ~$ j8 n- a3 x- R' {0 c% @## dev:每个节点的偏差, q1 ~ g% ?1 m& o) K5 e3 L
reg.tre$frame$dev / P- u7 @. @: u2 I6 {0 _0 w1 t## yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平;+ [! r9 c! r6 W$ n# ^9 [, H
#mean(train.env[reg.tre$where == 4,3]) # 第四个节点包含样本的因变量均值。- Y! `7 F' E7 k5 _- t
reg.tre$frame$yval7 z: e; y# }# ?
2 O- W1 ]1 T1 P: F' b) |
## split: 节点拆分,2列分别是属于左侧或右侧的标签; , ?4 |% [5 {3 Y# u" D" ireg.tre$frame$splits6 G; l# w7 e# Q
## yprob:回归树,此为NULL;分类树则为因变量各水平的拟合比率,此数据有5个处理,所以有5列。% u, }- K. D& @
reg.tre$frame$yprob4 u1 m% s* N0 a/ c- K
9 z( E4 \. _) ]0 ?2 o* H6 B## output,需要输出行名,则设置row.names=TRUE。( U4 X u# ~2 t0 x. w( Q
write.table(reg.tre$frame,"reg_tre_res.txt",sep="\t",quote = FALSE,row.names = FALSE) 1 d* e$ }7 X: K2 ]0 n* `' {8 e6 S: F, U* }$ l. z7 H
## 每个样本所属节点 ( K" L! J* _; P/ Ereg.tre$where 1 P! p5 B8 B; Y5 c/ y## formul形式' q* R8 ~' O+ }3 N. u
reg.tre$terms & {% K0 R( c* [ V* R* ^- ]## 自变量数据,x=FALSE则不会返回此数据 1 E1 n; R2 z" _" Kreg.tre$x2 q8 j9 k8 B- f) B
## 因变量,y=FALSE则不会返回此数据. g6 G D1 d) s7 j, h
reg.tre$y0 V$ a( w2 _0 ~5 s, ?
## 样本权重,未设置则均为1,权重值可以为分数形式。 9 o1 B1 v# q( j c. d( ?* rreg.tre$weights , C$ k! E1 Z: \' p5 R( j$ U4 E/ ~6 I1 A b* p* R% o7 ], O
## 结果描述统计* e: b1 t; i, L/ Y3 q6 u6 s( w
reg.tre.res = summary(reg.tre) - z$ l3 Q2 F* [) b) M1 d* N6 Rreg.tre.res 0 y( b+ y9 a- Xreg.tre.res$used # 用于构建回归决策树的自变量+ D; s P! p: Y
reg.tre.res$dev # 偏差,决策树的残差平方和。; ]/ @1 Y( L& D& e) G
reg.tre.res$df # 训练样本数减去终端节点数3 A1 T, C p+ P4 T" H4 _9 y; y
reg.tre.res$residuals# 每个训练样本因变量的残差 , b! K8 O: W4 A' f0 L# Q: T / s: Y3 {3 R- n# ~ `## 简单绘图 . b4 j# |" q% M( [* Splot(reg.tre)2 `! k' V0 v1 @" g, ]) U
text(reg.tre,pretty = 0) ; v3 q: M$ a( I. L) R, P; R; O7 i3 O# F0 n" h) j( q( G0 D1 N
# 2.1.4 预测测试集数据 ! {6 a: Z) o7 Mreg.pred = predict(reg.tre,newdata = test.spe[,-c(1:3)])# k( Z. j1 n/ X/ f$ l
reg.pred : x& }9 u& Q% ]' S( Y5 y
## 预测结果与原始结果绘图+ e9 G7 T7 q# X, L) R
plot(reg.pred,test.env$env1) ! ^4 L6 e, ]" c# ~abline(0,1)" \7 n' k( A! \' Y' k. b8 E
0 {! }" ?$ M% l4 v9 S; p## 计算残差平方和(MSE)和标准化均方误差(NMSE)8 I0 i) j+ Y: T3 Y- N
MSE0 = mean((reg.pred-test.env$env1)^2) : Q( g/ [' d* @4 Y( H2 L: z; XMSE09 Z1 M% r( z' I' _! s* O: d; y
NMSE0 = mean((test.env$env1-reg.pred)^2)/mean((test.env$env1-mean(test.env$env1))^2) 0 o- m4 ~9 n2 ^& B& QNMSE06 K x7 B) i, Z+ v% I* n( d
% V U6 u+ B Z! _. v, A3 K3 y ; _" M9 U- [% ?3 W 8 J5 s2 G' H( g4 N/ o% S图3|回归树构建结果,reg.tre。每个节点以整数标注,tree()默认树最大生长数值为31。因子变量的分类水平不能超过32。" A! ?" o* ]+ D- c9 o1 B
0 m; \7 e! b4 x# F; g& Z# X
) N( i( A% M; A# p P2 T" g+ k. }7 d+ m$ w3 v1 A L3 y
图4|回归树输出结果,reg_tre_res.txt。var:用于拆分节点的变量及终端节点(<leaf>);n:每个节点的样本数量;dev:每个节点的偏差;yval:拟合结果,回归树为节点包含样本的因变量均值,分类树为该节点样本最多属于的分类水平;split: 节点拆分,2列分别是属于左侧或右侧的标签。6 h9 L) N& B( ?! C7 c, T, [
+ C! W0 [8 ] K3 v ' r; |' I( ~) B2 a' o: w6 F% a) ]* y. ?+ ~, ]; F7 g6 H
图5|回归树输出结果描述统计,reg.tre.res。包括终端节点数、拆分使用变量和残差平方和均值等信息。 " p- z/ m1 X! q8 a' g8 b * l% M w: R+ `" G/ t5 \+ x6 n+ i+ A3 Y7 x( X" J. ^
. H" L' n `2 j
图6|简单回归树绘图。回归决策树的每个节点上的数值是该节点处因变量的均值。 9 P5 R1 Q4 Q. s. ?$ r2 C0 G" U [; ~7 x6 V3 P
% X8 B) V+ F) |/ q" g& P# e4 k- ^: _8 }
图7|测试数据因变量实际值与模型预测值散点图并添加趋势线。/ X5 g* p4 }2 V I. `) {- s$ g2 N
% e6 {$ N! T$ a+ d$ y( c2 [9 Z8 N* M
8 c' K1 Q- V2 r* h; V* w图8|均方误差与标准化均方误差。评价模型预测好坏的一个准则为标准化均方误差(normalized mean squares error,NMSE)。, j8 ?$ z) K, o0 Z
8 o/ z% U1 E. g! M. T5 ?. ?8 s ]7 p
$ N2 q2 U5 ^# I A) W. P2 y
) m: v3 j4 n |/ E8 g8 a分母表示用最简单的算术平均来预测y的残差平方和。分子为该模型拟合后的残差平方和。此模型的NMSE不小于1,说明此回归模型没有任何意义(NMSE≥1)。此处是虚构数据,只讲使用方法,产生的模型没有任何意义,也没有影响。/ D* H5 `2 _% N
1 `; ^7 R0 ]8 c" u# ?* _. I4 Z2.2 优化模型-剪枝(Tree Pruning)0 ^3 Q( |/ [0 r* L0 H' h% ~
经过上述过程,构建的模型可能会过拟合,导致模型对训练集数据有很好的预测能力,但对测试集数据的预测能力较差。可能的原因是生成的决策树模型过于复杂。- {% g( q A- B3 T3 S9 B+ I
/ a& v% W5 W3 x5 F0 \
解决的方法之一是仅在拆分能使RSS降低值超过某个阈值的情况下才继续进行拆分。但是低于某个阈值的拆分点的之后的拆分点可能降低RSS的能力很强,所以不能随意剪枝。所以更好的优化模型的方式是先不设定RSS阈值,构建一个较大的决策树T0,然后根据某种方法对其进行剪枝获得子树。 . [( a# W; x1 ~' ^( U v m& E6 x5 b4 l! D
剪枝的标准主要为获得的子树的错误率最低,常用交叉验证选择具有最低错误率的子树。但是子树集一般很大,所以一般限定在一个更小的子树集中进行交叉验证。这里引入一个新的概念复杂性代价剪枝(Cost complexity pruning,或最弱链接剪枝(weakest link pruning))。此时剪枝不考虑每棵子树,而只考虑由非负调整参数α索引的树序列。然后基于交叉验证选择α。α控制着树的复杂性及树与训练数据的适配性之间的权衡,当α=0时,子树T就是T0;当α的值逐渐增大时,表示拥有许多终端节点的树要付出的复杂性代价,因此终端节点越少,α值将会越小。4 |3 e- x0 I8 @) n4 M
% Q& x* [) l1 y1 e0 f
这里介绍两个定义方差(variance)与偏差(bias):1)方差是训练数据集的预测值或预测分类水平相对于其他数据集的预测值或预测分类水平的离散程度,代表了模型的泛化能力。2)偏差是模型的预测值或预测分类水平与训练数据中的实际值或实际分类水平之间的差别,代表了模型的预测准确性。模型构建要在方差与偏差之间权衡,使总体误差(偏差+方差)最小。5 u! F( q( m) c# C/ T. }