我们在实际数据分析中，经常会遇到这样的问题，比如，我想同时考虑多个因子，如X1, X2, X3等等，如何对比这些因子对因变量Y的影响的大小呢？事实上，这个问题并没有想象中的简单，处理这个问题，可以先弄清以下问题。

1）首先必须明确，多因素对Y的共同影响，并不一定等于单因素对Y的影响的作用之和。事实上，只有在不同的X彼此独立，即X之间完全没有相关性的前提下，多个X对Y的总影响才等于单个X对Y的影响之和。

比如，现在有一个模型

M0 <- lm(y~x1+x2+x3, data=d1)；

假定其R²=30%, 那么30%是否就等于以下三个模型的R²之和呢？

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的总体变异，也就是总方差V_T, 而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)的R²

A+C,    m2<-lm(Y~X1) 的R²

B+C,    m3<- lm(Y~X2) 的R²

那么A, B, C 对应的计算方法为：

首先求C, C=(A+C)+(B+C)-(A+B+C)

然后求B,  B=(B+C)-C

最后求A,  A=(A+C)-C

当然，以上各项也可以直接通过另一种方法求得，比如我们要求A,

那么首先，求方程：

m3<-lm(Y~X2) 的R²，

然后用m1的R² 减去m3的R² ，就是A

同样原理假设你有三个因素，如，环境，资源和微气候，对Y都有影响，那么和区别其各自的作用呢？

如图2案例(Heikkinen et al. 2005)：

 

首先三个因素对Y的总的解释量是R²=1-41.4%=58.6%, 也就是模型

M1<-lm(y~habitat + resources +microclimate)的R²=58.6%,

那么接下来，比如我们要求a，也就是habitat单独对Y的影响，从原理上，这个影响就是我们同时控制了resources 和microclimate之后，habitat对Y造成的影响。所以我们首先构造一个模型

M2<-lm(y~ resources +microclimate)

然后求出该模型的R²

那么用M1的R²减去M2的R²就是a,

同样道理，我们可以依次求出，b, c.

a, b, c求出之后，我们便可求，d, e, f, g的值，

比如，求d, 首先求模型

M4<-lm(y~ microclimate) 的R²

然后，以模型M1的R²减去M4的R²，其值即a+b+d,

因为a, b已知，则：

d=(a+b+d)-(a+b),

同理可以求出e, f,

最后，求g, 因为我们已知模型

M6<-lm(y ~ habitat) 的R² =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&ouml;yry. 2005. New insights into butterfly&#x2013;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&ntilde;o, J. Bosch, J. A. Barrientos,    A. Ricarte, M. &Aacute;. Marcos-Garcia, J. Mi&ntilde;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.

附案例及计算代码（与vegan包结果对比）：

####产生模拟数据

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

A<-r1$adj.r.squared-C
A

###同理，图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)

##其他几项的求法大家可以根据以上原理自行计算。

更多应用统计内容，欢迎大家关注我的微信公众号，二傻统计，就应用统计问题与大家交流

转载本文请联系原作者获取授权，同时请注明本文来自张霜科学网博客。
链接地址：https://blog.sciencenet.cn/blog-3442043-1242874.html


下一篇：关于Meta分析，你需要了解的十个技术问题