真香预警!想发临床预测模型高分SCI?还得看这个方法!
混合效应模型第二讲_线性混合效应模型
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:表示了我们混合效应模型的关系式,我们把海滩作为随机截距添加到模型的关系中,简单理解,随机截距其实就相当于我们以前线性回归中未曾考虑过的混杂因素
标记2:模型的一些评价指标,我们重点关注AIC即可
标记3:代表了随机效应的系数,我们可以很清楚的看到,Beach作为随机截距对响应变量造成的方差为7.507,然后模型不可以解释的信息(残差)对响应变量造成的方差为9.11
标记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'splotthe
ggplot(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) foreachbeach
geom_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 得到的估计值之间的差异也随之增大(存在分歧)
好,那么到这里,线性混合效应模型就给大家介绍到了这里,很基础的内容,希望可以帮助大家踏入混合效应模型的大门,然后大家有什么不明白的也可以在评论区提问,欢迎各位小伙伴积极讨论~
我是晨曦,我们下期再见
参考教程:
- Zuur, A. et al. 2009. Mixed effects models and extensions in ecology with R. Springer
- Bates, D. et al. 2015. Fitting linear mixed-effects models using lme4. Journal of Statistical Software 67: 1-48.
- Fitzpatrick, C. R., Mustafa, Z., and Viliunas, J. Soil microbes alter plant fitness under competition and drought. Journal of Evolutionary Biology 32: 438-450.
- Crossed vs. Nested random effects
- Fixed vs. Random effects
- ML vs. REML
- Linear mixed-effects models (uoftcoders.github.io)
- methods-supported-by-effects.pdf (r-project.org)
晨曦的混合效应模型系列传送门
晨曦的空间转录组笔记系列传送门
晨曦单细胞文献阅读系列传送门
1. 非肿瘤单细胞分析模板已到位!眼馋单细胞的小伙伴快来看!手把手教你产出第一篇单细胞SCI!
晨曦单细胞笔记系列传送门
晨曦单细胞数据库系列传送门
—
END—
撰文丨晨 曦
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
关键词
就是
数据集
生信
结果
模型
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。