基于mRMRe的最大相关最小冗余特征选择

基于mRMRe的特征选择

参考：

1. 论文：De Jay N, Papillon-Cavanagh S, Olsen C, et al. mRMRe: an R package for parallelized mRMR ensemble feature selection[J]. Bioinformatics, 2013, 29(18): 2365-2368.
2. 代码：https://cran.r-project.org/web/packages/mRMRe/index.html

5.1.1 方法介绍

• 在过去的十年中，通过将不太准确的分类器组合在一起来产生高精度的“集合”分类器称为了一种新的趋势，其理由是在不影响偏差的情况下降低了模型方差(Kittler et al.， 1998)。 然而，mRMR方法的集成变体的计算成本很高，因为需要进行多次mRMR特征选择。 这是开发一个新的R包mRMRe的主要动机，它实现了mRMR的一个集成变体，其中构建了多个特征集，而不是一个单一的特征列表。 软件包中还包括一个函数，用于计算基于每个变量类型(即连续、离散和生存数据)的适当估计的互信息矩阵(mutual information matrix, MIM)。 这两个包特性都经过了调整，以充分利用多核平台。

• mim使用互信息(mutual information, MI)等基于关联的线性近似来计算互信息矩阵(mutual information matrix, MIM)。互信息通过如下式子定义：

$$I(x,y)=-\frac{1}{2}ln(1-\rho (x-y{)}^2)$$

其中$I$和$\rho$分别表示变量$x$和$y$之间的互信息和相关系数。连续变量之间的相关系数可以通过皮尔逊(Pearson’s)或者斯皮尔曼(Spearman’s)估计量来计算，而Cramer’s V可用于计算离散变量之间的相关性，Somers’ Dxy index可用于计算连续型变量与生存数据之间的关联。

• mRMR技术，通过mRMR.classic函数来实现，允许进行关联和非关联特征的有效选择(Ding and Peng, 2005)。令$X=x_1,...,x_n$为$n$个输入特征集合，$y$是输出变量。这个方法通过最大特征与$y$的互信息（最大相关）和最小化先前选择的特征变量的平均互信息（最小冗余）来对$X$排序。令$x_i$是与输出变量的互信息最高的特征，那么应该先选择：
  $$x_i=arg\underset{x_i\in X}{\max }I(x_i,y)$$
  选择的特征集合定义为$S$, 使用$x_i$来初始化。接着，通过选择与输出变量相关性最高的特征且与先前选择的特征冗余性最低的特征来选择特征，第$j$步的最大化打分$q$为：
  $$q_j=I(x_j,y)-\frac{1}{|S|}\sum _{x_k\in S}I(x_j,x_k)$$
  重复这个步骤，直到达到理想的求解长度。我们在连续性的或生存性的以及离散性的变量上可以实现这个方法，分别称之为F-test Correlation Difference (FCD)和Mutual Information Difference (MID) 框架，如表1所示。

• 尽管mRMR是一种快速的贪婪的启发式算法，但是其不保证找到可能存在的最优全局解。可能存在其他质量相当或者更好的特征子集。另外，单个mRMR运行选择的特征不太可能充分说明研究中与表型相关的生物过程的多样性。

• 为了解决这些问题，我们实现了两个集成方法来并行地产生多个mRMR解；这辆种方法分别称之为$exhaustive$(详尽的)和$bootstrap$(自助的)集成mRMR。这种详尽的变体通过使用$k>1$个最相关的特征来初始化多个特征选择过程，拓展了传统mRMR启发式算法。后续中，k mRMR求解流程是并行的，其中第一个选择的特征保证是不同的。实现广泛集成mRMR特征选择算法的伪代码在Supplementary Material的Algorithm1中提供。

• 自助的mRMR变体对原始数据进行有放回地重采样，来产生k个bootstraps, 然后再每个bootstrapped datasets中并行化地进行传统的mRMR特征选择，进而产生k个mRMR解。描述自助式集成mRMR特征选择算法的伪代码在Supplementary Material的Algorithm2中提供。

