伟大的革命导师鲁迅他老人家早就教导我们“尽信R, 则不如无R”！软件本身是傻的，你输入命令，他就输出结果，至于模型是否合适，结果是否可信，甚至是对还是错，大家可能并不真正了解。

最近跟植物所赖江山老师测试他新开发的glmm.hp包（这个包主要用于混合模型中多元回归时R2在不同因子间的分解），自己也再次学习了各种混合模型R2计算原理。真是不究不知道，一究吓一跳。目前，混合模型(主要是广义线性混合模型)的R2有多种算法，不同的R包函数给出的结果，有时有明显的差别，甚至一些函数的计算的结果都是错的（R包由不同用户各自开发，没有专业人士统一严格审核，函数出错也是常有的事儿）。鄙人试图以几篇博文力争给大家介绍清楚混合模型目前的主流算法以及一些相关的包存在的局限乃至错误，也算是对自己学习心得的一个小结。

首先要明确，不同的混合模型具有不同的R2算法。比如线性混合模型和广义线性混合模型的算法不同。在广义线性混合模型中，针对不同数据分布类型的模型，比如二项、泊松、gamma, 伪泊松，负二项等等其算法也各不相同，所以这就大大增加了混合模型R2计算的复杂性和迷惑性。我们接下来从最简单的情况开始，先看下线性混合模型中R2的算法。

如果你用的是线性混合模型，也就是linear mixed effect model（lmm）, 那么恭喜你，你的问题就很简单，很容易处理了。你只需要知道，lmm用于因变量（以下都以y代替）的残差理论上符合正态分布数据就可以了。目前lmm模型R2最常用的是日本人定义的算法(Nakagawa and Schielzeth 2013), 把lmm模型的R2分成两部分， 一个叫marginal R2, 也就是固定效应造成的R2, 一个叫conditional R2, 也就是固定效应+随机效应共同造成的R2。

这里多说几句，日本人Nakagawa是应用统计领域真正的当红大师，他发表的一系列论文引用率都高的吓人。比如2013年这篇关于混合模型R2计算的论文，目前引用率已经达到了7369次。但也可能是狂热粉丝太多，这位大师似乎不大愿意回复咨询问题的邮件。毕竟大师，都有自己的性格啊！闲话少说，言归正传，我们以下面这个模型为例：

m1<- lmer(Biomass~Year+Temperature+Precipitation+SpeciesDiversity+

(1|Population), data = biomass)

  他的R2该如何计算呢，这个并不复杂！首先我们提取该模型固定效应的参数估计值和自变量矩阵，然后二者相乘，依据固定效应，对应的Y的期望值fm1：

fm1<- model.matrix(m1)%*%fixef(m1)

那么fm1的方差(var(fm1))，就是固定效应（模型中四个因素一起）总体上对Y造成的方差。除了固定部分的方差，混合模型还包括随机效应方差和残差的方差，其实这两项方差的值，模型已直接给出，我们用summary()函数查看模型结果，得到下图，其中上下红色箭头所指部分，即是lmm模型中随机效应造成的方差和残差的方差。

我们可以用VarCorr函数可以直接提取随机效应的方差和残差的方差。

d1<-as.data.frame(VarCorr(m1))

Y的总方差=固定效应方差+随机效应方差+残差方差。用固定效应造成的方差除以Y的总方差，就得了日本大师定义的所谓marginal R2， 即：

用固定效应方差与随机效应方差之和除以总方差，就得到了conditional R2

当然，我们可以用MuMIn包下的r.squaredGLMM函数直接得到lmm模型的R2

     也可用partR2包下的partR2命令直接得到marginal R2，
 

这些结果与手动计算的结果完全一样。partR2函数其实是用来进行方差分解的，我们可以用它来衡量模型中各个因子各自独立的R2，比如我们运行以下代码，就会得到：

由这一结果可知，物种多样性对R2的贡献远大于年份、温度、降雨这些指标，这说明生物量主要由物种多样性决定。那么这个R2（确切的讲，叫part R2, 注意，不是partial R2）是如何得来的呢？

其实也不复杂！比如，我们想求speciesdiversity的part R2, 那么我们在模型m1（full mdoel）的基础上，再构造一个模型m2(reduced model), 这个模型其实就是在m1中把speciesdiversity删掉，但其他因素和模型结构都不变。即：

m2<- lmer(Biomass ~Year+Temperature+Precipitation+(1|Population),data = biomass)

然后我们求模型m2中，依据固定效应参数估计值得到的Y的期望值fm2：

fm2<-model.matrix(m2)%*% fixef(m2)

其方差为var(fm2), 现在，我们以full model中固定效应的方差var(fm1)减去reduced model中固定效应的方差var(fm2)，然后除以full model的总方差，即speciesdiversity的part R2, 也就是运行以下代码：

结果与partR2中结果完全相同。这里需要说明的一点，理论上，full model的总方差=reduce model 的总方差= y的原始方差var(y)。但由于混合模型中，随机效应的方差是由于参数估计得来的，而不是由var(y)-固定效应方差-残差方差得来的，所以实际应用中，full model的总方差，reduce model 的总方差，及y的原始方差 var(y)会有差异，但一般较小。比如本案例中三者差异如下

    如果你仔细看partR2方差分解的结果，就会发现，四个因素各自的part R2之和，并不等于模型的marginal R2。这是因为，四个自变量对Y的影响并不是彼此完全独立的，他们对Y有共同影响的部分，没有算在各自的partR2里。 这也正是赖江山老师最近开发的glmm.hp包所要解决的问题，关于这个包的原理与应用，赖老师后续也会有专门的文章介绍。

    到这里，我们会觉得对于lmm模型，R2的计算与分解似乎都很简单，有现成的命令可以一键搞定。但是如果我们进入广义线性混合模型领域，那么R2的计算就没这么简单了。现成的R2包无论是r.squaredGLMM或是partR2都有局限甚至是错误，我们也很有必要厘清这些问题。那么，在后续的博文中，我们将陆续给大家介绍二项、泊松、伪泊松、负二项、gamma等glmm模型的R2计算方法与问题，希望能增进大家对混合模型的理解，敬请大家持续关注！

     最后敬请大家关注本人的微信公众号“二傻统计”（本人头像二维码），与大家互相交流，共同提高!

参考文献：

Nakagawa, S., and H. Schielzeth. 2013. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4:133-142.