R-hier.part包的层次划分方法及重大bug
R-hier.part包的层次划分方法及重大bug
近年来，多元回归分析(如广义线性模型，GLMs)在生态学和保护生物学中得到了广泛的应用。然而解释变量之间的多重共线性(相关性)，会严重影响这种统计建模方法。
共线性带来的问题可通过在模型创建过程中去除高度相关的解释变量，或者使用主成分分析(PCAs)将PCA导出的因子作为预测变量来优化。
采用层次划分的分析方法可以有效地缓解共线性问题。层次划分(hierarchicalpartitioning，HP)是一种多元回归分析方法，它在解决多重共线性问题的同时，还能识别出最可能的影响因素。
HP通过确定每个解释变量对响应变量的独立贡献来减少共线性问题，并将其从与其他变量的相关性产生的联合贡献中分离出来。这允许在独立于其他协变量解释响应变量时对协变量的重要性进行排序。
由于其对多元回归分析的补充作用，它在生态和保护方面的应用正在增加。
hier.part这个包就是专门进行此分析的。对128篇使用该包的文献进行分析，其中大部分属于“生态学”主题，而在“生物多样性保护”和“环境科学”主题中也有部分论文。
该包于2004年首发，最近一次更新在2013年，版本为1.0_4。
图1 使用hier.part的论文数
hier.part的使用
#安装
>install.packages("hier.part")
>library(hier.part)
#最核心的函数为hier.part。用法：
#hier.part(y, xcan, family = "gaussian", gof = "RMSPE",barplot = TRUE)
#y: 因变量向量
#xcan: 包含n个独立变量的数据
#family: glm广义线性模型的参数。其选择有：
#binomial(link = "logit")
#gaussian(link = "identity")
#Gamma(link = "inverse")
#inverse.gaussian(link = "1/mu^2")
#poisson(link = "log")
#quasi(link = "identity", variance = "constant")
#quasibinomial(link = "logit")
#quasipoisson(link = "log")
#gof :拟合优度度量.默认RMSPE为均方根。其他还有logLik对数似然函数；Rsqu决定系数（R2)
#barplot: TRUE会对每个变量单独和综合解释方差的百分比画图
##结果包含3个内容：
#gfs：每个独立变量的组合情况；以及拟合优度的度量
#IJ：I为变量独立的贡献；J为变量综合的贡献
#I.perc：I在总解释方差中的百分比
##Example
#描述流域特征的七个自变量对河流中电导率的线性回归
>data(urbanwq)
>env <- urbanwq[,2:5];env
fimp sconn sdensep unsealden
1 0.5115957 0.7522154 7.266361 0.005821574
2 0.6835963 0.9453745 6.148170 0.000295345
3 0.3976425 0.4129818 7.842194 0.013312857
4 0.5258169 0.6886951 6.618157 0.012302183
5 0.4328468 0.3457947 7.905694 0.017589067
6 0.3860738 0.0000000 7.655064 0.010000018
7 0.5895874 0.8877404 8.966605 0.004386429
8 0.8290842 0.9909369 0.000000 0.000000000
9 0.1642358 0.0000000 0.000000 0.017359571
10 0.1667202 0.0000000 1.183216 0.004718112
11 0.4433829 0.1735454 9.391486 0.015060110
12 0.4788387 0.4442015 8.654479 0.007291771
13 0.5322383 0.2879531 11.874342 0.012958804
14 0.7892755 0.9934153 0.000000 0.000000000
15 0.5668667 0.7643291 6.920983 0.006726084
>hier.part(urbanwq$lec, env, fam = "gaussian", gof = "Rsqu")
$gfs
[1] 0.00000000 0.59047318 0.82682075 0.01098174 0.39635905 0.83410378 0.59758915
[8] 0.62061433 0.82683793 0.82691037 0.41798147 0.83447868 0.83411192 0.62061622
[15] 0.82691040 0.83451931
$IJ
I J Total
fimp 0.235822420 0.354650761 0.59047318
sconn 0.455960475 0.370860277 0.82682075
sdensep 0.005274981 0.005706756 0.01098174
unsealden 0.137461437 0.258897613 0.39635905
$I.perc
I
fimp 28.2584736
sconn 54.6374983
sdensep 0.6320981
unsealden 16.4719300
然而，2010年的一篇Plos one却发现了此包中存在的一个严重缺陷。
此研究基于该包的前两个较老的版本。作者发现当变量超过9个时，变量的顺序对预测因子解释的独立方差量有较大影响，即我们输入的变量顺序不同最后同一因子对方差的解释结果也不同。
理论上相同的一组变量在不同的顺序下应该产生相同的结果。然而本研究的结果发现对于包含9个以上变量的模型，90%以上的情况结果都不同。
目前还不清楚导致HP结果出现这种意外变化的原因。但是已经有20.3%的使用HP包的研究在模型中使用了9个以上的解释变量(图1)。因此一些研究的结果可能存在错误。
#随便增加两列，使得变量数目超过9个。
>envvv = cbind(urbanwq[,3:11],new1=1:15,new2=3:17)
>hier.part(urbanwq$fimp, envvv, fam = "gaussian", gof = "Rsqu")
$I.perc
I
sconn 22.4477628
sdensep 1.1590317
unsealden 12.1985380
fcarea 1.2249424
selev 9.0407159
amgeast 22.9262141
ldoc 7.7136738
lec 17.4076980
lnox 0.8542297
new1 3.0697276
new2 1.9574661
##调整一下顺序的结果：
>envvv = cbind(new1=1:15,new2=3:17,urbanwq[,3:11])
>hier.part(urbanwq$fimp, envvv, fam = "gaussian", gof = "Rsqu")
$I.perc
I
new1 1.564391
new2 2.237838
sconn 22.930206
sdensep 3.045801
unsealden 11.549721
fcarea 2.902707
selev 9.558250
amgeast 21.031273
ldoc 7.676860
lec 14.688723
lnox 2.814230发现没有，真的不一样了！！！
而且我尝试了好几次，拿不同的变量当做因变量向量y，最后结果不但不一样，而且并不是每次都会有warming信息：
Warningmessage:
hier.part produces a rounding error if number of variables >9
See documentation.
显然目前的版本仍没有解决这个问题。
虽然该包的作者声称对于包含9个以上自变量的分析，该函数会产生“较小的舍入误差”。但得到的结果却产生了一个相当不一致的定量和定性的结果。
所以在运行9个变量以上的数据时，一定要小心！！！
作者建议运用HP模型时最好不要超过9个变量。如果变量超过了9个，建议至少运行100次。
Reference
Olea, P. P., P. Mateo-Tomas and A. de Frutos (2010). "Estimating and Modelling Bias of the Hierarchical Partitioning Public-Domain Software: Implications in Environmental Management and Conservation." Plos One 5(7).
腾讯云开发者
