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

博文

广义线性混合模型(glmm)的R2计算之二项回归

已有 4569 次阅读 2022-7-4 22:46 |系统分类:科研笔记

首先必须明确,二项回归(binomial regression) 有两种形式,一种是针对因变量Y为0,1数据,也就是传说中的logistic 回归, 这里的自变量Y必须不是0,就是1,而不能是0-1之间的某个值,这其实是一种伯努利回归,抛硬币,记录每次的结果,就是这种模型的典型。

比如下面这个数据:这里的Colour这个变量,就是个0,1 变量,所以如果我们想要分析habitat和treatment对Colour的影响,就要用到binomial regression了,或者,更确切的说是logistic regression。

图片1.png

对于这个数据,我们可以用glmm模型来分析,

m7 <- glmer(Colour ~ Treatment + Habitat + (1 | Population) + (1 |Container),

family = binomial(link ="logit"), data = DataMale)

我们可以直接用partR2函数来求其R2

图片2.png

同时我们也可以用,MuMIn 下面的r.squaredGLMM函数来求其R2

图片3.png

partR2这里的固定效应R2默认计算方法和r.squaredGLMM下的delta计算方法相同。

我们也可以通过在partR2中指定expct="liability"求得与r.squaredGLMM相同的theoretical相同的结果:

partR2(m7,expct="liability")

图片4.png

所以,对于0,1数据或者logistic 模型来讲,partR2 与r.squaredGLMM都能给出固定效应模型正确的R2。这是最简单的binomial regression 模型,没什么争议,也没什么误区,在论文撰写中,你只需要说清楚你是用的哪个函数求得R2就可以了。

上述算法的原理也简单明了,比如theoretical R2的算法,其原理如下,

图片5.png

而其中的delta算法则稍显复杂,先把M7中的固定效应全部删掉,构造一个新的模型M8:

m8 <- glmer(Colour ~ 1 + (1 | Population) + (1 | Container),

                        family = binomial(link = "logit"),

                        data = DataMale)

Delta R2的算法与theoretical R2算法的根本区别同样在于残差的方差计算方法不同。其具体算法如下:

图片6.png

好了,下面我们看二项回归的另一种形式: 因变量Y为计数数据转换而来的比率数据。看下面这个数据:比如我们每次都测量一定数量的个体,然后记录其颜色,其中有些个体为绿色,有些为棕色,那么这时候,我们分析的对象Y其实是每次实验结果中不同颜色个体所占的比率,但这个比率每次都是由计数数据得来的,比如第一观测数据的比率是绿色个体所占比率为2/(25+2)。所以这里,每次分析的数值都是0,1之间的一个值,而不是0,1本身。这是该模型和上面logistic模型的本质区别。

图片7.png

对于这个数据,比如我们想分析不同的生物因子对颜色比率的影响,那么可以:

m9 <-glmer(cbind(nGreen, nBrown) ~ Bio7 + Bio14 + Bio17 + Bio19 +

              (1|SiteID) + (1|OLRE), data=Grasshoppers, family="binomial")

注意,这里还加入了一个观测水平的随机效应以处理可能出现的过度离散。对于这个模型的R2, 我们同样可以用partR2与r.squaredGLMM来求。结果如下图所示

图片8.png 


那么,问题来了,这两个函数给出来的结果差别非常大。这是怎么回事呢?

真相是,partR2的结果竟然是错的partR2虽有日本大师亲自上阵开发(Stoffel et al. 2021),但对于比率数据(仅指由计数数据转化而来的比率数据),partR2给出的结果确实是错误的。错在哪里呢,其实也是低级失误,就是没有正确的计算这种比率数据的样本量,而是错误的认为这种比率数据的样本量和0,1数据一样,每个结果都是基于一个样本得到的(比如一次掷硬币得到的结果)。但事实上,这种比率数据的每一个比率并不是由一个样本得到的,比如这个案例数据里的第一个观测样本,就是由2+25=27个样本得到的,而不是像投硬币一样是一次得到的。

事实上,统计上这种数据的方差有明确的定义即var=p(1-p)*n, 这里的n就是样本量,对于0,1数据或伯努利数据,显然n=1。但对于比率数据来讲n是多少呢?其实n就是所有比率的平均样本量。比如第一个比率值是基于10个样本得到的,第二个比率值时基于20个样本得到的,那么平均样本量就是(10+20)/2=15。真相,就是这么简单。

比如针对上面的案例,partR2给出来的错误结果其实是这么得来的:

图片9.png

正确的算法该如何呢?其实只需要把正确的样本量加入到计算过程中就可以了,如以下代码所示:

图片10.png

所以,对于用glmm模型分析由计数数据转化而来的比率数据,记住一点,partR2给出来的R2结果是错误的,而r.squaredGLMM给出的R2结果则是正确的。

对于delta R2的计算原理,也不复杂,代码如下。相信看完了上面和鄙人其他连续的几篇博文后,读者自己也能理解。

图片11.png

特别说明,本博文中案例数据和原始代码来自于(Nakagawa et al. 2017),处于演示和对比目的,鄙人对原始代码有所改动,再次也向日本大师Nakagawa表示致谢,对其在应用统计领域做出的杰出贡献致敬!

参考文献:

Nakagawa, S., P. C. D. Johnson, and H. Schielzeth. 2017. The coefficient of determination R(2) and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. Journal of the Royal Society Interface 14:20170213.

Stoffel, M. A., S. Nakagawa, and H. J. P. Schielzeth. 2021. partR2: Partitioning R2 in generalized linear mixed models.  9:e11414.

 




https://m.sciencenet.cn/blog-3442043-1345865.html

上一篇:复杂回归模型在生态学中的应用 PPT分享
下一篇:Meta 分析中多个处理共用一个对照的解决方案

0

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

数据加载中...

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

GMT+8, 2024-3-29 05:28

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部