||
1)首先必须明确,多因素对Y的共同影响,并不一定等于单因素对Y的影响的作用之和。事实上,只有在不同的X彼此独立,即X之间完全没有相关性的前提下,多个X对Y的总影响才等于单个X对Y的影响之和。
比如,现在有一个模型
M0 <- lm(y~x1+x2+x3, data=d1);
假定其R2=30%, 那么30%是否就等于以下三个模型的R2之和呢?
M1 <- lm(y~x1, data=d1)
M2 <- lm(y~x2, data=d1)
M3 <- lm(y~x3, data=d1)
答案是,不一定,只有在x1, x2, x3这三个因素彼此完全独立,相关性为0的情况下,才成立。但通常情况下,这个假设都不成立,因为x1, x2, x3彼此之间或多或少,都存在一定的相关性。
2)上面提到,X之间或多或少,通常都会有一定的相关性,但这种相关性如果太强,即某些X彼此高度关联,这时候,模型对某因素效应的判定就会出现统计上的不可识别性。即对于高度相关的X(如X1, X2), 模型会无法区分二者各自的效应。通常认为,两个自变量的相关系数如果高于0.7,则二者不适宜存在于同一个模型中,所以对于相关系数高于0.7的变量,通常只在模型中保留一个即可(这时,也可以首先求自变量的主成分,然后拿主成分对Y进行拟合,但这稍显复杂,结果的直观性也较差,所以这并非上策)。
假设我们要用x1, x2两个变量来解其各自对Y的影响大小,可以用长方形表示Y的总体变异,也就是总方差VT, 而X1, X2对Y的影响,则由三部分构成:
A: X1单独对Y造成的影响
B: X2单独对Y造成的影响
C: X1, X2的协同效应对Y造成的影响
如前面所论述,如果,X1, X2,彼此完全独立,没有相关性,那么C=0,A, B完全分离。
图1中,A, B, C三部分可以由以下三个回归方程求得:
A+B+C, m1<-lm(Y~X1+X2)的R2
A+C, m2<-lm(Y~X1) 的R2
B+C, m3<- lm(Y~X2) 的R2
首先求C, C=(A+C)+(B+C)-(A+B+C)
然后求B, B=(B+C)-C
最后求A, A=(A+C)-C
当然,以上各项也可以直接通过另一种方法求得,比如我们要求A,
那么首先,求方程:
m3<-lm(Y~X2) 的R2,
然后用m1的R2 减去m3的R2 ,就是A
同样原理假设你有三个因素,如,环境,资源和微气候,对Y都有影响,那么和区别其各自的作用呢?
如图2案例(Heikkinen et al. 2005):
首先三个因素对Y的总的解释量是R2=1-41.4%=58.6%, 也就是模型
M1<-lm(y~habitat + resources +microclimate)的R2=58.6%,
那么接下来,比如我们要求a,也就是habitat单独对Y的影响,从原理上,这个影响就是我们同时控制了resources 和microclimate之后,habitat对Y造成的影响。所以我们首先构造一个模型
M2<-lm(y~ resources +microclimate)
然后求出该模型的R2
那么用M1的R2减去M2的R2就是a,
同样道理,我们可以依次求出,b, c.
a, b, c求出之后,我们便可求,d, e, f, g的值,
比如,求d, 首先求模型
M4<-lm(y~ microclimate) 的R2
然后,以模型M1的R2减去M4的R2,其值即a+b+d,
因为a, b已知,则:
d=(a+b+d)-(a+b),
同理可以求出e, f,
最后,求g, 因为我们已知模型
M6<-lm(y ~ habitat) 的R2 =H= a+d+g+e,
所以, d+g+e=H-a,
基于前面的运算,d, e已知,
则: g=(d+g+e)-(d+e)
这种方法目前应用较为常见,比如北大唐志尧老师2018年发表在PNAS上的论文,即采用这种方差分解方法(Tang et al. 2018)。如下图所示:
同时大家也能感到,这种方差分解方法计算比较繁琐(与此原理相同,r中vegan包已具有专门针对群落数据的方差分解方法,可以直接分析和出图十分方便),影响因素(或影响因素的类别)较多的时候,结果就会很难表达和解释。近年来,还有一些更为简洁明了的方差分解方法也被大家所认可接受,如以下两案例所示(Gross et al. 2017, Sirami et al. 2019),关于这一方法的详细原理和应用,下期博文将给大家带来详细讲解。
案例来源(Gross et al. 2017)
案例来源(Sirami et al. 2019)
参考文献:
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 & Evolution 1:0132.
Heikkinen, R. K., M. Luoto, M. Kuussaari, and J. Pöyry. 2005. New insights into butterfly–environment relationships using partitioning methods. Proceedings of the Royal Society B: Biological Sciences 272:2203-2210.
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.
Tang, Z., W. Xu, G. Zhou, Y. Bai, and Z. Xie. 2018. Patterns of plant carbon, nitrogen, and phosphorus concentration in relation to productivity in China's terrestrial ecosystems. Proceedings of the National Academy of Sciences 115:4033-4038.
####产生模拟数据
X1 <-c(2.75,2.5,2.5,2.5,2.5,2.5,2.5,2.25,2.25,2.25,2,2,2,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75)
X2 <- c(5.3,5.3,5.3,5.3,5.4,5.6,5.5,5.5,5.5,5.6,5.7,5.9,6,5.9,5.8,6.1,6.2,6.1,6.1,6.1,5.9,6.2,6.2,6.1)
X4 <- c(2017,2017,2017,2017,2017,2017,2017,2017,2017,2017,2017,2017,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016)
X3 <- c(12, 11,10,9,8,7,6,5,4,3,2,1,12,11,10,9,8,7,6,5,4,3,2,1)
Y <- c(1464,1394,1357,1293,1256,1254,1234,1195,1159,1167,1130,1075,1047,965,943,958,971,949,884,866,876,822,704,719)
###构建全模型, 此模型的R2即为整体的adj R2
lm0<-lm(Y ~ X1 + X2)
r0<-summary(lm0)
###只含X1的模型
lm1<-lm(Y ~ X1)
r1<-summary(lm1)
####只含X2的模型
lm2<-lm(Y ~ X2)
r2<-summary(lm2)
###则,图1其中C的求法:lm1的adj R2 加 lm2的adj R2 减去 lm0的adj R2
C<-r1$adj.r.squared+r2$adj.r.squared-r0$adj.r.squared
C
###同理,图1中B的求法
B<-r2$adj.r.squared-C
B
###与vegan包的结果做对比,完全一致
library(vegan)
res <- varpart(Y, X1, X2)
res
plot(res, bg=2:3, digits=3)
####如果你有三个自变量X1, X2, X3, 那么也一样
### 首先构造全模型
###构建全模型
tm0<-lm(Y ~ X1 + X2+ X3)
###那么如何求图2中的a
###构造一个只包含X2, X3的模型
tm1<-lm(Y ~ X2+ X3)
###然后,用tm0的adj R2 减去 tm1的adj R2, 就得到图2中的a
r_tm0<-summary(tm0)
r_tm1<-summary(tm1)
### a等于
r_tm0$adj.r.squared-r_tm1$adj.r.squared
### 与vegan结果作对比, 完全一致
library(vegan)
res <- varpart(Y, X1, X2, X3)
plot(res, bg=2:3, digits=3)
##其他几项的求法大家可以根据以上原理自行计算。
更多应用统计内容,欢迎大家关注我的微信公众号,二傻统计,就应用统计问题与大家交流
1/1 | 总计:4 | 首页 | 上一页 | 下一页 | 末页 | 跳转 |
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2025-2-22 04:34
Powered by ScienceNet.cn
Copyright © 2007-2025 中国科学报社