zhang2sha的个人博客分享 http://blog.sciencenet.cn/u/zhang2sha

博文

一种简单易行的方差分解方法

已有 14842 次阅读 2020-10-10 00:09 |系统分类:科研笔记| 方差分解

此前的博文中,给大家介绍了如何用几个圆圈来进行方差分解的方法,但事实上该方法依然有一定的复杂性,且几个圈交汇处的方差并不太容易解释。

那么有没有更简单一点的方法,能够方便的对比回归模型

Y=b0+b1*x1+b2*x2+…

中x1,  x2,  x3…对Y贡献的相对大小,以及不同X所属的因素类别(如生物因素,非生物因素)对Y的贡献大小呢?

当然有,本文就试图介绍这样一种方法,供大家在实际应用中参考。

要深入理解这种方法,首先要理解回归方程的本质意义:即代表了X变化一个单位时,Y会相应的变化多少。

理解了这一点后,我们引入标准化这个概念,所谓数据的标准化,不过是把原始数据的尺度进行转换,将其范围锁定在一个区间之内(如0-1之间)。通常我们用到的标准化是把某一变量的均值转变为0, SD转变为1. 比如在一个最简单的模型中:

Y=b0+b1*x1

如果我们对x1进行了标准化(mean=0,  SD=1),我们记为standard_x,那么这时模型变为:

Y=b0+b1*standard_x

这时,由于x1的变化而造成Y发生的变化幅度该如何衡量呢?我们可以直接用SD来表示这种变化,因为b0为常数,所以SD(b0)=0,  所以这一项不予考虑。那么Y的变异只来自于b1*standard_x这一项。则:

SD(Y) = SD(b1*standard_x) = b1*SD(standard_x)

前面我们提到了,standard_x是标准化后的x,  所以,

SD(standard_x) = 1

那么非常简单:上述公式变成了

SD(Y) = SD(b1*standard_x) = b1*SD(standard_x) = b1*1 = b1

这时候,X对Y的影响的大小就变成X变化一个SD,  对应的Y会变化多少, 所以这种方法消除了量纲,不同X的回归系数具有了直接可比性。比如以下模型中如果x1,  x2, x3…都是标准化后的值,

Y=b0+b1*x1+b2*x2+b3*x3+…

那么b1,  b2,  b3直接表示的就是各个x变化一个SD,  Y对应的变化多少。

所以这时候,我们便可以直接用b1,  b2, b3…的大小作为不同x对Y贡献大小的依据,且其值具备了直接可比性(消除了X的量纲),这就是依据这种方法进行方差分解的核心原理。到这里,大家体会到了对X进行标准化的妙处了吧(当然,标准化的妙处不止这一点)。

但是,还是要强调一下,这种方法的应用的前提是,不同的X之间的相关性必须较低,否则首先还是要处理共线性问题。

好了,原理讲完,通过一个案例来演示一下这种方法的具体应用。案例来源(Gross et al. 2017,  PNAS),其具体研究内容不做赘述,我们看他的问题:

分析24个变量对Y(logM5)的影响,并判定这24个变量对Y的影响的相对大小(即对R2)的贡献程度。

这24个变量是:LAT,  SINLONG,  COSLONG,  ARIDITY,   SLO,  SAND,  PH,  I(PH^2),  ELEVATION,  CWM_logSLA,  I(CWM_logSLA^2),    CWV_logSLA,  I(CWV_logSLA^2),  CWS_logSLA,  CWK_logSLA    I(CWK_logSLA^2),  CWM_logH,  I(CWM_logH^2),  CWV_logH    I(CWV_logH^2),  CWS_logH,  CWK_logH,  I(CWK_logH^2) ,  SR

而这24个变量, 又可以归为Rrichness,  Abiotic,  Geo等五大类。其最终结果如下图所示:

插图1.png

 案例来源(Gross et al. 2017,  PNAS)


那么这个结果是如何得到的呢?

首先我们判定这些X对Y的总体贡献是多少。

其实很简单,构造一个fullmodel,  包括所有的变量(当然要把共线性高的因素先去除掉)

mod12<-lm(logM5 ~ LAT + SINLONG + COSLONG +  

            ARIDITY + SLO + SAND + PH + I(PH^2) + ELEVATION+

            CWM_logSLA + I(CWM_logSLA^2)+ CWV_logSLA + I(CWV_logSLA^2) +  CWS_logSLA + CWK_logSLA + I(CWK_logSLA^2) +CWM_logH + I(CWM_logH^2)+ CWV_logH + I(CWV_logH^2) +  CWS_logH + CWK_logH + I(CWK_logH^2) +

            SR ,  data=datatotal)

此全模型的r2即为论文中fig4a(即上图)中报道的r2, 即0.66。

但是细心的读者会发现, 上图中x也不够24个, 其他的弄哪儿了?是我搞错了么?别着急, 我们先做一个模型选择, 把那些跟Y的关系比较弱的X通通删除(模型选择也是个技术活, 感兴趣的我们以后慢慢道来)。一定要注意, 这里的每个X都是标准化之后的X。

模型选择的核心代码如下:

########模型选择

