专栏首页小白鱼的生统笔记基于相似或相异度矩阵的多元回归(MRM)及R语言实例

基于相似或相异度矩阵的多元回归(MRM)及R语言实例

基于相似或相异度矩阵的多元回归(MRM)及R语言实例

基于相似或相异度矩阵的多元回归(Multiple regression on (dis)similarity matrices,MRM),是一种通过计算对象间相似或相异矩阵,代替原始变量组分进行回归的方法。在生物学数据分析中,使用频率非常高。

本篇通过一个非常实用的例子,介绍这种回归方法,及其在R语言中的计算过程。

文中涉及数据和R代码的百度盘链接(提取码,gnva):

https://pan.baidu.com/s/1J0Kcl4pf8QBSQELIGUAVdQ

若百度盘失效,也可在GitHub的备份中获取:

https://github.com/lyao222lll/sheng-xin-xiao-bai-yu

什么是MRM,及其特点

传统的基于变量的回归可能存在的局限

传统的多元回归基于独特的自变量建立与响应变量的关系模型。然而在生物学多变量数据分析中的一个难点是,许多数据集具有非常高维的特征,例如DNA序列数据、基因或蛋白表达谱、群落物种多度等,变量的数目远大于样本数量。在这些情况下,传统的多元回归程序将难以适用,很难在如此多的变量中建立响应关系。例如在确定试验条件对基因表达水平的影响时,如果拟合各个条件与每个基因表达值的关系,既繁琐又难以描述。再如,解释群落构建的成因时,若分别建立各环境变量与各物种丰度的回归,将也是很难解释效应。

另一方面,对于部分线性模型而言,例如多元线性回归通常基于F统计量估计模型的显著性,因此为了保证较高的统计功效,自变量要尽可能满足正态性以及方差齐性。然而大多数实际情况中,很难保证所有自变量都满足这些条件。尽管仍然可以继续使用这些自变量执行多元线性回归,只是回归中的参数估计和统计检验可能不再那么精确。

MRM的特点

一种解决方法是,根据多变量数据计算所有样本之间的相似度或距离,并构造N×N(N=样本的数量)的相似或相异度矩阵替代原始的多变量数据矩阵,它代表了各样本中一组变量组成差异的整体概况。随后,可通过基于相似或相异度矩阵的多元回归(Multiple regression on (dis)similarity matrices,MRM)在两个或多个相似或相异矩阵的基础上执行回归分析,以解释一组变量对另一组变量的整体效应,并使用非参数方法的置换检验来确定回归系数的显著性。

MRM因其具有较高的灵活性,广泛应用于解决高维的生物学数据问题。就以生物地理学研究中的一个常见问题为例,响应变量矩阵可以为通过物种多度数据计算出的相似或相异矩阵,自变量矩阵代表了预测或解释性数据(地理空间测量,环境因子组成等)计算出的相似或相异矩阵。通过相似或相异矩阵代替原始的多元数据,建立自变量和响应变量的多元回归,评估地理或环境因素对群落组建的效应,解释群落构建的成因,以及寻找显著贡献的地理或环境组合。

MRM使用时的一些注意事项

(1)必须为每组待分析的数据类型(自变量或响应变量矩阵)选择一个合适的相似或相异指数,以获得对象间的相似或相异测度。例如上述提到的生物地理学有关分析中,物种组成通常计算Bray-curtis距离测度或相似度,环境数据通常在标准化后计算欧式距离,地理位置通常根据经纬度坐标计算地理距离。

(2)多元回归分析不需要一定线性的,允许使用非参数或非线性多元回归方法。但是在方法选择时,注意选择方法的前提假设。

(3)不同对象之间的极端相似或相异测度(即接近0或1)可能严重影响置信区间估计。

(4)数据转换后会改变相似或相异矩阵中的值,进而可能会对回归系数产生很大影响,因此涉及到数据转换时要慎重选择合适的转换方法。

MRM应用的一个具体案例解读

我们不妨来看一下具体的MRM分析应用的实例吧。白鱼同学找了一个很有代表性的例子,一个涉及空间尺度的微生物地理学研究,节选了文章中的部分内容描述以帮助大家理解。

来源文献:Wang X, Lu X, Yao J, et al. Habitat-specific patterns and drivers of bacterial β-diversity in China’s drylands. The ISME Journal, 2017, 11(6): 1345-1358.

细菌群落组成相似度与地理距离和环境异质性有关

