混合效应模型第二讲_线性混合效应模型

Hi,大家好,我是晨曦
今天这期继续我们的混合效应模型专题,我们这回将从基础开始介绍混合效应模型,相信跟完本期教程的小伙伴至少可以对混合效应模型有一个大概的印象
引言
既然谈到了线性混合效应模型,那么我们自然就要和单纯的线性回归进行一下前后的对比,让我们花费一点时间来快速回归一下线性回归的数据假设:
1. 残差符合正态性
2. 残差符合方差齐性
3. 样本之间彼此独立
4. 预测变量和响应变量之间存在线性关系
但是我们的生活中的数据大部分都是嵌套数据,也就是说样本与样本之间彼此并不独立,比如学校和班级、医院和科室、国家以及多中心医院等等
对于没有掌握混合效应模型的我们只能忽略这种嵌套形式,强制性认为我们的数据是彼此独立的,但是既然我们知道了混合效应模型,我们自然就不能这么简单的认为了
晨曦解读
简单来说混合效应模型其实就是考虑了对结果有影响的随机效应(混杂因素),而普通的线性回归压根就没有考虑这点
当然理论的讲解可能略显空洞,接下来我们将通过代码来讲解整个混合效应模型的构建过程
输入数据集介绍
The RIKZ dataset——不同海滩物种以及某些指标的结果
晨曦解读
数据集的具体内容不需要我们了解,我们只需要知道数据集是这种嵌套形式的,这个数据不就是很像我们多中心临床研究吗?
不同地区中的不同位置测量相同的指标
不同医院中的不同患者测量相同的指标(多中心研究某种治疗因素的影响)
准备工作
library(tidyverse)#数据操作library(lme4)#混合效应模型library(lmerTest)#模型检验
准备输入数据
rikz_data <- "https://uoftcoders.github.io/rcourse/data/rikz_data.txt"download.file(rikz_data, "rikz_data.txt")rikz_data <- read_delim("rikz_data.txt", col_names = TRUE,                       delim = "\t")rikz_data <- rikz_data %>% mutate(Beach = as.factor(Beach))#分类变量转换为因子的形式
晨曦解读
我们可以很清楚的看到,这个数据就是一个典型的嵌套数据,不同海滩分别取了不同的五个位置然后测量一些相关数据
拟合简单的线性回归方程
既然我们知道要拟合混合效应模型,那么好奇同样的数据在普通的线性回归方程中的表现自然是一种很常见的心理
#拟合线性模型basic.lm <-lm(Richness~ NAP, data = rikz_data)summary(basic.lm)#可视化ggplot(rikz_data, aes(x = NAP,y = Richness)) +geom_point() +geom_smooth(method = "lm") +theme_classic()
晨曦解读
当我们相信这个线性回归的结果的时候,我们需要对线性回归模型进行模型诊断,好让我们对这个模型的表现或者说是模型的假设有一个详细的了解
#模型诊断par(mfrow=c(2,2))plot(basic.lm)
晨曦解读
左上角的可视化表示代表残差和拟合值的关系,正常来说应该是没有关系的,我们可以很清楚的看到这个可视化结果是比较糟糕的,因为随着拟合值的增加明显和残差有着一定的关系
右上角的可视化表示残差是否符合正态,正常来说应该是点的分布沿着对角线或者上下偏移幅度不大,我们可以很清楚的看到这个可视化结果明显是不符合正态的
左下角的可视化则是对数据进行了转换,表达的内容和左上角的是类似的
右下角的可视化则是展示了离群值的,当然这个图其实我们也很少看
其实看到这里我们就已经知道,我们的数据肯定不可以去使用普通的线性回归模型,而且最重要的一点则是,数据本身是嵌套数据,违反了数据独立性原则,所以我们需要使用混合效应模型
拟合混合效应模型
#随机截距mixed_model_IntOnly <- lmer(Richness ~ NAP + (1|Beach), data = rikz_data, REML = FALSE)summary(mixed_model_IntOnly)
晨曦解读
标记1:表示了我们混合效应模型的关系式,我们把海滩作为随机截距添加到模型的关系中,简单理解,随机截距其实就相当于我们以前线性回归中未曾考虑过的混杂因素
随机效应:水平数最好大于5,然后具有可以互相替换的性质,比如说位置1、位置2、位置3...,本身可以随意替换,但是下雨、刮风就不可以
标记2:模型的一些评价指标,我们重点关注AIC即可
标记3:代表了随机效应的系数,我们可以很清楚的看到,Beach作为随机截距对响应变量造成的方差为7.507,然后模型不可以解释的信息(残差)对响应变量造成的方差为9.11
我们可以查看一下上面我们没有纳入随机效应的模型(线性回归)发现残差的标准差为4.16,则方差为17.3056,大致等于混合效应模型中的随机效应+残差,所以我们可以大致得到一个简单的结论,在线性回归中考虑了混杂因素(随机效应)就是混合效应模型
标记4:固定效应的系数(结果类似于线性回归)
标记5:预测变量之间的相关性,如果等于±1则存在多重共线性,要考虑删掉一个预测变量
再次划重点!!!
1. Note that it’s considered best practice that random effects have at least 5 levels, otherwise it should be used as a fixed effect(作为随机效应至少包含5个水平甚至以上,否则应该把其设置为固定效应)
2. 关于随机效应的方差可以简单理解为:与海滩效应相关的方差是7.507。换句话说,考虑到模型中的固定效应,海滩之间的差异占残差方差的45% (7.507/7.507 + 9.111) * 100 = 45% 。注意,这里的分母是总方差(即所有随机效应(包括残差)的方差之和)
#可视化rikz_data$fit_InterceptOnly <-predict(mixed_model_IntOnly)# Let'splottheggplot(rikz_data, aes(x = NAP,y = Richness,colour = Beach)) + # Addfixedeffectregressionline (i.e.NAP)geom_abline(aes(intercept = `(Intercept)`, slope = NAP),size = 2,as.data.frame(t(fixef(mixed_model_IntOnly)))) + # Addfittedvalues (i.e.regression) foreachbeachgeom_line(aes(y = fit_InterceptOnly),size = 1) +geom_point(size = 3) +theme_classic() +theme(legend.position = "none") +scale_colour_brewer(palette="Set1")
晨曦解读
可以很清楚的看到,我们每条之间斜率是相同的,但是截距不同,所以这个就是随机截距模型(混合效应模型)
粗黑线对应于与模型的固定效应分量(即6.58 -2.58(x))相关的拟合值。细的彩色线条对应于每个海滩的估计拟合值。如您所见,正如预期的那样,它们都有单独的截距。随着随机效应的估计方差的增加,这些线会变得更加分布在粗黑线周围。如果方差为0,则所有的有色线条将与粗黑线条重合
前面我们讨论的是基于不同的Beach,NAP和Richness之间的关系,但是我们这回通过对数据以及生物学背景知道,针对不同的Beach,NAP和Richness之间的关系是不同的,前面我们假设的是不同Beach,但是NAP和Richness之间的关系是相同的,比如说某些Beach上NAP和Richness之间的关系就很差,有些Beach则是很好,那么这种斜率不同的自然就需要另外一种混合效应模型
#随机斜率|随机截距模型mixed_model_IntSlope<- lmer(Richness ~ NAP + (1 + NAP|Beach), data = rikz_data, REML = FALSE)#当然我们换一种写法可能会更加直观一点(结果是一样的)mixed_model_IntSlope<- lmer(Richness ~ NAP + (NAP|Beach), data = rikz_data, REML = FALSE)
晨曦解读
代码运行到这里会出现一个警告:boundary (singular) fit: see help('isSingular')
这个是警告我们模型未能完全拟合(畸形拟合)或者说是因为我们的随机效应过小而导致的,网络上对此也是没有一个统一的看法,但是大部分人并没有关注这里,所以我们可以选择略过
summary(mixed_model_IntSlope)
晨曦解读
结果大致和前面的类似,这里唯一的区别是随机效应中的附加斜率的方差,它估计了海滩上NAP的方差。它还包括一个 Cor 项,用于估计截距和斜率方差之间的相关性。在 -1,这意味着如果某一个Beach造成了响应变量的波动较大,那么此时其也拥有更大负斜率的NAP
#可视化rikz_data$fit_IntSlope <-predict(mixed_model_IntSlope)ggplot(rikz_data, aes(x = NAP,y = Richness,colour = Beach)) +geom_abline(aes(intercept = `(Intercept)`, slope = NAP),size = 2,as.data.frame(t(fixef(mixed_model_IntSlope)))) +geom_line(aes(y = fit_IntSlope),size = 1) +geom_point(size = 3) +theme_classic() +theme(legend.position = "none") +scale_colour_brewer(palette="Set1")
晨曦解读
可以很清楚的看到固定效应还有截距和斜率的关系变化
这里我们需要注意的是模型的构建并不是一个固定流程,比如说我们也可以没有固定效应而直接构建模型
#没有固定效应mixed_model_NoFix <- lmer(Richness ~ 1 + (1|Beach), data = rikz_data, REML = TRUE)summary(mixed_model_NoFix)
那么到这里混合效应模型或者说是线性混合效应模型其实就是这些内容,但是我们仍然有点疑惑所以我简单整理了一些问题:
提问:什么是随机效应?什么是固定效应?
回答:固定效应简单理解可以看作是我们感兴趣的预测变量,我们感兴趣这些预测变量和响应变量之间的关系,而随机效应则是我们并不感兴趣但是确实存在并且影响着响应变量(最好随机效应内的水平大于5个)

