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

博文

广义线性混合模型(GLMM)伪泊松回归的R2计算

已有 4890 次阅读 2022-3-29 22:51 |系统分类:科研笔记

前面的博文给大家介绍了一般线性混合模型和广义线性混合模型中泊松回归时R2的计算方法,今天跟大家聊下GLMM模型中的伪泊松回归中的R2计算问题。首先我们必须清楚两点,

1)伪泊松是用来分析自变量对计数数据的影响的,所以因变量Y必须是计数数据。

2)对于计数数据,理论上,泊松模型假设其均值=方差。但事实上,现实数据中,计数数据的方差往往大于均值(即出现了所谓过度离散),这时候,伪泊松(以及负二项)回归就要出场了。

好了明白了什么时候用泊松回归之后,我们看一个案例,这个案例也来自于日本大师的文章(Nakagawa et al. 2017)。

比如我们要分析Treatment和Habitat两个因素对,egg数量的影响,考虑其数据特点和结构之后后决定采用GLMM-伪泊松回归模型,我们称这个模型为full model,记为fecmodPQLf。

library(MASS)

fecmodPQLf <- glmmPQL(Egg ~ Treatment + Habitat, random = list(~1 | Population,  ~1 | Container),  family = "quasipoisson", data = DataFemale)

那么这个模型的R2该如何求呢?按下列步骤就可以了:

1) 重新构造一个null model, 这个模型不含任何固定效应,记为fecmodPQLr,即:

图片1.png

2)依据null model,求eggs理论上的期望值的均值lambda

构造null model的目的,就是为了求eggs理论上的期望值的均值lambda。因为伪泊松回归的link函数也是log link, 所以我们这里要用指数函数exp将其转换到原始尺度上来:

图片2.png

其中,fixed函数提取模型截距的参数估计值,VarCorr(fecmodPQLr)[2,1]提取随机效应population的方差,VarCorr(fecmodPQLr)[4, 1]提取随机效应container的方差,

其实我们也可以用VarCorr函数直接查看这两部分方差,即箭头所指数字。如果你不理解为什么要采用上述公式,那么请参看鄙人上一篇博文,总之这是日本大师改进后的算法。

图片3.png

3)求full model的过度离散系数omegaF

得到了lambda之后,我们还需要求得另外一个参数:omegaF,即模型的过度离散系数

omegaF <- as.numeric(VarCorr(fecmodPQLf)[5, 1])

事实上,omegaF 就是VarCorr函数显示的full model(重点,这里是full model,而不是null model里)残差的方差(这里只是这么显示,这并不是模型真正的残差的方差),即下图红色箭头所指数字。

图片4.png

 

得到了lambda和omegaF之后,我们就可以由三种方法得到full model真正意义上的的残差的方差,即delta法,lognormal法和trigamma法,其计算方法如下:

图片5.png

得到了full model残差的方差之后,full modelR2的计算就相当简单了。其中:

marginal R2=固定效应的方差 / (固定效应的方差+随机效应的方差+残差的方差)

conditional R2=(固定效应的方差+随机效应的方差) / (固定效应的方差+随机效应的方差+残差的方差)。

本案例中,固定效应的方差为:

图片6.png

所以三种R2的值分别为:

图片7.png

 

图片8.png

 

图片9.png

其中,delta R2 是日本大师推荐的GLMM R2。但是,如果我们用MuMIn包(Version: 1.45.1及更老版本)中的r.squaredGLMM函数来求fecmodPQLf模型的R2,得到的结果跟我们上述计算结果明显不同

图片10.png

这是为什么呢?原来r.squaredGLMM在计算omegaF的时候出现了错误。经过鄙人查证,
原来这个函数错误地把omegaF的平方的倒数,即omegaF的-2次方,当做omegaF来进行运算,
所以出现了错误。比如我们调用r.squaredGLMM的核心函数r2glmm(这个函数要自己找到源代码,
并在R中运行以下才可以生效,源代码地址: 
如果我们正确的输入omegaF值,那么:


图片11.png

怎们样?跟我们手动计算的完全一致吧!这里的sigma2.glmmPQL函数可以帮我们直接提取模型的omegaF值,也来自于r.squaredGLMM函数的源代码中。

但是如果我们错误的把omegaF的-2次方当做omegaF输入到r2glmm函数中,那么就得到跟r.squaredGLMM函数得到的一模一样的结果,但很遗憾,这是个由于低级失误造成的错误结果。

图片12.png

图片14.png

鄙人于春节假期期间发现这个错误后,就发邮件告知了MuMIn这个包的作者Kamil Bartoń大师(Institute of Nature Conservation PAS, 这哥们的单位PAS是啥我也没看懂)。Kamil Bartoń大师已经修正了这些错误,并于2022-02-23日对MuMIn进行了升级,即目前最新的1.46.1版,这个版本的r.squaredGLMM的函数,终于产出了完全没有问题的结果,各位同仁可以大胆使用了!

图片15.png

Kamil Bartoń大师在其更新日志中,也对鄙人进行了致谢,如下:
MuMIn package news:
Version: 1.46.0 (2022-02-23)
(fixed) r.squaredGLMM: now gives correct values for negative binomial models (thanks to Shuang Zhang for research)


提及这个致谢,除了不能免俗的自我吹捧之外,请注意我给加粗的部分,这里Kamil Bartoń大师为什么说的是对negative binomial models的修正,而不是对伪泊松模型呢?这是因为原来的MuMIn包GLMM负二项模型R2的计算也存在同样的问题,且负二项和伪泊松有相通之处哈哈,这个问题,我们下次再聊!

最后需要指出的是,MuMIn是生态学领域引用最多的R包之一,根据植物所赖江山老师研究,该包的总引用次数在生态学领域所有的R包中排第五(Lai et al. 2019)。在这些引用中,又有多少论文是得到了错误的R2并且成功发表呢?所以没有什么比搞清楚统计的基本原理最重要了。总之,人在江湖,还是要多加强学习,多多交流啊!

最后欢迎大家关注鄙人的个人微信公众号“二傻统计”(可通过鄙人头像二维码扫码关注),大家互相学习共同进步!

 

参考文献:

Lai, J., C. J. Lortie, R. A. Muenchen, J. Yang, and K. Ma. 2019. Evaluating the popularity of R in ecology. Ecosphere 10:e02567.

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.




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

上一篇:广义线性混合模型之泊松回归模型的R2计算
下一篇:复杂回归模型在生态学中的应用 PPT分享

1 郭新磊

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

数据加载中...

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

GMT+8, 2024-3-29 21:13

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部