作者研究了中国北方干旱地区不同生境之间土壤细菌的β多样性。通过初步分析细菌群落组成相似度的距离衰减模式,发现在空间尺度上,细菌群落的组成相似度随着地理距离的增加而显著下降。但通过线性回归估算的距离衰减斜率在不同的生境中存在不同,表明潜在的驱动因素不同。Mantel分析表明,细菌群落组成的相似度与地理距离和环境异质性密切相关。

左上图,样本采样地点。

右上图,细菌群落组成相似度在不同空间尺度上的距离衰减模式分析的部分截图。横轴为采样地点间的地理距离,已经过ln(x+1)转化;纵轴为细菌群落组成相似度,计算为1-Bray-curtis。

表1,Mantel分析细菌群落组成的相似度与地理距离和环境异质性的关系。

MRM分析确定地理距离和环境因素对细菌群落组成相似度的效应

随后,作者为了进一步确定地理距离或环境因素对细菌群落组成相似度的相对贡献,使用了基于矩阵的多元回归(MRM),具体执行过程如下。

首先由于使用MRM模型的目的重在解释变量的效应,因此在回归之前,作者为了避免多重共线性的影响,使用Hmisc包varclus()函数评估环境变量的冗余度。作者发现,总氮(TN)、总有机碳(TOC)、氮磷比(NP)、年平均温度(MAT)等相互之间具有非常高的相关强度(Spearman ρ2 > 0.7),因此考虑将它们去除,只保留一个年平均降水量(MAP);ANPP和PR之间选择保留PR;其余变量间的相关性较低,可以忽略。在手动识别并过滤了高度相关的变量后,将剩余变量用作MRM分析。

MRM使用ecodist包中的MRM()执行。分析前,地理距离通过采样地点的经纬度坐标换算,并进行ln(x+1)转换;细菌群落的组成相似度计算为1-Bray-curtis距离,并进行ln转换;环境变量计算欧式距离。MRM基于多元线性回归,将所有自变量矩阵标准化后(这样可以使所有自变量的回归系数处于可比较的水平下,以评估相对重要性),分两次执行MRM:第一次使用全部的自变量矩阵拟合模型,并从中去除不显著的自变量后,将剩余的显著自变量用于拟合第二次的回归(这类似手动执行了一种后向选择)。

结果表格中展示了第二次回归的结果,因为第二次的回归中排除了不显著变量的影响。在高寒草原(Alpine grassland),MRM模型可解释高达75%的群落相似度变异(P <0.001),土壤水分(Soil moisture)和植物物种丰度(Plant richness)是解释群落相似度的最重要变量。MRM模型仅解释了沙漠(Desert)中细菌群落相似度变异的5%(P<0.001),地理距离是唯一的影响因素。在典型的草原(Typical grassland)中,土壤pH(Soil pH)和总磷(Soil total P)贡献更大,归因于它们更高的偏回归系数,而其它因素例如地理距离、年平均降水量(MAP)、海拔(Altitude)和植物物种丰度(Plant richness)则较小。总的来说,对于几乎所有的生境而言,年平均降水量(MAT)和地理距离是对细菌群落组成相似度影响最大的因素。

R语言的MRM分析过程展示

接下来,我们就展示MRM在R中的分析过程。白鱼同学觉得上文的示例非常好,整个过程中涉及了多重共线性问题的解决、多种相似度或距离矩阵的计算、数据转化、MRM分析以及可能用到的变量选择,总之考虑的非常全面。因此,就找了个类似的数据,模仿它的思路执行MRM。

备注:因为很多文献(包括上文示例)的数据都未公开,要么未提供物种组成数据,要么未提供环境因子数据…….因此,同时找到提供有物种多度数据(OTU丰度表)、环境因子数据以及地理位置数据共3个完整矩阵的合适例子,白鱼同学花了不少心思……

示例数据概要

数据来源文献:Milici M, Tomasch J, Wosoxley M L, et al. Bacterioplankton Biogeography of the Atlantic Ocean: A Case Study of the Distance-Decay Relationship. Frontiers in Microbiology, 2016: 590-590.

作者在大西洋表层深度(20-200m)采集水样,采样地点横跨南北半球覆盖了约12000 km的距离,并对浮游细菌进行生物地理学相关研究。文章中的部分内容提到,大西洋表层浮游细菌群落的相似度显示出明显的距离衰减模式,但是对于北大西洋、南大西洋、水域不同深度以及不同类群的浮游细菌等而言,群落组成相似度的潜在驱动因素不同。

左图,样本采样地点。