dd12<-dredge(mod12,  subset = ~ LAT & SINLONG & COSLONG & ARIDITY & SLO & SAND & PH &SR & ELEVATION & dc(CWM_logSLA, I(CWM_logSLA^2)) & dc(CWV_logSLA, I(CWV_logSLA^2)) & dc(CWK_logSLA, I(CWK_logSLA^2)) & dc(CWM_logH, I(CWM_logH^2)) & dc(CWV_logH, I(CWV_logH^2)) & dc(CWK_logH, I(CWK_logH^2)),  options(na.action = "na.fail"))

 

######提取最优模型集, 以AIC值差值小于2为临界点

subset(dd12, delta<2)

求模型平均后的参数估计值,  即figure 4a右边的参数估计图

de12<-model.avg(dd12,  subset = delta < 2)

summary(de12)

即得到以下结果:

插图2.png

我们提取这些参数估计值, 进行对应的作图:

就得到了跟文献中一模一样的结果,  这里的X跟文章图中的X数量就相同了吧!

插图3.png插图4.png

              原图结果                                                                      重现结果

100%结果完全一样(细节大家自己调吧,  我要哄娃,  就不费这劲了)!


那么,  最核心的方差分解,  也就是下图这个结果如何得到呢?

插图5.png

首先,  R2=0.66这个值,  我们刚才已经通过fullmodel得到了,  剩下的就是如何分析五类因素各自的贡献了, 那么怎么做呢?

非常简单, 把以上模型平均后的结果中属于每一类的X的回归系数相加, 除以所有X的回归系数之和就可以了。

在R中,可以:

########求各类型参数贡献的百分比

d1<-summary(de12)

d2<-d1$coefficients

d3<-d2[1, ]

d4<-d3[c(2:22)]

d5<-abs(d4)

sum(d5)

如我们要求LAT + SINLONG + COSLONG+ELEVATION这四个因素一起在总R2中所占的比率(也就是图中的geo这一类因子), 那么:

插图6.png

即为geo这个因素在模型解释的总R2中所占的比重为6.6%,  

Richness的贡献为:

插图7.png

即richness这个因素在模型解释的总R2中所占的比重为5.2%。

所以用各类回归系数的和除以所有回归系数的和, 就得到了各类因素对的贡献比率, 而所有的这五类因素的贡献加到一起, 当然就是100%, 这个案例里, 他们一起贡献了Y的变异的66%。

 

这种方法简单易行, 应用也较为常见, 如以下三篇PNAS论文都采用了这一方法。

插图8.png

案例来源:(Le Provost et al. 2020,  PNAS)

插图9.png

案例来源(Sirami et al. 2019,  PNAS)

 

插图10.png

案例来源:(García-Palacios et al. 2018,  PNAS)

此外做个小广告, 鄙人将于今年10月底和11月底, 分别做一期meta分析和混合效应模型的线上培训课程(其他老师还有结构方程模型, R统计制图与地图制作, Invest模型等), 敬请关注, 详情请扫描以下二维码

插图11.png

参考文献:

1)     García-Palacios,  P.,  N. Gross,  J. Gaitán,  and F. T. Maestre. 2018. Climate mediates the biodiversity–ecosystem stability relationship globally. Proceedings of the National Academy of Sciences.

2)     Gross,  N.,  Y. L. Bagousse-Pinguet,  P. Liancourt,  M. Berdugo,  N. J. Gotelli,  and F. T. Maestre. 2017. Functional trait diversity maximizes ecosystem multifunctionality. Nature Ecology &Amp; Evolution 1:0132.

3)     Le Provost,  G.,  I. Badenhausser,  Y. Le Bagousse-Pinguet,  Y. Clough,  L. Henckel,  C. Violle,  V. Bretagnolle,  M. Roncoroni,  P. Manning,  and N. Gross. 2020. Land-use history impacts functional diversity across multiple trophic groups. Proceedings of the National Academy of Sciences 117:1573-1579.

4)     Sirami,  C.,  N. Gross,  A. B. Baillod,  C. Bertrand,  R. Carrié,  A. Hass,  L. Henckel,  P. Miguet,  C. Vuillot,  A. Alignier,  J. Girard,  P. Batáry,  Y. Clough,  C. Violle,  D. Giralt,  G. Bota,  I. Badenhausser,  G. Lefebvre,  B. Gauffre,  A. Vialatte,  F. Calatayud,  A. Gil-Tena,  L. Tischendorf,  S. Mitchell,  K. Lindsay,  R. Georges,  S. Hilaire,  J. Recasens,  X. O. Solé-Senan,  I. Robleño,  J. Bosch,  J. A. Barrientos,  A. Ricarte,  M. á. Marcos-Garcia,  J. Miñano,  R. Mathevet,  A. Gibon,  J. Baudry,  G. Balent,  B. Poulin,  F. Burel,  T. Tscharntke,  V. Bretagnolle,  G. Siriwardena,  A. Ouin,  L. Brotons,  J.-L. Martin,  and L. Fahrig. 2019. Increasing crop heterogeneity enhances multitrophic diversity across agricultural regions. Proceedings of the National Academy of Sciences 116:16442-16447.

 

关注以下公众号,发送“方差分解数据及代码”,可获取本文案例数据及详细代码

二傻统计公众号二维码.jpg



https://blog.sciencenet.cn/blog-3442043-1253775.html

上一篇:关于Meta分析,你需要了解的十个技术问题
下一篇:PPT分享:Meta分析的原理与要点
收藏 IP: 14.131.244.*| 热度|

3 杨卫东 檀成龙 黄河宁

该博文允许注册用户评论 请点击登录 评论 (3 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-16 19:12

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部