多元回归树分析Multivariate Regression Trees，MRT

多元回归树（Multivariate Regression Trees，MRT）是单元回归树的拓展，是一种对一系列连续型变量递归划分成多个类群的聚类方法，是在决策树(decision-trees)基础上发展起来的一种较新的分类技术。同一般回归模型一样,MRT也需要因变量(响应变量，群落学中一般是物种数据)和自变量(预测变量，一般为环境因子数据)。不同的是, MRT不需要在响应变量和预测变量之间建立参数估计的回归关系,而是以预测变量为分类节点,利用二歧式的分割法(binary split), 将由响应变量定义的空间(样方)逐步划分为尽可能同质的类别,每一次划分都由某一预测变量(环境因子)的一个最佳划分值来完成, 将响应变量定义的空间(样方)分成两个部分(也叫节点, node),最佳划分原则是让两个节点内部的差异尽量小,节点间的差异尽量大。

MRT是一种强大而可靠的分类方法，即使被划分的变量含有缺失值，或者响应变量与解释变量是非线性关系，或解释变量之间存在高阶相互关系，经过交叉验证等一系列筛选过程，多元回归树都能够发挥很好的预测作用。

多元回归树结果通常用一个倒立的树状结构图表示，上面的节点为分支点，分枝点是能够使得两分枝的响应变量的变异最大的预测变量的某个值，使得各节点（分类群）内样本的总方差最小或由于样本数量过少无法继续分层，这里的终节点为叶，而分枝开始的节点被称为根。在大量的类群划分方案中，通常保留“最具预测能力的”划分方案。

回归树建模的原理与方法

多元回归树的计算通常分两步来完成，最初生成一颗较大的树，然后通过统计估计来删除底部的一些节点对树进行修剪，以防止过度拟合并保留最佳的分类方案。

1.  数据约束划分

在预测变量（连续型或分类变量）的控制下，通过持续的(或递归的)分层将响应变量（连续型变量）不断分类(亦即分枝)，其划分的依据为这一系列预测变量的解释变量。使得各节点（分类群）内样本的总方差最小或由于样本数量过少无法继续分层。此时，保留的是相对误差（relative error， RE)最小的回归树，但是，这种情况下，回归树只具有解释功能，而缺乏预测功能。相对误差为所有叶子组哪平方和除以原始数据的平方和，也就是回归树不能解释的方差比例。

2.  交叉验证和回归树的剪枝

为获得最最具预测能力的分类方案，第一步产生的回归树还需要通过剪枝处理。回归树的预测能力可以用其预测误差进行评估。评价的标准是既要保证回归树包含了足够的信息,又要把并不重要的枝节去掉。比较著名的规则就是“1SE”（1标准差）准则，其意思是: 首先要保证交叉验证误差（CVRE，通过交叉验证获得）尽量小，但不一定要取最小值, 而是允许它在“最小误差”加一个相应标准差的范围内，然后在此范围内选取尽量小的复杂性参量，进而以它为依据进行剪枝。这个规则体现了兼顾树的规模（复杂性）和误差大小的思想。

交叉验证通常是利用原始数据的一部分（称为训练集 training set）构建一颗树，剩下另一部分（称为验证集 test set）验证训练集构建的树的准确性，具有良好预测能力的回归树会将验证集合中的各数据点划分到合适的类群中，即新分配的响应变量总是接近所在组的形心（centroid）。交叉验证误差（CVRE）可以被用来评估回归树的预测能力。公示为：

图1. 交叉验证误差的公式

因此，CVRE可以定义为验证组未能被树解释的方差除以响应变量的总方差。当然，CVRE计算公式的分子会随着分组情况的变化而变化。理想的预测情况下，CVRE值为0是最理想的预测结果，CVRE的值越接近于1，预测结果越差。

mvpart程序包中MRT的运算流程

1. 将数据随机划分成k组（默认为10组，xval＝10）。

2. 从k组中随机选取一组作为“验证组”（testing set），剩余k-1组（训练组，training set）重现混合后通过约束分析，按照组内平方和最小的原则，建立回归树。

3. 将以上过程重复k-1次，即依次剔除一组数据。

4. 共产生k个回归树，对于每个回归树的不同分类方案，将验证组（一组数据）内的对象分配到分组结果中。计算每个回归树分类方案的CVRE。

5. 对回归树进行剪枝：可以保留CVRE最小的分类方案。也可以根据“1SE”准则，保留CVRE在最小的CVRE加1个标准差范围内最小的分类方案。

6. 为了获得上面运行过程的误差估计值，需要重复多次（100次或500次）将对象随机分配为k组。

7. 置换检验保留具有最小CVRE误差（或者是1SE值最小，此法最常用）的回归树。

实战：

后台回复“mvpart”获取安装包windows和mac版下载链接

#计算需要用到的主要函数为mvpart程序包的mvpart()。该函数要求输入的响应变量为matrix，解释变量为data.frame(). 
#安装程序包，因mvpart在R网站已不在更新，因此必须从本地安装；
# 后台回复“mvpart”获取安装包windows和mac版下载链接
#windos请mvpart_1.6-1.zip解压缩后拷贝mvpart文件夹到R程序安装文件夹library目录下，例如：D:\Program Files\R\R-3.1.1\library (找不到安装包，后台回复 mvpart 或 mrt 获得下载链接)
#for mac https://cran.r-project.org/src/contrib/Archive/mvpart/
#install.packages("~/training/MRT/mvpart_1.6-2.tar", repos = NULL,type="source")