中图,三种不同的浮游细菌群落(FL、SPA和LPA,文中通过滤膜过滤按其体积大小划分),在北半球大西洋、南半球大西洋以及同时考虑南北南球大西洋时的群落相似度的距离衰减模式。

右图,FL浮游细菌群落,在北半球大西洋、南半球大西洋以及同时考虑南北南球大西洋时,不同水域深度中的群落相似度的距离衰减模式。

因此,接下来我们使用该文章中的数据,模仿上文第一个示例文献中的方法,尝试通过MRM探索影响浮游细菌群落组成变化的因素,包括地理位置、环境变量等。

备注:其实在原文中,作者使用了一种称为基于距离的多元线性模型(Distance-based multivariate analysis for a linear model, distLM),量化了地理位置、环境变量等对浮游细菌群落组成的相对贡献。distLM的过程和MRM是相似的,将变量转换为距离矩阵后,执行多元线性回归。文献中看到的distLM分析基本都是用PRIMER软件来做的,但这个软件买不起……所以也就不重现原文的distLM的过程了,而套用这个数据演示MEM。好了多的别问了,白鱼同学自己也没搞明白distLM和MRM的联系或区别......

然后回过头来继续看数据具体内容。浮游细菌OTU丰度表、地理以及环境数据等均可在原文附录中找到。下文为了方便演示,白鱼同学只从中节选了一部分来自北半球大西洋的物种和环境数据,可在网盘附件中找到。

“FL.otu_table.txt”为FL(文中对一类小型浮游细菌类群的统称)浮游细菌OTU丰度表。每一行为一种微生物OTU,每一列代表一个样本,列名称即采样地的名称。

“site_env.txt“中记录了采样地点的地理位置和水体环境。第一列是样本名称;Lat是纬度(正值为北纬N,负值南纬S,由于节选的数据均来自北半球大西洋,所以均为正值),Lon是经度(正值为东经E,负值为西经W);Depth是采样水体深度(m);Temperature为水体温度(℃);Salinity为水体盐度(psu);Chlorophyll_a为测量的水体中叶绿素a含量(µg/L),其代表了光合微型真核生物含量。

环境变量的多重共线性分析

由于本次执行回归分析的目的是了解地理或环境因素等对FL浮游细菌群落组成的影响,也就是作为“解释模型”而非“预测模型”来使用,因此需要避免变量间的多重共线性。多元回归中的多重共线性问题已经在前文有过简介,并提到了一些评估多重共线性的方法。在这里,仿照上文第一篇文献中的方法,使用Hmisc包varclus()函数评估环境变量的冗余度。

#读取地理坐标和环境组成数据
site_env <- read.delim('site_env.txt', sep = '\t', row.names = 1, check.names = FALSE)
 
#除了经纬度外,其它均定义为环境因素
env <- site_env[c('Depth', 'Temperature', 'Salinity', 'Chlorophyll_a')]
 
#评估环境变量间的多重共线性
#Hmisc 包 varclus() 的方法,计算变量间的 spearman 相关性进行评估
library(Hmisc)
 
env_varclus <- varclus(as.matrix(env), imilarity = 'spearman', method = 'complete')
plot(env_varclus)

结果显示,采样水域深度(Depth)和水域中叶绿素a含量(Chlorophyll_a)之间共线性程度较高(Spearman ρ2 接近0.5),暗示二者中保留一个就可以了。二者共线性的原因可能归因于,太阳光照的穿透强度有限,因此随着海洋深度的增加,水体中的光照强度会越来越弱,不利于光合生物的生长,因此在水体中测量的叶绿素a含量会随深度的增加而减弱。

其余变量间的共线性程度较低,可以忽略。因此,最终选择采样水域深度(Depth)、温度(Temperature)和盐度(Salinity)共三个环境因素用于后续建模。

备注:这里通过Spearman ρ2评估共线性程度的高低,最终结果的选择存在一定的主观性。例如我这里认为Depth和Chlorophyll_a的共线性程度较高,但如果您认为Spearman ρ2不超过0.5就不代表多重共线性明显。将二者都保留也是可以的。

相似度或距离矩阵的计算

为了执行MRM,需要基于变量组成计算样本间的相似度或相异度矩阵。

仿照上文,FL浮游细菌群落的组成相似度计算为1-Bray-curtis距离,并进行ln转换;地理距离通过采样地点的经纬度坐标换算,并进行ln(x+1)转换;环境变量计算欧式距离。

#读取浮游细菌群落物种组成数据
spe <- read.delim('FL.otu_table.txt', sep = '\t', row.names = 1, check.names = FALSE)
spe <- data.frame(t(spe))
 