• 现有的mRMR包的相当大一部计算复杂性是由于互信息矩阵(mutual information matrix, MIM)的低效计算。例如minet(Meyer et al., 2008)和sideChannelAttack([Lerman et al., 2011](javascript:😉)包，在进行mRMR特征选择前要计算互信息矩阵。然而， 在特征选择过程中，只需要计算一小部分互信息矩阵来计算mRMR打分。我们通过lazy-evaluation的方式（只有当需要时才计算）来计算特征之间的互信息打分来加速mRMR的实现。这显著降低了mRMR的运行复杂度，是控制集成变体运算时间的一个关键点。

5.1.2 mRMRe的教程

参考： https://cran.r-project.org/web/packages/mRMRe/vignettes/mRMRe.pdf

代码位于： F:\Zhaolab2020\gut-brain-axis\metaAD\scripts\06predict\02mRMRe_feature_selection.R

#install.packages("mRMRe")
library(mRMRe)

set.thread.count(4)
data(cgps)
data.annot <- data.frame(cgps.annot)
View(cgps.ic50) ##
View(cgps.ge)
data.cgps <- data.frame(cgps.ic50, cgps.ge)  
View(data.cgps)
##data_cgps, samples/observations in rows and features/variables in columns

data_cgps如下图所示，第1列为ic50, 后面的若干列是基因。

## Test on the 'cgps' dataset, where the
## variables are all of continuous type.
dd <- mRMR.data(data = data.cgps)
dd <- subsetData(dd, 1:10, 1:10)
# Uses Spearman as correlation estimator
spearman_mim <- mim(dd, continuous_estimator = "spearman")
print(spearman_mim[1:4, 1:4])

# Uses Pearson as correlation estimator
pearson_mim <- mim(dd, continuous_estimator = "pearson")
print(pearson_mim[1:4, 1:4])

##<--------------------------2.2 Correlations--------------------------->##
# Compute c-index between feature 1 and 2
correlate(cgps.ge[, 1], cgps.ge[, 2], method = "cindex")
# Compute Cramer's V
x <- sample(factor(c("CAT_1", "CAT_2", "CAT_3"),
                   ordered = TRUE), 100, replace = TRUE)
y <- sample(factor(c("CAT_1", "CAT_2"),
                   ordered = TRUE), 100, replace = TRUE)
correlate(x, y, method = "cramersv")
# Compute Pearson coefficient with random strata and
# sample weights between features 1 and 2
strata <- sample(factor(c("STRATUM_1", "STRATUM_2",
                          "STRATUM_3"),
                        ordered = TRUE),
                 nrow(cgps.ge), replace = TRUE)
weights <- runif(nrow(cgps.ge))
correlate(cgps.ge[, 1], cgps.ge[, 2], strata = strata,
          weights = weights, method = "pearson")

• 计算相关性:

  correlate(X, Y, method = c("pearson", "spearman", "kendall", "frequency", "cramersv",
  		"cindex"),strata, weights, outX = TRUE, bootstrap_count = 0, alpha = 0.05,
   		alternative = c("two.sided", "less", "greater"))
  

##<-------------------------3.mRMR Feature Selection-------------------->##
## 3.1 Classic mRMR
dd <- mRMR.data(data = data.cgps)
dd_filtered <- mRMR.classic(data = dd, target_indices = c(1), feature_count = 30)   ##dd中包括200个样本，1000个特征。
###问题：如何知道提取的是哪些特征呢？

## 3.2 Ensemble mRMR
dd <- mRMR.data(data = data.cgps)
# For mRMR.classic-like results
mRMR.ensemble(data = dd, target_indices = c(1), solution_count = 1, feature_count = 30)
# For mRMR.ensemble-like results
mRMR.ensemble(data = dd, target_indices = c(1), solution_count = 5, feature_count = 30)