#加载mvpart包
library(mvpart)
#调用程序包自带数据集spider
#spider有28行18列，前12列为不同种蜘蛛的多度数据，剩余的为环境数据
data(spider)    
#defaults,"."代表env中所有变量
fit<-mvpart(data.matrix(spider[,1:12])~herbs+reft+moss+sand+twigs+water,spider)

图2. 内置蜘蛛数据的多元回归树分析分组结果

从图中可以看出蜘蛛数据所在的样方被分成3组，首先被herb变量分成2组，然后第一组又被twigs变量分成了两组。图下边Error为误差，CV Error为交叉验证误差，SE为标准差

注意：因MRT交叉验证存在随机过程，每次分类可能存在不同结果，这时最好多运行几次或者选择交互模式，选择合适的分组方案。

#设定xv="1se"，根据“1SE”准则自动选择最优分类方案，与默认结果相同，因为默认选择1SE fit<-mvpart(data.matrix(spider[,1:12])~herbs+reft+moss+sand+twigs+water,spider,xv="1se")

图3. 根据“1SE”准则自动选择最优分类方案

#xv="min",选择具有最小CVRE值的回归树 fit<-mvpart(data.matrix(spider[,1:12])~herbs+reft+moss+sand+twigs+water,spider,xv="min")

图4. 选择具有最小CVRE值的回归树

从图中可以看出，数据被分成了7组，最小CVRE的树所分组数通常会比“1SE”准则选择的最优分类方案的组数多。分组越多，CVRE越小，从这里我们可以理解剪枝的概念，尽管此时，CVRE最小，但是我们有时候要兼顾分组的组数，有时候分组太多，并不利于我们对数据的分析。这也是函数默认选择“1SE”准则的原因。

#xv="pick",允许通过人机交互方式从函数提供的误差图中选择自己认为合适的分组 fit<-mvpart(data.matrix(spider[,1:12])~herbs+reft+moss+sand+twigs+water,spider,xv="pick") #运行函数生成误差图，绿色为相对误差，蓝色为CVRE，红色水平线指示最小CVRE（大红点）的1个标准差范围。红点为具有最小CVRE的分类方案。橙点为最接近1SE准则的分组。

图5. 人机交互方式从函数提供的误差图中选择自己认为合适的分组

通常，红点和橙点之间的分组方案都是可以接受的分组方案，我们用鼠标左键点击合适的分组处的点，就会生产多元回归树的树形图。例如，我们点击橙点，生成和前边xv=”1se”相同的分类结果。

图6. 点击合适的分组处的点，就会生产多元回归树的树形图

#还可以设置交叉验证的迭代次数，默认xvmult = 0，上方绿色条形图指出获得最佳分类方案交叉验证迭代的次数 fit<-mvpart(data.matrix(spider[,1:12])~herbs+reft+moss+sand+twigs+water,spider,xv="pick",xvmult =100)

图7. 绿色条形图指出获得最佳分类方案交叉验证迭代的次数

#根据“1SE”准则自动选择最优分类方案，以pca图的形式展示分组 fit <- mvpart(data.matrix(spider[,1:12])~herbs+reft+moss+sand+twigs+water,spider,xv="1se",pca=TRUE)

点击生产的图，将以pca图的形式展示分组

图8. 以pca图的形式展示分组

直观显示多元回归树应用的实例

#例子来源于数量生态学，并对代码中错误的地方进行了修改 #改变工作路径 setwd("D:/training/MRT") #导入实例数据 #法国和瑞士边境的Jura山脉的Doubs河30个取样点的数据，spe为各种鱼类在各取样点的多度，env为11个与河流相关的环境变量，spa为样方的地理坐标 #Species (community) data frame (fish abundances) spe <- read.csv("DoubsSpe.csv", row.names=1) #Environmental data frame env <- read.csv("DoubsEnv.csv", row.names=1) #Spatial data frame spa <- read.csv("DoubsSpa.csv", row.names=1) #"."代表env中所有变量,可以选择你感兴趣的变量也可以都选择 spe.ch.mvpart <- mvpart(data.matrix(spe) ~ ., env, xv="1se") summary(spe.ch.mvpart)

图9. 河流取样点分组。多元回归树展示分组情况，样品可以被分成3组，首先das（距离源头距离）变量将样品分为含18个样品和12个样品的两组，然后amm（铵浓度）变量又将12个样品分为2组。

#Group composition (labels of terminal nodes) where指示样品分组情况 (gr <- spe.ch.mvpart$where) #Renumber clusters sequentially gr2<-gr for (i in 1:length(unique(gr2))){  gr2[gr2==unique(gr2)[i]]<-i }     #Plot of the clusters on a map of the Doubs river     plot(spa, asp=1, type="n", main="MRT groups", xlab="x coordinate (km)", ylab="y coordinate (km)") lines(spa, col="light blue") text(50, 10, "Upstream", cex=1.2, col="red") text(25, 115, "Downstream", 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("Group",1:k), pch=(1:k)+20, col=2:(k+1), pt.bg=2:(k+1), pt.cex=2, bty="n")

图10. 分组结果按坐标展示，按组着色

与样品地理分布信息结合，可以清楚显示我们样品的分组情况，然后可结合地理生物学知识进行解释

参考文献

1. 数量生态学—R语言的应用，作者: 博卡德 (Daniel Borcard) / 吉莱 (Franqois Gillet) / 勒让德 (Pierre Legendre) ，赖江山译

2. 什么是多元回归树？作者：张金龙，http://blog.sciencenet.cn/home.php?mod=space&uid=255662&do=blog&id=417115

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：生信宝典 学术图表 高分文章 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外150+ PI，1500+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读