提问:线性模型估计参数的方法是OLS(最小二乘法),那么线性混合效应模型估计参数的方法是什么?
回答:线性混合效应模型估计参数的方法是最大似然法(ML)和限制性最大似然法(REML),简单来说,如果你对准确估计随机效应感兴趣,你应该用 REML 来拟合模型,而如果你对估计固定效应感兴趣,你应该用 ML 来拟合模型,鉴于我们通常对固定效应最为感兴趣,所以ML可能是我们用的比较多的,随着模型中参数数量的增加,ML 和 REML 得到的估计值之间的差异也随之增大(存在分歧)
好,那么到这里,线性混合效应模型就给大家介绍到了这里,很基础的内容,希望可以帮助大家踏入混合效应模型的大门,然后大家有什么不明白的也可以在评论区提问,欢迎各位小伙伴积极讨论~
我是晨曦,我们下期再见
参考教程:
  1. Zuur, A. et al. 2009. Mixed effects models and extensions in ecology with R. Springer
  2. Bates, D. et al. 2015. Fitting linear mixed-effects models using lme4. Journal of Statistical Software 67: 1-48.
  3. Fitzpatrick, C. R., Mustafa, Z., and Viliunas, J. Soil microbes alter plant fitness under competition and drought. Journal of Evolutionary Biology 32: 438-450.
  4. Crossed vs. Nested random effects
  5. Fixed vs. Random effects
  6. ML vs. REML
  7. Linear mixed-effects models (uoftcoders.github.io)
  8. methods-supported-by-effects.pdf (r-project.org)