#vegan 包 vegdist() 计算群落间物种组成的 Bray-curtis 相异度
#并通过 1-Bray-curtis 转换为相似度后,再通过 ln 转化
library(vegan)
 
comm_dis <- vegdist(spe, method = 'bray')
comm_sim <- log(1 - comm_dis, exp(1))
 
#地理坐标和环境组成数据已经在上文读取进来了,继续执行
 
#geosphere 包 distm() 根据经纬度计算采样点间的地理距离
#默认距离单位是米,除以 1000 可转换为千米单位,并再通过 ln(x+1) 转化
library(geosphere)
 
site_dis <- distm(site_env[c('Lon', 'Lat')]) / 1000  #distm() 要求两列数据,第一列是经度,第二列是纬度
site_dis <- log(site_dis+1, exp(1))
rownames(site_dis) <- rownames(site_env)
colnames(site_dis) <- rownames(site_env)
site_dis <- as.dist(site_dis)
 
#对于各环境变量,均计算为欧氏距离测度
depth_dis <- vegdist(site_env['Depth'], method = 'euclidean')  #深度
temperature_dis <- vegdist(site_env['Temperature'], method = 'euclidean')  #温度
salinity_dis <- vegdist(site_env['Salinity'], method = 'euclidean')  #盐度

MRM分析和结果解读

最后执行MRM,了解地理或环境因素等对FL浮游细菌群落组成的影响。类似地,考虑到线性关系最容易理解和解读,我们基于地理和环境的相异矩阵执行一个多元线性回归模型。

首先对所有相似度或距离测度进行标准化,以便在后续的回归中,能够使所有自变量的回归系数能够处于可比较的水平下,以评估相对重要性。MRM使用ecodist包中的MRM()执行。

#由于使用地理或环境的距离测度作为自变量拟合多元线性回归
#如果期望后续能够通过回归系数比较各自变量的相对重要性,需要统一对它们进行标准化处理
#MuMIn 包 stdize() 可以实现对地理或环境距离测度的标准化功能,参考自:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4961435/
library(MuMIn)
 
site_dis <- stdize(site_dis)
depth_dis <- stdize(depth_dis)
temperature_dis <- stdize(temperature_dis)
salinity_dis <- stdize(salinity_dis)
 
##MRM 分析
library(ecodist)
 
#详情 ?MRM,以下基于 9999 次置换估计 p 值
 
#综合多个自变量矩阵的多元线性回归,在多元回归中,各自变量的回归系数即“偏回归系数”
#“偏回归系数”的含义,对于某个自变量而言,当控制其它自变量保持恒定时,该自变量对响应变量的相对效应
fit_MRM <- MRM(comm_sim~site_dis+depth_dis+temperature_dis+salinity_dis,
    nperm = 10000, method = 'linear')
fit_MRM
 
#如果希望获得各个自变量独立的回归系数,则单独执行各自变量与响应变量的一元回归即可
#例如考虑地理距离对群落组成相似度独立的效应
fit_MRM_site <- MRM(comm_sim~site_dis, nperm = 10000, method = 'linear')
fit_MRM_site

这种情况下,通常不怎么关注变量独立的效应(一元回归关系),更多关注变量组合的整体效应及其相对效应(多元回归关系)。因此,下文中的“回归系数”一词均是指的多元回归中的回归系数(偏回归系数)。

基于相似度或距离矩阵的多元线性回归中,如果存在不显著的变量矩阵时,建议将它们手动剔除,之后再使用剩余显著变量矩阵拟合第二次回归,以排除不显著变量的影响,获得一个易于解释但不失准确度的模型。这也是上文第一篇文献中为什么执行两次MRM的原因(这类似手动执行了一种后向选择)。在这个结果中,由于各变量矩阵都是都是显著的,因此不再执行第二次回归。

这个基于相似度或距离矩阵的多元线性回归中,所有回归系数均为负值,且都是显著的,表明FL细菌群落组成相似度随着地理距离、环境异质性的增加而呈现显著线性降低的趋势。不过,MRM模型整体解释了约25%的群落相似度变异,不是很高,可能归因于尚有其它因素未考虑在内,或者复杂的生物学关系难以通过线性模型进行解释。

通过比较回归系数的绝对值大小,可以确定地理或各环境因素对群落组成相似度的相对重要性。其中,水域深度(Depth)的效应最大,盐度(Salinity)次之,再次为地理距离,温度(Temperature)效应更小。

