MRM中进行变量筛选
写在前面的2+1件事：
1.最近总有人加我好友称呼我的时候把我的姓写错。我的姓是雷厉风行的厉！厉行节约的厉！不明觉厉的厉！不是日历的历啊！
2.从2017.12.4开始至今，已经坚持写了一年多，目前公众号的关注人数每周都在以个位数稳定增长。虽然现在不做任何宣传与推广，全靠读者自发的分享，但还是不断地有新人在关注。我深感欣慰。
+1.本文因为一些原因需要重新发一下。之前的那篇删掉了。本文主体内容和之前相同，只是在最后多加了一点说明。
回归正题：
前几天有人在R——ecodist&MRM methods文章中提问：
如何在MRM中对变量进行forward-selection？
其实我之前对forward-selection是啥意思也不太懂。这篇文章也是在自己学习了相关概念之后的一点粗浅简介。权当班门弄斧抛砖引玉，欢迎讨论~
写之前以为本文会很短，结果写完发现内容很多，因此先列一个提纲。本文内容包括：
1.变量筛选(Variable Selection)
2. 两大类变量筛选的方法
3. R中实现变量筛选
4. 其他方法检验因子之间的相关性及贡献
5. MRM中如何实现变量筛选
变量筛选(Variable Selection)
变量筛选是在所有的变量中选出解释或预测能力最好的子集的过程。这样做的原因有以下几个：
根据奥卡姆剃刀原理(Occam’s Razor),越简单的模型就是越好的模型，因此需要去掉冗余的变量；
不必要的变量会给其他变量增加噪声，由这些变量贡献的自由度也会浪费；
太多变量做一件事会增加共线性；
节约时间和金钱，减少成本。
两大类变量筛选的方法
01Stepwise Procedures
Backward Elimination:
先将所有变量都放入模型中进行拟合，计算每个变量的p值；
移除得到最大p值的变量，并对模型进行重新拟合；
反复进行，直到所有的p值都小于你所设定的阈值(如0.05)，不再删除变量。
Forward Selection:
和Backward Elimination思想正好相反：在模型中分别放入每一个变量并分别计算p值；
挑出最小p值的变量放入模型，再计算放入其他变量之后的p值；
反复进行，直到剩下变量中最小的p值大于你所设定的阈值(如0.05)，不再加入变量。
Backward Elimination相较于Forward Selection的优势在于允许一些低贡献度的变量进入模型，而Forward Selection不能。且这些低贡献度的变量与其他变量组合之后可能贡献度变大。
Stepwise Regression：
这种方法是前两种方法的综合。每一次拟合的过程都可以添加或删除一个变量。但这种方法的缺点也很明显：
同时增加和删除一个变量可能会错过最优的模型；
删除不太重要的变量往往会增加剩余变量的显著性，这样就会高估剩余变量的重要性；
作为一种统计上的方法，这种对变量的筛选并不能直接解决实际的问题，模型的选择不能脱离实验目的。且被删除的变量仍可以与关注的因子相关。
1
使模型的Akaike Information Criterion (AIC)或者Bayes Information Criterion (BIC)尽可能的小。过程和p值的筛选相同，不再赘述。
2
Adjusted R2。R2 = 1- RSS /TSS。增加变量总会减小RSS，最终使得R2增加。因此这种方法并不是很好。
3
使预测的模型残差平方和(Predicted Residual Sum of Squares (PRESS))最小。
4
Mallow’s Cp Statistic：
Criterion-based的方法搜索范围更广，且搜索方式更好。推荐采用这种方法进行变量筛选。
Reference
http://www.biostat.jhsph.edu/~iruczins/teaching/jf/ch10.pdf
R中实现Variable Selection
1.最简单的方法就是根据定义，变量一个一个的增加或者删除。
2.现有工具。Adjusted R2和Mallow’s Cp Statistic在R中有多种方式进行计算：
leaps package
leaps()对于给定的变量进行穷举搜索，找到线性拟合效果最优的子集。
Usage：
leaps(x=, y=, wt=rep(1, NROW(x)), int=TRUE, method=c("Cp", "adjr2", "r2"), nbest=10,names=NULL, df=NROW(x), strictly.compatible=TRUE)
x:预测的矩阵
y:响应变量
wt:权重
int:模型中增加一个截距
method:计算Cp，Adjusted R2，R2。默认计算Cp
df:自由度
strictly.compatible=TRUE条件下x不能超过31列，即最多31个变量。
> library(leaps)
> x<-matrix(rnorm(100),ncol=4)
> y<-rnorm(25)
> leaps(x,y)
$which
1 2 3 4
1 TRUE FALSE FALSE FALSE
1 FALSE FALSE TRUE FALSE
1 FALSE FALSE FALSE TRUE
1 FALSE TRUE FALSE FALSE
2 TRUE FALSE TRUE FALSE
2 TRUE FALSE FALSE TRUE
2 TRUE TRUE FALSE FALSE
2 FALSE FALSE TRUE TRUE
2 FALSE TRUE TRUE FALSE
2 FALSE TRUE FALSE TRUE
3 TRUE FALSE TRUE TRUE
3 TRUE TRUE TRUE FALSE
3 TRUE TRUE FALSE TRUE
3 FALSE TRUE TRUE TRUE
4 TRUE TRUE TRUE TRUE
$label
[1] "(Intercept)" "1" "2" "3" "4"
$size
[1] 2 2 2 2 3 3 3 3 3 3 4 4 4 4 5
$Cp
[1] 0.7316891 0.9694651 1.3730445 1.6584606 1.0795534 2.6059937 2.6543183 2.8720554
[9] 2.9684635 3.3729559 3.0031894 3.0774385 4.5859802 4.8652879 5.0000000
结果中which的每一行为一个模型，穷举共15个。TRUE表示保留x中对应列的变量。
Cp即为15个模型的Cp统计量。
regsubsets
regsubsets是leaps包中另一个函数，其不仅可以进行穷举，还可以用forward,backward stepwise, sequential replacement的方法进行。用法类似，但功能更强大。
MASS package
stepAIC函数可以根据AIC值进行变量筛选。其参数direction可以选择both，forward，backward三种。
olsrr package
> model <- lm(mpg ~ disp + hp + wt + qsec, data = mtcars)
> k <- ols_step_all_possible(model);k
#结果给出Cp，Adjusted R2，R2
# A tibble: 15 x 6
Index N Predictors `R-Square` `Adj. R-Square` `Mallow's Cp`
* <int> <int> <chr> <dbl> <dbl> <dbl>
1 1 1 wt 0.753 0.745 12.5
2 2 1 disp 0.718 0.709 18.1
3 3 1 hp 0.602 0.589 37.1
4 4 1 qsec 0.175 0.148 107.
5 5 2 hp wt 0.827 0.815 2.37
6 6 2 wt qsec 0.826 0.814 2.43
7 7 2 disp wt 0.781 0.766 9.88
8 8 2 disp hp 0.748 0.731 15.2
9 9 2 disp qsec 0.722 0.702 19.6
10 10 2 hp qsec 0.637 0.612 33.5
11 11 3 hp wt qsec 0.835 0.817 3.06
12 12 3 disp hp wt 0.827 0.808 4.36
13 13 3 disp wt qsec 0.826 0.808 4.43
14 14 3 disp hp qsec 0.754 0.728 16.3
15 15 4 disp hp wt qsec 0.835 0.811 5
> plot(k)
> p <- ols_step_best_subset(model);p
Best Subsets Regression
------------------------------
Model Index Predictors
------------------------------
1 wt
2 hp wt
3 hp wt qsec
4 disp hp wt qsec
------------------------------
Subsets Regression Summary
-------------------------------------------------------------------------------------------------------------------------------
Adj. Pred
Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
-------------------------------------------------------------------------------------------------------------------------------
1 0.7528 0.7446 0.7087 12.4809 166.0294 74.2916 170.4266 9.8972 9.8572 0.3199 0.2801
2 0.8268 0.8148 0.7811 2.3690 156.6523 66.5755 162.5153 7.4314 7.3563 0.2402 0.2091
3 0.8348 0.8171 0.782 3.0617 157.1426 67.7238 164.4713 7.6140 7.4756 0.2461 0.2124
4 0.8351 0.8107 0.771 5.0000 159.0696 70.0408 167.8640 8.1810 7.9497 0.2644 0.2259
-------------------------------------------------------------------------------------------------------------------------------
AIC: Akaike Information Criteria
SBIC: Sawa's Bayesian Information Criteria
SBC: Schwarz Bayesian Criteria
MSEP: Estimated error of prediction, assuming multivariate normality
FPE: Final Prediction Error
HSP: Hocking's Sp
APC: Amemiya Prediction Criteria
> plot(p)
> model <- lm(y ~ ., data = surgical)
> k<- ols_step_forward_p(model)
Forward Selection Method
---------------------------
Candidate Terms:
1. bcs
2. pindex
3. enzyme_test
4. liver_test
5. age
6. gender
7. alc_mod
8. alc_heavy
We are selecting variables based on p value...
Variables Entered:
✔ liver_test
✔ alc_heavy
✔ enzyme_test
✔ pindex
✔ bcs
No more variables to be added.
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.884 RMSE 195.454
R-Squared 0.781 Coef. Var 27.839
Adj. R-Squared 0.758 MSE 38202.426
Pred R-Squared 0.700 MAE 137.656
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
ANOVA
-----------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
-----------------------------------------------------------------------
Regression 6535804.090 5 1307160.818 34.217 0.0000
Residual 1833716.447 48 38202.426
Total 8369520.537 53
-----------------------------------------------------------------------
Parameter Estimates
------------------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
------------------------------------------------------------------------------------------------
(Intercept) -1178.330 208.682 -5.647 0.000 -1597.914 -758.746
liver_test 58.064 40.144 0.156 1.446 0.155 -22.652 138.779
alc_heavy 317.848 71.634 0.314 4.437 0.000 173.818 461.878
enzyme_test 9.748 1.656 0.521 5.887 0.000 6.419 13.077
pindex 8.924 1.808 0.380 4.935 0.000 5.288 12.559
bcs 59.864 23.060 0.241 2.596 0.012 13.498 106.230
------------------------------------------------------------------------------------------------
> plot(k)
> model <- lm(y ~ ., data = surgical)
> k <- ols_step_backward_p(model)
Backward Elimination Method
---------------------------
Candidate Terms:
1 . bcs
2 . pindex
3 . enzyme_test
4 . liver_test
5 . age
6 . gender
7 . alc_mod
8 . alc_heavy
We are eliminating variables based on p value...
Variables Removed:
✖ alc_mod
✖ gender
✖ age
No more variables satisfy the condition of p value = 0.3
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.884 RMSE 195.454
R-Squared 0.781 Coef. Var 27.839
Adj. R-Squared 0.758 MSE 38202.426
Pred R-Squared 0.700 MAE 137.656
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
ANOVA
-----------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
-----------------------------------------------------------------------
Regression 6535804.090 5 1307160.818 34.217 0.0000
Residual 1833716.447 48 38202.426
Total 8369520.537 53
-----------------------------------------------------------------------
Parameter Estimates
------------------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
------------------------------------------------------------------------------------------------
(Intercept) -1178.330 208.682 -5.647 0.000 -1597.914 -758.746
bcs 59.864 23.060 0.241 2.596 0.012 13.498 106.230
pindex 8.924 1.808 0.380 4.935 0.000 5.288 12.559
enzyme_test 9.748 1.656 0.521 5.887 0.000 6.419 13.077
liver_test 58.064 40.144 0.156 1.446 0.155 -22.652 138.779
alc_heavy 317.848 71.634 0.314 4.437 0.000 173.818 461.878
------------------------------------------------------------------------------------------------
> plot(k)
> model <- lm(y ~ ., data = surgical)
> k <- ols_step_both_p(model)
>plot(k)
> model <- lm(y ~ ., data = surgical)
> k <- ols_step_forward_aic(model)
Forward Selection Method
------------------------
Candidate Terms:
1 . bcs
2 . pindex
3 . enzyme_test
4 . liver_test
5 . age
6 . gender
7 . alc_mod
8 . alc_heavy
Variables Entered:
✔ liver_test
✔ alc_heavy
✔ enzyme_test
✔ pindex
✔ bcs
No more variables to be added.
Selection Summary
----------------------------------------------------------------------------
Variable AIC Sum Sq RSS R-Sq Adj. R-Sq
----------------------------------------------------------------------------
liver_test 771.875 3804272.477 4565248.060 0.45454 0.44405
alc_heavy 761.439 4743349.776 3626170.761 0.56674 0.54975
enzyme_test 750.509 5515514.136 2854006.401 0.65900 0.63854
pindex 735.715 6278360.060 2091160.477 0.75015 0.72975
bcs 730.620 6535804.090 1833716.447 0.78091 0.75808
----------------------------------------------------------------------------
> plot(k)
Stepwise AIC Backward Regression 和 Stepwise AIC Regression类似，不再赘述。
Reference
https://cran.r-project.org/web/packages/olsrr/vignettes/variable_selection.html
https://stackoverflow.com/questions/21099840/variable-selection-methods
其他方法检验因子之间的相关性及贡献
要考察很多因子之间的相关性及他们对另一个因子的影响，方法目前也有很多，我之前很多文章都有提及。
效应量可以表征因子之间的相关性，详见效应量的计算——Cohen's d statistic
PCA，DCA，NMDS等降维方法也可以挑出其主要作用的成分，详见R-三种做PCA函数的差异：princomp，prcomp及rda
在计算CCA之前通常要计算膨胀系数(R中的vif.cca函数)来判断环境因子中的共线性，并删除高度共线性的因子。这部分以后可以简单写一下。
R包Hmisc中的varclus函数可对因子进行相关性聚类，找到共线性的因子，详见R package Hmisc——varclus函数简介
R包hier.part可用来计算多种因子对响应因子的贡献度，详见R-hier.part包的层次划分方法及重大bug
MRM中实现变量筛选
以上说了这么多，终于可以试着回答读者的提问了。需要注意的是上文第三部分提及的R中实现的变量筛选都是针对向量进行的，而MRM的输入数据是矩阵。
> library(ecodist)
>data(graze)
>head(graze)
sitelocation forestpct ACMI2 ANOD ASSY BRIN2 CIAR4 CIIN CIVU DAGL ELRE4
1.1.2001 12.187743 63.88 0 0 0.0 0 0 0.0 0.0 20.5 4.7
1.2.2001 12.186077 71.33 0 0 0.5 0 0 0.0 5.1 0.0 6.4
2.1.2001 12.406362 72.45 0 0 0.0 0 0 0.0 0.0 2.2 21.5
2.2.2001 12.416265 77.13 0 0 0.0 0 0 0.0 0.0 1.5 3.4
3.1.1998 8.409213 18.35 0 0 0.0 0 0 0.0 0.2 0.0 0.0
3.1.1999 8.409213 18.35 0 0 0.0 0 0 0.5 0.7 0.0 0.0
GAMO LOAR10 LOCO6 LOPE OXST PLMA2 POPR PRVU RAAC3 RUCR SORU2 STGR TAOF
1.1.2001 0 1.5 0 0 0.0 4.0 35.0 1.7 0.9 0.2 0 22.5 58.0
1.2.2001 0 8.9 0 0 0.0 0.0 54.5 1.0 0.5 0.0 0 10.8 39.0
2.1.2001 0 7.0 0 0 0.7 0.0 52.0 0.0 0.2 0.0 0 8.7 7.4
2.2.2001 0 14.5 0 32 0.2 0.0 11.0 0.0 0.0 0.0 0 11.4 5.7
3.1.1998 0 36.5 0 0 0.5 6.9 0.1 0.0 0.0 3.7 0 0.0 2.9
3.1.1999 0 71.5 0 4 0.0 3.6 0.0 0.0 0.0 1.4 0 0.0 4.4
TRPR2 TRRE3 VEOF2
1.1.2001 21.5 53.5 0
1.2.2001 3.7 39.5 0
2.1.2001 32.0 40.0 0
2.2.2001 5.2 72.5 0
3.1.1998 0.5 16.5 0
3.1.1999 0.0 4.8 0
>MRM(dist(LOAR10) ~ dist(sitelocation), data=graze, nperm=100)
$coef
dist(LOAR10) pval
Int 10.4452883 0.22
dist(sitelocation) -0.4174337 0.42
$r.squared
R2 pval
0.003978668 0.420000000
$F.test
F F.pval
4.885349 0.420000
结果中$coef里的pval是所有输入数据转化为矩阵的p值。
对于Backward Elimination思想，可将环境因子依次去掉一个做MRM，扔掉去掉因子后得到最小p值的因子，并重复这个过程。直到最大的p小于设定的阈值。
对于Forward Selection思想，可将环境因子依次与响应变量做MRM，保留最小的p值进入模型，并重复这个过程。直到最小的p大于设定的阈值。
手动挑选这种方法对因子比较少时适用。当因子多的时候还是要想点可以偷懒的策略。但是我在ecodist中也没有看到相关的函数。那就只能自己写代码了！
#Backward Elimination方法
>data(graze)
>test = graze[,-c(1:3)]
>head(test)
ANOD ASSY BRIN2 CIAR4 CIIN CIVU DAGL ELRE4 GAMO LOAR10 LOCO6 LOPE OXST PLMA2
1.1.2001 0 0.0 0 0 0.0 0.0 20.5 4.7 0 1.5 0 0 0.0 4.0
1.2.2001 0 0.5 0 0 0.0 5.1 0.0 6.4 0 8.9 0 0 0.0 0.0
2.1.2001 0 0.0 0 0 0.0 0.0 2.2 21.5 0 7.0 0 0 0.7 0.0
2.2.2001 0 0.0 0 0 0.0 0.0 1.5 3.4 0 14.5 0 32 0.2 0.0
3.1.1998 0 0.0 0 0 0.0 0.2 0.0 0.0 0 36.5 0 0 0.5 6.9
3.1.1999 0 0.0 0 0 0.5 0.7 0.0 0.0 0 71.5 0 4 0.0 3.6
POPR PRVU RAAC3 RUCR SORU2 STGR TAOF TRPR2 TRRE3 VEOF2
1.1.2001 35.0 1.7 0.9 0.2 0 22.5 58.0 21.5 53.5 0
1.2.2001 54.5 1.0 0.5 0.0 0 10.8 39.0 3.7 39.5 0
2.1.2001 52.0 0.0 0.2 0.0 0 8.7 7.4 32.0 40.0 0
2.2.2001 11.0 0.0 0.0 0.0 0 11.4 5.7 5.2 72.5 0
3.1.1998 0.1 0.0 0.0 3.7 0 0.0 2.9 0.5 16.5 0
3.1.1999 0.0 0.0 0.0 1.4 0 0.0 4.4 0.0 4.8 0
>response = test[,1:5] #前5列为响应变量
>predictor = test[,-c(1:5)] #剩余变量作为输入变量
>all = MRM(dist(response) ~ dist(predictor), nperm=10)
>origin.p = all$coef[2,2];origin.p # p = 0.276
[1] 0.276
#设定一个阈值，如p=0.2时终止
>p.set = 0.2
>p.last = origin.p
>while (p.last > p.set){
#每一列分别进行MRM，最大的p去掉后再重新对整体做MRM
p.each = c()
for (i in 1:ncol(predictor)) {
each = MRM(dist(response) ~ dist(predictor[,-i]), nperm=10)
p.each = cbind(p.each,each$coef[2,2])
}
p.each ##去掉最大p值的变量，再重新对剩下来的整体做MRM。这里有个小缺陷，如果同时有多个最大值，则一次全去掉了。
predictor = predictor[,-c(which(p.each == max(p.each)))]
all.new = MRM(dist(response) ~ dist(predictor), nperm=10)
p.last = all.new$coef[2,2]
print(ncol(predictor))
}
[1] 18
>p.last
[1] 0.185
>ncol(predictor)
[1] 18
Backward Elimination的方法只去掉一个因子就满足了条件。最终p为0.185。
###############
#Forward Selection
>data(graze)
>test = (graze[,-c(1:3)])
>head(test)
>response = test[,1:5] #前5列为响应变量
>predictor = test[,-c(1:5)] #剩余变量作为输入变量
#设定一个阈值，如p=0.2时终止
>p.set = 0.2
>p.last = 0.01
>predictor.new = matrix(data=NA,nrow=nrow(test))
>while (p.last < p.set){
#每一列分别进行MRM，取最小的p值加入predictor后重新对整体做MRM
p.each = c()
for (i in 1:ncol(predictor)) {
each = MRM(dist(response) ~ dist(cbind(predictor.new,predictor[,i])), nperm=10)
p.each = cbind(p.each,each$coef[2,2])
}
p.each ##取最小p值的变量，再重新对剩下来的整体做MRM。这里有个小缺陷，如果同时有多个最小值，则一次全加入模型。
predictor.new = cbind(predictor.new, predictor[,c(which(p.each == min(p.each)))] )
predictor = predictor[,-c(which(p.each == min(p.each)))]
all = MRM(dist(response) ~ dist(predictor.new), nperm=10)
p.last = all$coef[2,2]
print(ncol(predictor.new))
}
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
>p.last
[1] 0.216
>ncol(predictor.new)
[1] 8
Forward Selection进行了7次循环。最终变量为8个,p为0.216。
此结果再次证明了第二部分中，Backward Elimination更有优势的观点，其能够保留更多的变量。
这两种方法的代码在小细节上有不少差别，需要仔细看哦~
这两段代码很多地方还可以优化，如一次只严格的去掉或增加一个变量；for循环的优化等。以后用到了再做吧。
还有一个细节需要注意，如果nperm次数太少每次的结果都会差挺大的。nperm=1000每次结果也会有差别。
Plus:
最后的代码只是我自己在学习了基本概念之后进行的练习，并不一定完全正确，仅供参考。
Plus plus:
向别人提问是学习的一种正常途径。对问题的不同看法更可以让彼此多了解一种观点。如果有问题欢迎讨论。
一个环境工程专业却做生信分析的深井冰博士，深受拖延症的困扰。想给自己一点压力，争取能够不定期分享学到的生信小技能，亦或看文献过程中的一些笔记与小收获，记录生活中的杂七杂八。
欢迎大家扫描下方二维码关注我的公众号，若有问题也可直接加我的微信：水岸风堤(lii32703)。
欢迎分享，转载请联系我。