在生态学研究领域，广义线性混合模型（Generalized Linear Mixed Models，简称GLMMs）是一种强大的统计工具，能够同时处理固定效应和随机效应，从而更准确地揭示生态系统中复杂关系的本质（点击文末“阅读原文”获取完整代码数据）。
随着数据分析技术的不断发展，R语言已成为生态学家们进行数据分析的首选工具之一，而GLMMs在R语言中的实现与应用也日益受到关注。
相关视频
本文旨在通过2个实例，帮助客户展示R语言中广义线性混合模型在生态学中的应用及其可视化方法。
1.广义线性混合模型（Generalized Linear Mixed Models，GLMMs）在生态学中的应用
广义线性混合模型（Generalized Linear Mixed Models，GLMMs）在生态学中的应用以及如何在R中实现它们是一个广泛且深入的主题。这类模型是线性混合模型（Linear Mixed Models，LMMs）和广义线性模型（Generalized Linear Models，GLMs）的自然延伸，它们提供了在存在随机效应时，对响应变量与预测变量之间关系进行建模的灵活性。
在生态学中，广义线性混合模型的应用特别广泛，这主要得益于它们在处理复杂数据结构、考虑随机效应以及解释非线性关系方面的能力。例如，在种群动态、群落组成、物种分布等研究中，广义线性混合模型经常被用来解释和预测各种生态现象。
这篇文章主要是为了展示如何拟合GLMM、如何评估GLMM假设、何时在固定效应模型和混合效应模型之间做出选择、如何在GLMM中进行模型选择以及如何从GLMM中得出推论的R脚本。
使用数据（查看文末了解数据免费获取方式）如下：
以下是一个R脚本的示例，用于展示如何在广义线性混合模型（GLMM）中演示GLMM的拟合、假设检验、模型选择以及结果推断。我们将使用lme4和arm包进行混合模型的分析，并使用RCurl包来下载示例数据集。
最终得到以下的模型预测图
# 展示数据的前几行以确认数据加载正确 head(data) # 第一部分：拟合GLMM # 随机截距模型 summary(mod_lmer1) # 查看模型摘要 # 随机斜率和截距模型 summary(mod_lmer2) # 查看模型摘要 # 第二部分：GLMM假设检验 # 检查随机效应的方差成分是否显著 Anova(mod_lmer1, type="II Wald") # 使用Anova函数进行方差分析 # 检查残差的正态性、同方差性等假设 plot(mod_lmer1) # 绘制模型诊断图 # 第三部分：模型选择 # 使用AIC进行模型选择 AIC(mod_lmer1, mod_lmer2) # 第四部分：从GLMM中得出推论 # 获取固定效应的系数估计和置信区间 # 获取随机效应的方差估计 VarCorr(mod_lmer2) # 绘制模型预测图 library(ggplot2) ggplot
geom_smooth函数在ggplot2中默认不支持lmer模型，你可能需要手动计算预测值并添加到数据框中，或者使用其他包（如ggeffects或effects）来生成预测值并绘制图形。
另外，关于嵌套和交叉随机效应的问题，lme4包中的lmer函数支持拟合这些复杂的随机效应结构。你可以通过在公式中指定适当的随机效应项来实现这一点。例如，对于嵌套随机效应，你可以使用(1|Group1/Group2)；对于交叉随机效应，你可以使用(1|Group1) + (1|Group2)。具体的随机效应结构取决于你的研究设计和假设。
展示了如何模拟固定效应和混合效应模型的数据，并比较了这两种模型的拟合结果。您的目标是展示固定效应模型和混合效应模型之间的区别，并通过绘图来可视化这些差异。然而，从您提供的模拟数据和结果来看，没有明显的差异是可见的，这可能是因为模拟的数据没有表现出强烈的随机效应。
添加一些额外的解释和可视化步骤，以帮助更好地理解固定效应和混合效应模型之间的区别。
plot(fitted(m_lme), resid(m_lme)) qqnorm(resid(m_lme)) qqline(resid(m_lme)) # 添加参考线 qqnorm(ranef(m_lme)[, 1]) qqline(ranef(m_lme)[, 1]) # 添加参考线 plot(x, resid(m_lme)) # 固定效应模型拟合图 pal <- brewer.pal(6, "Set1") plot(y ~ x, col = pal[f], pch = 16, main = "Fixed Effects Model") # 混合效应模型拟合图 xyplot(y ~ x | f # 比较固定效应和混合效应的拟合结果 # 可以计算模型的AIC、BIC等指标，或者通过交叉验证来评估模型性能 AIC(m_lm) AIC(m_lme) # 注意：混合效应模型通常更适用于存在随机效应的情况， # 如果随机效应不显著，固定效应模型可能更简洁且易于解释。 # 结论： # 在这个模拟例子中，固定效应和混合效应模型的拟合结果没有明显的视觉差异， # 这可能是因为随机效应的影响较小。在实际应用中，需要根据具体数据和问题来选择适当的模型。 # 通过混合效应技术生成随机斜率和截距模型 # 检查模型假设 par(mfrow=c(2,2)) plot(m_lm2) plot(fitted(m_lme2), resid(m_lme2)) abline(h=0, lty=2, col="red") qqnorm(resid(m_lme2)) qqnorm(ranef(m_lme2)[,1]) qqnorm(ranef(m_lme2)[,2]) p for(i in 1:length(levels(f))){ lines(x, fixef(m_lme2)[1] + (fixef(m_lme2)[2] + ranef(m_lme2)[i,2])*x + ranef(m_lme2)[i,1], col=pal[i], lwd=1.5) }
在这个修订后的脚本中，我添加了qqline函数来在QQ图上绘制参考线，以便更清晰地查看残差和随机效应的正态性。我还使用了lattice包的xyplot函数来绘制混合效应模型的拟合图，其中每个组（f）的拟合线被单独绘制。
请注意，为了清楚地看到固定效应和混合效应模型之间的差异，您可能需要模拟更强的随机效应，或者在实际数据集上应用这些模型，这些数据集通常包含更复杂的结构和随机性。最后，我还添加了AIC值的计算，这是一个常见的模型选择指标。通过比较不同模型的AIC值，您可以获得关于哪个模型更适合数据的额外信息。然而，请注意，AIC只是模型选择的一个方面，还需要考虑其他因素，如模型的假设合理性、解释性等。
Kaizong Ye
分析师
检查模型假设 qqline(ranef(m_lme2)[,2]) # 在qqnorm图上添加参考线 scatter.smooth(fitted(m_lme2), sqrt(abs(resid(m_lme2)))) # 绘制拟合值与残差绝对值的平方根的散点平滑图 # 错误数据 m_wrg <- lme(y ~ x, random = ~x|f) # 使用错误数据构建混合效应模型 scatter.smooth(fitted(m_wrg), sqrt(abs(resid(m_wrg)))) # 绘制拟合值与残差绝对值的平方根的散点平滑图 # 绘制拟合值与残差的关系图，对残差和所有随机效应进行qqnorm检验
- 在qqnorm图上添加qqline可以更容易地判断数据是否符合正态分布。
- scatter.smooth函数用于绘制散点图并添加平滑曲线，用于观察变量之间的关系。
- 在实践2中，我故意制造了一些错误数据，用来展示当数据不符合模型假设时，混合效应模型的表现。通过比较正确数据和错误数据的模型结果，可以更好地理解模型假设的重要性。
这段代码主要是进行模型选择，它使用了RIKZ数据集，并对随机效应进行了测试。以下是这段代码的详细解释：
# 实践2测试随机效应 # 第一个模型 # 第二个模型没有随机项，使用gls是因为它也可以通过REML拟合模型 # 似然比检验，对于小样本量来说可能不够精确 anova(mod1, mod2) # 参数自助法 lrt.sim <- numeric(n.sim) dattemp <- data for(i in 1:n.sim){ # 从零模型中模拟新的观测值 dattemp$ysim <- simulate(lm(Richness ~ NAP + Exposure, data = dattemp))$sim_1 # 保存似然比检验统计量 lrt.sim[i] <- anova(modnullsim, modaltsim)$L.Ratio[2] } # 计算p值 (sum(lrt.sim >= lrt.obs) + 1) / (n.sim + 1)
解释：
- 读取数据：从指定路径读取RIKZ数据集，数据由空格分隔，并且包含表头。
- 测试随机效应：
mod1：使用lme函数拟合一个混合效应模型，其中Richness（丰富度）是响应变量，NAP和Exposure是固定效应，Beach是随机效应的分组变量。mod2：使用gls函数拟合一个广义最小二乘模型，该模型没有随机效应。
- 似然比检验：使用
anova函数比较两个模型，但请注意，对于小样本量，似然比检验可能不够精确。 - 参数自助法：这是一种估计模型选择检验p值的方法，通过模拟数据来估计检验统计量的分布。
- 从零模型中模拟新的观测值。
- 拟合零模型和替代模型。
- 保存似然比检验统计量。
lrt.obs：保存观察到的似然比检验统计量。- 进行1000次模拟，每次：
- 使用模拟的似然比检验统计量来估计p值。
最终，代码返回了一个p值，该值基于参数自助法估计，用于评估随机效应是否显著。
这段代码主要进行了以下操作：
- 绘制直方图：绘制了模拟的似然比检验统计量的直方图，并在图上标出了观察到的似然比检验统计量。
# 绘制直方图 par(mfrow=c(1,1)) # 设置图形参数，使接下来的图形在一个单独的图形窗口中显示 hist(lrt.sim, xlim=c(0, max(c(lrt.sim, lrt.obs))), col="blue", abline(v=lrt.obs, col="orange", lwd=3) # 在直方图上标出观察到的似然比检验统计量
- 固定效应部分的模型选择：使用
lme和lmer函数拟合不同固定效应的混合效应模型，并比较这些模型。
# 使用最大似然法(ML)拟合混合效应模型 # 使用lmer函数拟合混合效应模型 # 显示模型摘要 summary(mod1_lmer) summary(mod1_ML) # 使用anova函数比较模型 anova(mod1_lmer, mod3_lmer)
R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码2：https://developer.aliyun.com/article/1501233