因此,对于北半球大西洋表层200m深度内的FL浮游细菌群落而言,环境选择过程占据主要因素,而地理因素的扩散限制相对微弱。这也是很容易理解的,海洋中地理扩散限制的影响较小,可解释于海洋的动态非常剧烈,使得海洋微生物群落中很难定义真正的地方性物种。例如有研究表明,每天每m2存在数百万种微生物的移动,并且可以在大约4天的时间内传播超过11000km。通过比较环境因素的相对重要性,深度的影响最大,体现了海洋中不同深度的分层效应明显。

已经获得了各变量的相对重要性,最后不妨根据它们的重要性排名,手动组建一个前向逐步回归模型,以评估这些变量矩阵(环境或地理因素)对FL浮游细菌群落组成相似度影响的累积效应。累积解释率以R2作为呈现。同时累积R2的越趋平缓的递增趋势也侧面反映了变量的相对重要性。

#根据重要性排名,水深>盐度>地理距离>温度
#手动组建一个前向逐步回归,评估重要环境或地理因素对群落组成相似度的累积解释率
fit_MRM1 <- MRM(comm_sim~depth_dis, nperm = 10000, method = 'linear')
fit_MRM2 <- MRM(comm_sim~depth_dis+salinity_dis, nperm = 10000, method = 'linear')
fit_MRM3 <- MRM(comm_sim~depth_dis+salinity_dis+site_dis, nperm = 10000, method = 'linear')
fit_MRM4 <- MRM(comm_sim~depth_dis+salinity_dis+site_dis+temperature_dis, nperm = 10000, method = 'linear')
 
fit_MRM1$r.squared  #水深的贡献
fit_MRM2$r.squared  #水深和盐度的综合贡献
fit_MRM3$r.squared  #水深、盐度和地理距离的综合贡献
fit_MRM4$r.squared  #水深、盐度、地理距离和温度的综合贡献

本文分享自微信公众号 - 小白鱼的生统笔记(gh_5f751e893315),作者:生信小白鱼 鲤小白

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-06-20

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 一个简单线性回归和多项式回归在R中的实现示例

    常见线性回归的原理就不多说了,大家都懂,就是普通最小二乘法(Ordinary Least Square,OLS)回归法,主要用于简单线性回归、多项式回归和多元线...

    用户7585161
  • R语言绘制雷达图的示例

    雷达图(radar charts)又叫蜘蛛网图。传统的雷达图被认为是一种表现多维(4 维以上)数据的图表。它将多个维度的数据量映射到坐标轴上,这些坐标轴起始于同...

    用户7585161
  • R包geosphere根据经纬度坐标计算地理距离

    其实正文大部分内容和标题无关......毕竟根据样地的经纬度坐标计算地理距离,这个方法也没什么可讲的,在R中就一条函数的事......因此白鱼同学额外引用了一篇...

    用户7585161
  • 剖析Web技术栈(二)

    TCP/IP是一种使用套接字(socket)的网络协议。Socket是包括IP地址(在网络中是唯一的)和端口(对于特定的IP地址是唯一的)的元组,计算机使用IP...

    老齐
  • POJ1258:Agri-Net-最小生成树

    Farmer John has been elected mayor of his town! One of his campaign promises was...

    ACM算法日常
  • 华为首次曝光量子计算成果:全幅模拟42量子比特电路并提供云服务

    全联接大会最后一天(12日),华为公布了在量子计算领域的最新进展:量子计算模拟器HiQ云服务平台问世,平台包括HiQ量子计算模拟器与基于模拟器开发的HiQ量子编...

    新智元
  • 05.记录合并&字段合并&字段匹配1.记录合并2.字段合并3.字段匹配3.1 默认只保留连接上的部分3.2 使用左连接3.3 使用右连接3.4 保留左右表所有数据行

    将两个结构相同的数据框合并成一个数据框。 函数concat([dataFrame1, dataFrame2, ...])

    用户1250179
  • MySQL使用on duplicate key update时导致主键不连续自增

    在做数据统计的时候,我们经常会用到mysql的on duplicate key update语法来自动更新数据,比如

    JouyPub
  • 2018 Wannafly summer camp Day2--New Game!

    New Game! 描述 题目描述: Eagle Jump公司正在开发一款新的游戏。泷本一二三作为其员工,获得了提前试玩的机会。现在她正在试图通过一...

    Enterprise_
  • jQuery 动态绑定

    这是在项目过程中所遇到的一个问题,给 JS 动态生成的元素绑定事件失效,代码如下所示:

    Nian糕

扫码关注云+社区

领取腾讯云代金券