R统计绘图-线性混合效应模型详解(理论、模型构建、检验、选择、方差分解及结果可视化)

���录

一、 基础理论

二、数据准备

三、构建线性混合效应模型(LMMs)

3.1 lme4线性混合效应模型formula

3.2 随机截距模型构建及检验

3.3 随机截距模型分析结果解释及可视化

3.4 随机斜率模型构建、检验及可视化

四、线性混合效应模型选择

4.1 多模型比较

4.2 模型最优子集筛选

一、 基础理论

各个数据模型都包含一些假设，当数据满足这些假设时，使用这些模型得到的结果才会更准确。比如广义线性模型要满足“线性”和“独立性”等条件。为了适应更多数据类型，开发出了各种模型，可根据数据情况，选择合适的模型。比如，如果数据不满足“线性”，则可考虑广义可加模型；如果不满足“独立性”，则可以考虑混合效应模型(mixed effect model)。

多水平模型同时包含多个分类因子数据，最主要特点就是非独立性，在多个水平上都存在残差，因此适合混合效应模型，总的来说，其思想就是把高水平上的差异估计出来，这就使得残差变小，估计的结果更为可靠。混合效应模型又称随机效应模型(Random Effect model)，分层线性模型，随机系数模型，方差成分模型等。

比如，一份包含不同施肥处理的多土层微生物数据，其中不同施肥处理与不同土壤深度即是嵌套数据(多水平数据)。

如果用线性模型分析不同施肥处理间微生物群落的差异，则没有考虑到不同土壤深度的微生物群落的差异，暗含了假设：不同土壤深度的微生物群落随不同施肥处理变化的截距与斜率都是相同的。不同土壤深度土壤的微生物群落差异就被线性模型归类到误差中去了。

如何解决这一问题呢？其中一种方法就是使用固定效应模型(Fixed Effect Model)，即将不同土壤深度的微生物群落结构随不同施肥处理变化的截距差异和斜率差异都估计出来。在土壤深度分类较少时可以这么做，将次级分类水平作为虚拟变量回归，3分类，则估计2个截距差异和斜率差异。当次级分类水平过多时，需要估计的参数太多，用虚拟变量就会消耗太多的自由度，估计结果不可靠，此时即可考虑使用混合效应模型。

固定效应与随机效应因子之间没有明显的分类，但显然随机效应因子要是分类数据，且数据跟随机效应存在相关性。有文章推荐随机效应的分类水平最好大于5类，因为当分类水平较少时，模型对方差的估计不太准确。当随机因子大于1个时，随机因子间的关系也需要注意，根据随机因子间关系，可能包括交叉(Crossed)、部分交叉或嵌套(nested)。

混合效应模型也分为线性混合效应模型、广义线性混合效应模型和非线性混合效应模型。此文章介绍使用R中lme4包进行线性混合效应模型分析。

线性混合效应模型假设：

1. 残差应该是正态分布的-可以适当放宽，前提是预测值与残差没有相关性，

2. 残差的方差是同质的-非常重要的假设；

3. 随机效应的观测是相互独立的；

4. 自变量与因变量存在线性关系；

5. 自变量间不存在相关性-方差膨胀因子筛选或者数据降维。

当然，对数据中异常点的检测也是很有必要的。

 二、数据准备

此套数据是嵌套数据，包含分类因子depth和tillage。嵌套数据中包含多个分类自变量，数据是非独立性的，模型构建时考虑使用混合效应模型。选择pH为因变量。

​# 设置工作路径#knitr::opts_knit$set(root.dir="D:\EcoEvoPhylo\MEM\LMM")# 使用Rmarkdown进行程序运行Sys.setlocale('LC_ALL','C') # Rmarkdown全局设置#options(stringsAsFactors=F)# R中环境变量设置，防止字符型变量转换为因子library(tidyverse)library(rstatix)library(car)library(lme4)library(lmerTest) # 模型固定效应系数检验library(multilevelTools) # 模型诊断检验library(extraoperators)library(JWileymisc) # 模型诊断及APAStyler以表格形式输出模型结果library(effectsize) # 计算标准化回归系数library(influence.ME) # 检测高影响观测点library(GGally)​# 2.1 导入数据env % head​# 2.2 计算均值与标准误env %>%  group_by(condition) %>%  get_summary_stats(names(env[4:ncol(env)]),type="mean_se") %>%  left_join(env[1:3],by = "condition") -> datadata %>% sample_n_by(tillage,size=2)

图1|原始数据表，env.csv。

图2|均值-标准误差计算结果，data。

# 2.3 数据探索-查看数据是否适合进行线性混合效应模型分析## 2.3.1 ICC-组内相关系数计算icc1  p1) # 绘图函数加上括号，赋值给变量的同时，会在屏幕上显示。​### 输出pdf图到本地ggsave("p1.pdf", p1, device = "pdf",width = unit(10,"cm"),height = unit(8,"cm"))

图6|long format格式数据表，data.p。

图7|可视化pH数据分布，p1.pdf。

## 2.3.5 绘制变量相关性图### 设置绘图颜色cols % group_by(`tillage`) %>% select(pH) %>% rstatix::get_summary_stats() %>% data.frame()

​

图29|估计边际均值结果，emmeans(fm1, ~tillage)。估计边际均值结果中的emmean列表示固定因子变量各水平以及固定定量变量的pH的最小二乘均值，tillage的结果中WL的emmean值最大，表示pH在NT和WL的值更大，pH均值在两者中也是更大的，只有些微差别。结果解释可参考：https://www.zhihu.com/question/370213292。

