R语言一种无偏变量选择的多元统计方法
导读
今天给大家介绍的是一款施琳大神（现就读于查尔默斯大学，postdoc）于几年发表的一个R包-
MUVR(Multivariate methods with Unbiased Variable selection in R)，主要是为计算预测变量和相应变量之间的关系提供了一种新的算法，那么这里的响应变量既可以是连续变量用于回归分析，也可以是因子变量用于分类分析。另外一个我之所以介绍这个软件的原因是，它可以用于处理大规模变量数目大于观测数的数据格式，所以还是很友好的一个软件。因为大多数科研单位都不太可能做到那么大的样本量，除非是豪门！
文章介绍
-
这个方法施琳师姐也是以一作和通讯的身份发表在了Bioinformatics杂志上，让我等有机会拜读。文章标题为Variable selection and validation in multivariate modelling
-
再给个文章截图，有需要的可以自行下载（Open Access，DOI:10.1093/bioinformatics/bty710）
文章截图
方法原理简介
- 那么从技术分析的角度，MUVR采用的是一种
双重复交叉验证（repeated double cross validation, rdCV）的递归式的变量选择方法。MUVR包括了偏最小二乘（partial least squares, PLS）和随机森林（random forest, RF）模型。它有两种变量的选择方法：
-
minimal-optimal variables: 最少-最优变量，主要用于发现具有预测性的生物标记物；
-
all-relevant variables: 全部相关的变量，主要用于生物学呈现和机制研究。
- MUVR的工作原理见下图
MUVR工作流程图
包的安装和案例演示
1.包的安装
- 这里建议大家在安装了R软件的同时，安装一个Rstudio软件，这个软件非常友好，提供了很多可视化的界面。然后大家的包的下载，安装以及工作都在Rstudio这个软件上。
library(devtools)
install_git("https://gitlab.com/CarlBrunius/MUVR.git")
library(MUVR)
2.对于Y变量是分类变量的数据
首先要对输入的数据做个说明，这里的预测变量，暂且用大写的X表示，以及响应变量，用Y表示，两个输入数据必须是观测数要要一致（不论是数目还是ID都要相符合）。另外X的变量不能重复，要每一个都是唯一的。另外的一点就是数据的前处理这一块，因为随机森林模型RF是一种缩放不变性（尺度不变性）的方法，所以对数据的转换不敏感，因此不论你用log对数、sqrt平方根，还是中心化centering、标度化scaling可能结果都没有什么影响；但是，对于PLS而言，对数据的转化或者标度化是非常敏感的，那么这个包也提供了一个数据前处理的函数
preProcess()。
-
好的，那么对于响应变量是分类变量的情况，这里也提供了一个可供练习的数据。是2016年Buck发表[2]的一篇关于三个不同村庄采集的疟蚊的微生物16S测序数据，数据集名叫“mosquito”，那么在R对话框中输入data(“mosquito”)就会加载进来两个数据框，其中一个分类的
Yout包含响应变量（三个水平，如，村庄），共29个样品；另外一个表Xotu是一个数值矩阵，包括1,678个OTUs（列）和29个样品（行）。 -
正式开始数据的分析
# Classification example using the "mosquito" data
# 加载相关的R包
library(doParallel) # 并行处理的包
library(MUVR) # MUVR包
# 加载MUVR包的mosquito数据
data("mosquito")
# 检查每个分类的样品数
table(Yotu) # 主要是为了确定后面的参数nOuter的设定，这里一般设为分类中样本数最小的那一个
# 设定方法参数
nCore=detectCores()-1 # 运用到线程的设置
nRep=nCore # 重复数的设定
nOuter=8 # 外部交叉验证的数目
varRatio=0.8 # 每次迭代包含的变量数目的比例
method='RF' # 根据需要选择模型，RF或者PLS
# 对数据进行前处理
Xotu <- preProcess(Xotu,trans = "sqrt")
# 用doParallel 设定并行处理
cl=makeCluster(nCore)
registerDoParallel(cl)
# 执行模型
classModel = MUVR(X=Xotu, Y=Yotu, nRep=nRep, nOuter=nOuter, varRatio=varRatio, method=method)
这里的运行速度主要是由参数的设置以及数据量的大小决定的，但是用我的笔记本也很快就出来结果了，还是蛮快的
- 结果解读
# 检查模型的性能和输出
classModel$miss # 在min, mid, max三种模型下分类错误的个数，可以看到区别不大，尽管用最少的变量
# min mid max
# 4 3 3
classModel$nVar # 在每种模型下选择的预测biomarker
# # min mid max
# 13 15 17
cbind(Yotu, classModel$yClass) # 将真实的分类和预测的进行合并
# Yotu min mid max
# 1_ID1 VK3 VK3 VK3 VK3
# 2_ID2 VK3 VK3 VK3 VK3
# 3_ID3 VK3 VK3 VK3 VK3
# 4_ID4 VK3 VK3 VK3 VK3
# … … … … … …
plotVAL(classModel)
在不同变量选择条件下，响应变量错误分配图
plotMV(classModel, model='min') # Look at the model of choice: min, mid or max
当选择最少变量时分类预测的可能性
plotVIP(classModel, model='min')
当选择min模型时，其选择的预测变量的准确性的箱线图
getVIP(classModel, model='min') # 变量的排名，排名较低的则为预测性最好的
# order name rank
# OTU_28 1 OTU_28 2.488095
# OTU_4 2 OTU_4 4.363095
# OTU_400 3 OTU_400 74.005952
# OTU_3454 4 OTU_3454 145.297619
# OTU_243 5 OTU_243 145.607143
# OTU_1208 6 OTU_1208 356.476190
# OTU_178 7 OTU_178 357.869048
# OTU_114 8 OTU_114 432.699405
# OTU_39 9 OTU_39 495.101190
# OTU_16 10 OTU_16 704.422619
# OTU_26 11 OTU_26 705.285714
# OTU_21 12 OTU_21 774.464286
# OTU_134 13 OTU_134 983.464286
2.1 小结
这里通过选择了min模型而得到13个变量可用于很好的预测响应变量的分类，因此可以作为区分三个村庄蚊子的biomarker。
3.对于Y变量是连续变量的数据
由于参数的设定基本上和前面的一致，因此这里就不在赘述，方法上选择"PLS"就行。这里也提供了一个数据供练习，那么这个是代谢组学数据，用的是2015年Hanhineva K发表在Mol Nutr Food Research上的"Discovery of urinary biomarkers of whole grain rye intake in free-living subjects using nontargeted LC-MS metabolite profiling"[3]，看能否从尿液代谢物中找到摄入全谷物黑麦的生物标记物。
- 结果解读
library(doParallel) # 并行处理的包
library(MUVR) # MUVR包
# 加载数据
data("freelive")
# 设定参数，这里和前面基本一致，请参考前面的解释
nCore=detectCores()-1
nRep=nCore
nOuter=8
varRatio=0.8
method='PLS' # 这里模型用的是PLS模型
cl=makeCluster(nCore)
registerDoParallel(cl)
# 开始计算
regrModel = MUVR(X=XRVIP, Y=YR, ID=IDR, nRep=nRep, nOuter=nOuter, varRatio=varRatio, method=method)
# 结束并行计算
stopCluster(cl)
# 检查模型的性能和输出
regrModel$fitMetric # 这里可以看到min,mid,max三种变量选择模式下的模型拟合度R2，可以看到min的拟合度是最好的
# $R2
# [1] 0.7978157 0.8066973 0.8591222
# $Q2 # Q2代表模型的可解释性，就是用于预测新加入变量的能力，三种模式差不多，max最好
# [1] 0.6772411 0.6735702 0.6616667
regrModel$nComp # 三种模式需要取多少主成分才能达到前面的拟合度和解释度
# min mid max
# 3 3 4
regrModel$nVar # 三种模式选择的变量的格式
# min mid max
# 19 27 38
cbind(YR, regrModel$yPred) # 合并真实的响应变量和三种模式预测的相应变量
# YR min mid max
# 1_ID1 11.0666667 10.736887 10.736887 10.736887
# 2_ID100 12.1000000 11.920450 10.124938 12.222567
# 3_ID101 57.7333333 51.636043 58.287970 78.334162
# 4_ID102 10.8000000 24.057478 24.057478 24.057478
# 5_ID104 6.6333333 22.543692 19.820767 17.172377
# 6_ID106 9.8333333 47.834552 56.527561 52.002799
# ... ... ... ... ...
- 相关图表
plotMV(regrModel, model='min')
选择min时R2和Q2的情况
PLSfit=regrModel$Fit$plsFitMin # 提取PLS模型
biplotPLS(PLSfit, comps=1:2, xCol = YR, labPlSc = FALSE, labPlLo = FALSE)
Loading 图
plotVIP(regrModel, model='min') ## 将min模式选择的变量按照排名作图
变量的排名，排名较低的则为预测性最好的
getVIP(regrModel, model='min') # Extract most informative variables: Lower rank is better
# order name rank
# HPX283.9974.1.2967925 1 HPX283.9974.1.2967925 5.255102
# RNX138.0677.2.7376492 2 RNX138.0677.2.7376492 8.239796
# RNX247.0152.2.2259395 3 RNX247.0152.2.2259395 8.538265
# RNX262.0145.2.1115243 4 RNX262.0145.2.1115243 11.617347
3.1 小结
回归模型可以通过一定数量的informative的变量去预测原响应变量的值，可以作为一组潜在的biomarker，用于将来的疾病诊断和无创检查
总结
该软件包集合了预测变量预测响应变量是分类变量或者连续变量的方法。可以看到能否比较高效的完成日常的宏基因组或者代谢组数据用于寻找不同分类情况下biomarker作为一种指示的任务，也可以找到能够比较准确预测响应变量的一组潜在的biomarker代谢物，用于将来的疾病诊断以及无创筛查等都是非常好的途径。
参考
[1] 本文出处：Variable selection and validation in multivariate modelling
[2] Bacterial associations reveal spatial population dynamics in Anopheles gambiae mosquitoes
[3] Discovery of urinary biomarkers of whole grain rye intake in free-living subjects using nontargeted LC-MS metabolite profiling
