ixxmu
commented
1 year ago 线性混合模型 by 芬菲随笔
本文内容:
1、模型原理
2、模型公式
3、结果解读
Part1 模型原理
实验设计导致数据嵌套方差解释的逐渐补充图
1-1 拟合自变量和因变量的关系
数据受到背景噪音的影响，数据间存在聚集 真实数据来源的复杂性
crossed设计。每个样本同时属于两个或多个或所有上一层级的类别，在实验中，图a的下面两层就类似两个被试都接受了三种实验条件。
nested设计。数据中每个最底层的样本是独立的，上一级是可以将下一层级的样本进行分类组合的模式；
当然nested和crossed设计也可以如c,d组合来组合去搞人脑子，不过为了头发，咱们进行实验设计的时候就搞搞好吧。
线性混合模型中三种方差和两种效应的关系图
在一般线性模型中，其自变量全部为固定效应自变量，而线性混合模型中，除了固定效应自变量外，还包含了随机效应自变量。
使用一般线性模型，需满足3点假设：
1、正态性，因变量y服从正态分布；
2、独立性，不同类别y的观察值之间相互独立，相关系数为零；
3、方差齐性，不同类别y的方差相等。
然而，在线性模型中，当存在数据分层时，可以分开处理分层数据，并且要求分层因素不能过多。因而，不如换一个模型，如线性混合模型（linear mixed effects model, LMM），该模型只要求因变量y服从正态分布。
1-2 固定效应和随机效应自变量分类
固定效应和随机效应变量的区别之处关键在自变量的类别。
如果抽样数据中，包含一个自变量的所有类别，则将该变量作为固定效应。比如性别，只要抽样的数据中同时包含了两种性别，就可以将性别作为固定效应自变量。
如果一个自变量在抽样数据中只包含部分类别，则将该变量作为随机效应自变量。简言之，抽样数据集中的自变量可以包含该自变量的所有情况，则作为固定效应，如果只能代表总体的一部分，则作为随机效应。
另外，如果一个效应在不同的参与者样本中会持续存在，那么该效应就被认为是固定的。
随机效应是从一些人群中抽样的离散单位，它们本身就是分类的。
固定效应对平均趋势进行建模，而随机效应则对这些趋势在某些分组因素（如参与者或项目）的不同水平上的变化程度进行建模。
总之，相比一般线性模型，混合线性模型的目的是通过引入随机效应，得到一个更好的解释数据的模型。
Part2 模型公式
2-1 混合效应模型的种类
根据截距和斜率是否随机，混合效应模型可以有四种情况：
①随机截距+随机斜率；②随机截距+固定斜率；③固定截距+随机斜率；④固定截距+固定斜率。
2-2 随机效应模型的选择方案
做加法
在线性模型中加入多少随机效应项而言，一些学者认为应该采用渐进递增的方式进行，即从每一个随机效应项的截距到斜率逐渐叠加到已有模型，并使用卡方检验比较AIC和BIC等参数是否有显著提升来判断模型是否为最优。做减法
然而，另一些学者提倡一次性将所有“符合逻辑”的随机效应项一次性纳入模型当中，再检验这些随机效应项对模型估计方差的影响。如果影响小，再剔除这些随机效应项。
2-3 模型的选择
REML criterion at convergence限制最大似然 (REML) 方法，指残差最大似然法，它是LME模型默认的参数估计方法。
常用的两个模型选择方法——赤池信息准则（Akaike Information Criterion，AIC）和贝叶斯信息准则（Bayesian Information Criterion，BIC）。
赤池信息量准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型，优先考虑的模型应是AIC值最小的那一个。
AIC和BIC的原理是不同的，都是选择值最小的模型。AIC是从预测角度，选择一个好的模型用来预测，BIC是从拟合角度，选择一个对现有数据拟合最好的模型，从贝叶斯因子的解释来讲，就是边际似然最大的那个模型。
anova()函数将模型对象作为自变量，并返回一个anova，测试更复杂的模型是否比更简单的模型更善于捕捉数据。如果得到的p值足够低（通常小于0.05），我们得出的结论是，更复杂的模型明显优于更简单的模型，因此有利于更复杂的模式。反之，如果p值不够低（通常大于0.05），我们应该选择更简单的模型。
2-4 混合线性模型表达式
fit = lmer(data = , formula = DV ~ Fixed_Factor + (Random_intercept + Random_Slope | Random_Factor))data - 要处理的数据集；
formula - 表达式；
DV - 观测值 （因变量）；Fixed_Factor - 固定因子（自变量）；
() - 完整的随机因子(包括截距和斜率)；Random_intercept - 随机截距，即认为不同群体的因变量的分布不同（可以理解成有些人出生就在终点，而你是在起点...）；
Random_Slope - 随机斜率，即认为不同群体受自变量的影响是不同的（可以理解成别人花两个小时能赚10000元，而你只能挣个被试费...）；Random_Factor - 随机因子。
括号内是随机效应，括号外是固定效应。g表示分组变量，x、y分别表示自变量和因变量
截距中，1表示随机截距，0表示固定截距，默认截距为1。
1 | g表示固定斜率；x | g表示随机斜率
线性混合效应模型可以有四种情况：
随机截距 + 随机斜率；y ~ x + (x | g)
随机截距 + 固定斜率；y ~ x + (1 | g)
固定截距 + 随机斜率；y ~ x + (0 + x | g)
固定截距 + 固定斜率。y ~ x
Part3 结果解读
建立模型步骤如下：
1、先建立全模型，查看模型结果是否过拟合；
2、做减法，生成第二个模型；
3、模型比较，确定最终模型。
模型过度拟合问题：随机效应中变量间相关性接近1或-1，一般在0.9以上就认为模型过拟合，可以将过度拟合的随机斜率去掉。
本次以单个随机因素为例，对于多个随机因素，可能存在随机因素间的交叉和嵌套，本文暂不讨论。
本次示例数据为R语言内置数据集ChickWeight。
ChickWeight是578行4列的小鸡体重数据集，记录了不同喂养和育龄下小鸡的体重数据，采集了50只小鸡12个时间点四种喂养方式下体重数据。
weight：小鸡体重，单位是公克gm；
Time：育龄，小鸡出生天数；
Chick：小鸡序号，按最后体重从最轻到最重排序；
Diet：喂养方式，分为4种方式。
# Load packages
library(tidyverse)
library(ggplot2)
library(ggsci)
library(ggpubr)
library(cowplot)
library(lmerTest)
library(ggpmisc)
# 理解数据集
data(ChickWeight)
dim(ChickWeight)
## [1] 578 4
head(ChickWeight)
## Grouped Data: weight ~ Time | Chick
## weight Time Chick Diet
## 1 42 0 1 1
## 2 51 2 1 1
## 3 59 4 1 1
## 4 64 6 1 1
## 5 76 8 1 1
## 6 93 10 1 1
sapply(ChickWeight, class)
## $weight
## [1] "numeric"
##
## $Time
## [1] "numeric"
##
## $Chick
## [1] "ordered" "factor"
##
## $Diet
## [1] "factor"
# 两种喂养方式各选3只共6只小鸡
smallSet <- ChickWeight %>% subset(Chick %in% c(1:3, 21:23))
dim(smallSet)
## [1] 72 4
探究小鸡体重与育龄的关系
# 探究小鸡体重与育龄的关系
my_formula <- y ~ x
p1 <- smallSet %>%
ggplot(aes(x = Time, y = weight)) +
geom_point() +
# 添加线条
geom_smooth(formula = my_formula, method = "lm") +
# 添加相关系数
stat_cor(method = 'spearman', label.y.npc = 0.9) +
# 添加公式
stat_poly_eq(formula = my_formula, aes(label = stat(eq.label)),
label.x.npc = "left", label.y.npc = 0.75) +
theme_pubr()
p1
# 考虑喂养方式
p1 + facet_wrap(vars(Diet))
# 考虑小鸡个体差异
p1 + facet_wrap(vars(Chick))
3-1 先建立全模型
# -----------------------------------------------
# 建立全模型
fit <- lmer(weight ~ Time + (Time | Chick), data = smallSet)
summary(fit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: weight ~ Time + (Time | Chick)
## Data: smallSet
##
## REML criterion at convergence: 560.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.11330 -0.72657 -0.02159 0.61926 2.62272
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Chick (Intercept) 61.20 7.823
## Time 11.48 3.388 -0.95
## Residual 99.73 9.986
## Number of obs: 72, groups: Chick, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 27.741 3.896 5.001 7.120 0.000847 ***
## Time 8.872 1.394 5.000 6.365 0.001415 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Time -0.832
查看模型结果
# -----------------------------------------------
# 查看模型结果
# 限制最大似然 (REML) 方法
# 标准化模型皮尔逊残差的分布情况
quantile(residuals(fit, type = "pearson", scaled = T))
## 0% 25% 50% 75% 100%
## -2.11330258 -0.72657146 -0.02159436 0.61925812 2.62271953
# 随机效应
# 模型对数据的解释能力
# 随机因子随机效应和显著性检验
ranef(fit)
## $Chick
## (Intercept) Time
## 3 0.302098064 -0.6820963
## 1 1.391325144 -1.1669241
## 2 0.002620769 -0.3373429
## 22 6.978608038 -2.6493179
## 23 5.257667357 -1.8413093
## 21 -13.932319371 6.6769905
##
## with conditional variances for "Chick"
ranova(fit)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## weight ~ Time + (Time | Chick)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 -280.05 572.10
## Time in (Time | Chick) 4 -335.84 679.69 111.58 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 固定效应，同一般的线性模型
# 固定效应部分和线性模型结果相似，包括截距估计值和标准误；斜率估计值和标准误。
# 查看固定效应和显著性检验
coef(fit)
## $Chick
## (Intercept) Time
## 3 28.04273 8.190263
## 1 29.13196 7.705435
## 2 27.74326 8.535016
## 22 34.71924 6.223041
## 23 32.99830 7.031050
## 21 13.80831 15.549350
##
## attr(,"class")
## [1] "coef.mer"
anova(fit)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 4040.2 4040.2 1 5.0004 40.511 0.001415 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 查看类型和适用函数
# class(fit)
# methods(class = "lmerModLmerTest")
# confint()函数查看随机效应标准差和固定效应系数的置信区间：
# confint(fit, level = 0.95)
3-2 做减法，生成第二个模型
# ----------------------------------------------------------------
# 随机截距 + 随机斜率
fit.1 <- lmer(weight ~ Time + (Time | Chick), data = smallSet)
# 随机截距 + 固定斜率
fit.2 <- lmer(weight ~ Time + (1 | Chick), data = smallSet)
# 固定截距 + 随机斜率
fit.3 <- lmer(weight ~ Time + (0 + Time| Chick), data = smallSet)
# 固定截距 + 固定斜率
fit.4 <- lm(weight ~ Time, data = smallSet)
# ----------------------------------------------------------------
可视化模型预测结果
分面绘图
# -----------------------------------------------
# 可视化
my_data <- cbind(smallSet, pred = predict(fit))
head(my_data)
## weight Time Chick Diet pred
## 1 42 0 1 1 29.13196
## 2 51 2 1 1 44.54283
## 3 59 4 1 1 59.95370
## 4 64 6 1 1 75.36457
## 5 76 8 1 1 90.77544
## 6 93 10 1 1 106.18631
my_formula <- y ~ x
# Facet by Chick
smallSet %>%
ggplot(aes(x = Time, y = weight)) +
geom_point() +
facet_wrap(vars(Chick)) +
# 添加线条
geom_line(data = my_data, aes(y = pred)) +
# 添加相关系数
stat_cor(data = my_data, aes(y = pred), method = 'spearman', label.y.npc = 0.9) +
# 添加公式
stat_poly_eq(data = my_data, formula = my_formula,
aes(y = pred, label = stat(eq.label)), parse=TRUE,
label.x.npc = "left", label.y.npc = 0.75) +
theme_pubr()
绘图线性混合效应模型四种情况
# as whole plot
DrawPreditLine <- function(model, label_text){
my_data <- cbind(smallSet, pred = predict(model))
my_formula <- y ~ x
smallSet %>%
ggplot(aes(x = Time, y = weight, group = Chick, color = Chick)) +
geom_point() +
# 添加线条
geom_line(data = my_data, aes(y = pred), size = 1) +
labs(title = label_text) +
scale_color_jco() +
theme_pubr() +
theme(legend.position = 'right',
plot.title = element_text(hjust = 0.5))
}
plot_List <- purrr::map2(list(fit.1, fit.2, fit.3, fit.4),
list('随机截距 + 随机斜率', '随机截距 + 固定斜率', '固定截距 + 随机斜率', '固定截距 + 固定斜率'),
DrawPreditLine)
plot_grid(plotlist = plot_List, labels = 'AUTO')
3-3 模型比较，确定最终模型
与全模型fit.1相比，观察到p < 0.05，说明全模型和简化后模型之间存在显著差别，因此，全模型为最终采用的模型。
# ----------------------------------------
# 比较模型
anova(fit.1, fit.2, fit.3)
## Data: smallSet
## Models:
## fit.2: weight ~ Time + (1 | Chick)
## fit.3: weight ~ Time + (0 + Time | Chick)
## fit.1: weight ~ Time + (Time | Chick)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit.2 4 686.53 695.64 -339.26 678.53
## fit.3 4 582.01 591.12 -287.01 574.01 104.517 0
## fit.1 6 577.79 591.45 -282.90 565.79 8.217 2 0.01643 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit.1, fit.2, fit.3, fit.4)
## Data: smallSet
## Models:
## fit.4: weight ~ Time
## fit.2: weight ~ Time + (1 | Chick)
## fit.3: weight ~ Time + (0 + Time | Chick)
## fit.1: weight ~ Time + (Time | Chick)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit.4 3 724.87 731.70 -359.43 718.87
## fit.2 4 686.53 695.64 -339.26 678.53 40.338 1 2.136e-10 ***
## fit.3 4 582.01 591.12 -287.01 574.01 104.517 0
## fit.1 6 577.79 591.45 -282.90 565.79 8.217 2 0.01643 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit.1, fit.2)
## Data: smallSet
## Models:
## fit.2: weight ~ Time + (1 | Chick)
## fit.1: weight ~ Time + (Time | Chick)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit.2 4 686.53 695.64 -339.26 678.53
## fit.1 6 577.79 591.45 -282.90 565.79 112.73 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit.1, fit.3)
## Data: smallSet
## Models:
## fit.3: weight ~ Time + (0 + Time | Chick)
## fit.1: weight ~ Time + (Time | Chick)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit.3 4 582.01 591.12 -287.01 574.01
## fit.1 6 577.79 591.45 -282.90 565.79 8.217 2 0.01643 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit.2, fit.3)
## Data: smallSet
## Models:
## fit.2: weight ~ Time + (1 | Chick)
## fit.3: weight ~ Time + (0 + Time | Chick)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit.2 4 686.53 695.64 -339.26 678.53
## fit.3 4 582.01 591.12 -287.01 574.01 104.52 0
# 与全模型fit.1相比，观察到p < 0.05，说明全模型和简化后模型之间存在显著差别，因此，全模型为最终采用的模型。
参考资料
[1] 线性混合模型（linear mixed model, LME）
https://zhuanlan.zhihu.com/p/477910086
[2] R语言练习的时候那些内置数据集
https://mp.weixin.qq.com/s/HeYb9bL5sTHCjswheh_fzA
[3] R数据分析：混合效应模型的可视化解释，再不懂就真没办法
https://zhuanlan.zhihu.com/p/353317296
[4] Linear Mixed Model：线性混合模型简介
https://blog.csdn.net/weixin_43569478/article/details/108079707
[5] lme4 | 在R中运行混合效应模型（多层模型）
https://zhuanlan.zhihu.com/p/567950808
[6] 混合线性模型的实现（更新20190607）
https://zhuanlan.zhihu.com/p/63092231
[7] R语言线性混合效应模型不同类型的比较
https://mp.weixin.qq.com/s/NPCu4k5tGYodqSntiLMHGg
[8] 线性混合模型（linear mixed model, LME）
https://zhuanlan.zhihu.com/p/477910086
[9] R使用线性混合效应模型一些常见的结果指标
https://zhuanlan.zhihu.com/p/649029208
[10] 线性混合效应模型入门之二 - 实例操作及结果解读
https://zhuanlan.zhihu.com/p/483173133
[11] INTRODUCTION TO LINEAR MIXED MODELS
https://ourcodingclub.github.io/tutorials/mixed-models/?nsukey=JfcXN4ZtaaEgpYJMyxe3oLJVjnlrrISjyCNa8jr%2FJxC5boovrpcf%2BTJJhmqhALfQNznAP1VPeUUSZghPYi93AZHNXDFsaXoP9oL%2FEzSgzpmo9VmN7oUwmYExM%2F5cQ5Y3BhDFhWxHSxbAeT3iYd6fwXLzR65TQThi239DVAazfetJiUJlMdBoxgvfGykUcx50Fnd8BemyGv%2FLMezfhVe0tA%3D%3D&from=timeline#FERE
https://mp.weixin.qq.com/s/7UWieBmszQAHBzzkx1XCiQ