晨曦的混合效应模型系列传送门
晨曦的空间转录组笔记系列传送门
晨曦碎碎念系列传送门(未完待续...)
1. 想白嫖单细胞生信文章?这五大源头数据库,是你发文章的源泉!高频预警!你一定要收藏!
2. 盘活国自然的新思路!你研究的热点真的是热点吗?大数据帮你定位!
3. 好家伙!90%以上审稿人都会问到的问题,今天帮你解决!就是这么齐齐整整!
4. 没想到!生信分组还有这个大坑!你被坑过吗?!
5. 关于富集分析这件事,我有话想说。。。
6. 好御好高级!CNS级别美图是如何炼成的?看这篇就懂了!
7. 化繁为简!一文帮你彻底搞懂机器学习!想发高分文章,这篇是基础!
8. 你不知道的机器学习算法!关键时候能救命!
9. 致命!芯片&测序的联合到底能不能联合分析?审稿人最爱用这刁难你!
10. 躲不过的树!80%的生信SCI中都见过它!你真的搞懂了吗?
11. Python or R? 哪个更适用于生信发文章?深入浅出给你讲透!
12. 生信和抖音是一样的算法原理?不仅让你成瘾,也能发高分文章!
13. 跟3-5分SCI相比,CNS里的生信玩的可太花了!其实简单的离谱!
14. 揭秘!小鼠和人的免疫浸润分析有何区别?看这篇就够了!
15. 临床预测模型中的宠儿!最常见的机器学习 算法,没有之一!直接拿来用 !
16. 临床预测模型评价,不只有ROC,这个指标你遗漏了吗?
17. 肺肿瘤机器学习模板奉上!还不赶快产出2022年的你的第一篇SCI?!
晨曦单细胞文献阅读系列传送门

1. 非肿瘤单细胞分析模板已到位!眼馋单细胞的小伙伴快来看!手把手教你产出第一篇单细胞SCI!

晨曦单细胞笔记系列传送门
晨曦从零开始学画图系列传送门
1. 看完这篇,彻底掌握生信画图精髓!超级实用,我不许你不知道!
2. 想让SCI看上去更高逼格?这些绘图技巧你一定要知道!
3. 3min掌握SCI配色神技,学会你就是组会汇报上最靓的仔!
晨曦单细胞数据库系列传送门

END

撰文丨晨   曦
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文