# 3.3.3  Post-hoc / Contrast Analysislibrary(knitr)## tillage水平间估计边际均值的两两比较emmeans::emmeans(fm1, pairwise ~ tillage, adjust = "bonferroni")$contrasts %>%   tidy() %>%   mutate_if(is.numeric, ~round(., 4)) %>%   kable()​## tillage与AP交互效应的估计边际均值的两两比较emmeans::emmeans(fm1, pairwise ~ tillage * AP, adjust = "bonferroni")$contrasts %>%   tidy() %>%   mutate_if(is.numeric, ~round(., 4)) ​## 估计边际均值事后检验结果绘图emmeans::pwpp(emmeans::emmeans(fm1, pairwise ~ tillage * AP, adjust = "bonferroni"), type = "response")+theme_minimal()

​

图30|tillage水平间估计边际均值的两两比较。

图31|tillage水平间估计边际均值的两两比较。

图32|估计边际均值事后检验结果绘图。

# 3.3.4 方差分解-不同方法计算出的R2结果会有差异，但是趋势都是一致的。## 方法一：R2计算模型整体的MarginalR2和ConditionalR2 R2(fm1)​## 方法二：modelTest()会计算模型整体和固定效应的MarginalR2和ConditionalR2。fm1.da %  select(Model,MarginalR2,ConditionalR2)### 固定效应的MarginalR2和ConditionalR2。fm1.da$EffectSizes %>%  select(Variable,MarginalR2,ConditionalR2)​## 方法三：r2()会计算模型整体的Nakagawa's MarginalR2和ConditionalR2。performance::r2(fm1,version=TRUE)​## 方法四：MuMln包的r.squaredGLMM()只提供模型整体的marginal R2和conditional R2。# install.packages("D://software//MuMIn_1.47.8.zip",'repos = NULL')MuMIn::r.squaredGLMM(fm1) #计算模型的边际及条件R方,程辑包‘MuMIn’需要>= 4.2.0。​## 方差分解：glmm.hp包的glmm.hp():"hierarchical partitioning"glmm.hp::glmm.hp(fm1)$hierarchical.partitioning​## 方差分解结果绘图### 方差分解绘图数据提取r2  var.decom)ggsave("var.decom.pdf",var.decom,device = "pdf",family="Times",width = 7,height = 3)

​

图33|方差分解结果表。此处选择的glmm.hp包的glmm.hp()方差分解的结果。hierarchical partitioning法的方差分解结果，Unique表示固定因子的解释的方差，与fm1.da$EffectSizes中计算的固定效应的边际R2值相近。Average.share表示层次分割计算的共享方差，Inpidual表示固定因子的单独MarginalR2，I.perc(%)是固定因子独自解释因变量的方差占模型MarginalR2的比例。r.squaredGLMM的方差与R2的计算结果有些微差异，但是趋势是一致的。方差分解有很多算法，选用方法的时候应该弄清楚背后的原理，理论文章是要引用的。

图34|方差分解结果绘图。柱子中的星号表示固定效应参数检验显著。

# 3.3.5 绘制效应图(effect plot)## 绘制随机效应水平间变化：绘制随机截距和斜率与模型总体截距和斜率的差异。lattice::dotplot(ranef(fm1, condVar=T))​## 绘制固定效应预测效应图。plot(effects::allEffects(fm1)) ​## 绘制模型预测效果图。performance::check_predictions(fm1)

​

图35|随机效应系数散点图。随机截距模型只有截距的随机效应参数。

图36|模型固定效应预测效应图。AP效应图横轴上的小竖线表示数据点。

图37|模型预测效果图。

# 3.3.6 输出模型文字描述report::report(fm1)

​

We fitted a linear mixed model (estimated using ML and nloptwrap optimizer) to predict pH with tillage and AP (formula: pH ~ tillage + AP). The model included depth as random effect (formula: ~1 | depth). The model's total explanatory power is substantial (conditional R2 =

0.88) and the part related to the fixed effects alone (marginal R2) is of 0.86. The model's intercept, corresponding to tillage = CK and AP = 0, is at 5.83 (95% CI [5.65, 6.01], t(29) = 67.05, p

- The effect of tillage [NT] is statistically significant and positive (beta = 0.35, 95% CI [0.22, 0.49], t(29) = 5.42, p

- The effect of tillage [PT] is statistically non-significant and positive (beta = 0.08, 95% CI [-0.05, 0.22], t(29) = 1.25, p = 0.223; Std. beta = 0.20, 95% CI [-0.13, 0.54])

- The effect of tillage [WL] is statistically significant and positive (beta = 0.38, 95% CI [0.25, 0.51], t(29) = 6.04, p

- The effect of AP is statistically significant and negative (beta = -0.02, 95% CI [-0.03, -0.02], t(29) = -10.61, p

Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation.

可参照此输出描述自己的结果。

3.4 随机斜率模型构建、检验及可视化

在文章中需要描述以下信息：1)模型描述(因变量、固定效应和随机效应)；2)数据是否满足线性混合模型假设；3)数据是否存在异常值及如何处理异常值；4)模型系数与模型检验结果；5)根据研究目的，有些还需要自变量的MarginalR2分解结果。

理论知识在随机截距模型部分已提及，此处用随机斜率模型演示如何快速得到上述结果。

此处因为fm1使用的ML法，为了方便后续比较，后续模型都使用的ML法，但是ML法可能会低估随机效应的方差，所以模型构建可优先选用REML法。特别当混合效应模型存在嵌套随机效应时，则不应该选用ML算法。

# 3.4.1 构建随机斜率模型control