R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码

全文链接 :https://tecdat.cn/?p=35607

原文出处：拓端数据部落公众号

在生态学研究领域，广义线性混合模型（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脚本。

生态学数据如下：

这个数据集包含了一些生态学相关的变量，我们可以逐一解释：

变量解释

Sample：样本编号，用于区分不同的样本或观察值。
Richness：物种丰富度，通常指的是一个特定区域或样本中物种的数量。这是一个常用的生态学指标，用于描述生态系统中生物多样性的一个方面。
Exposure：暴露程度，这指的是样本所在位置接受到的某种环境因素的暴露程度，如阳光、风、雨水等。这个变量可以影响物种的分布和丰富度。
NAP：这个变量代表某种环境参数，如土壤营养状况、盐分含量等。
Beach：海滩或区域的标识，是一个分类变量，用于区分不同的样本是来自哪个海滩或区域。在这个数据集中，有两个不同的值（1和2），可能表示两个不同的海滩或区域。

数据内容

从给出的数据样本中，我们可以看到：

有9个样本被记录，每个样本都有对应的物种丰富度、暴露程度、NAP值和海滩标识。
物种丰富度（Richness）的值在8到19之间，说明不同样本间的生物多样性存在差异。
暴露程度（Exposure）在样本1到样本5中都是10，在样本6到样本9中都是8，这意味着可能有两组不同暴露程度的样本被收集。
NAP值在样本间有所不同，既有正值也有负值，表示不同的环境参数状态或水平。
Beach值表明样本被分为两组，来自两个不同的海滩或区域。

以下是一个R脚本的示例，用于展示如何在广义线性混合模型（GLMM）中演示GLMM的拟合、假设检验、模型选择以及结果推断。我们将使用lme4和arm包进行混合模型的分析，并使用RCurl包来下载示例数据集。

最终得到以下的模型预测图

$output\_55\_0.png$

$glmm\_plot3.png$




    # 展示数据的前几行以确认数据加载正确  

    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只是模型选择的一个方面，还需要考虑其他因素，如模型的假设合理性、解释性等。



# 检查模型假设

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数据集，并对随机效应进行了测试。以下是这段代码的详细解释：



    # 测试随机效应  

    # 第一个模型  

    # 第二个模型没有随机项，使用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)

参数自助法似然比检验：对新的固定效应模型进行了参数自助法似然比检验。



    # 模拟过程  

    for(i in 1:n.sim){  

      # 从零模型中模拟新的观测值  

      dattemp$ysim <- unlist(simulate(mod3_lmer))  

        
      

    # 计算p值  

    (sum(lrt.sim >= lrt.obs) + 1) / (n.sim + 1)

解释：

绘制直方图：通过直方图，我们可以直观地看到模拟的似然比检验统计量的分布情况，并标出观察到的统计量，从而评估其显著性。
固定效应部分的模型选择这部分代码通过拟合不同的混合效应模型来比较固定效应部分的影响。lme函数用于拟合线性混合效应模型，而lmer函数用于拟合线性混合效应模型，但使用的是lme4包。两种模型都考虑了随机效应，但使用了不同的估计方法（lme使用最大似然法，lmer使用REML或ML）。
参数自助法似然比检验：与之前的自助法类似，但这次是针对固定效应部分的模型进行比较。代码从mod3_lmer（只包含NAP作为固定效应的模型）中模拟新的观测值，然后拟合零模型和替代模型，并计算似然比检验统计量。最后，基于模拟的统计量计算p值，以评估固定效应Exposure是否显著。

注意：在代码中，simulate函数用于从模型生成模拟数据，而anova函数用于比较模型的差异。此外，unlist函数用于将列表转换为向量，因为simulate函数返回的可能是一个列表

R复制代码
    # 使用蓝色绘制直方图，直方图的x轴范围为0到lrt.sim和lrt.obs中的最大值，并设置x轴和y轴的标签大小  

      

    # 在直方图上添加一条垂直于x轴的橙色线，表示lrt.obs的值，线宽为3  

    abline(v=lrt.obs, col="orange", lwd=3)  

      

    # 下一步将是先删除NAP，然后看似然比检验是否显著，以及是否总是先删除Exposure导致更高的LRT统计量  

    # 其他方法，如AIC...  

    # GLMM的R平方计算，参见Nakagawa 2013 MEE的补充材料  

    # 计算固定效应方差  



    # VarCorr()函数用于提取方差分量  

    # attr(VarCorr(lmer.model),’sc’)^2提取残差方差，VarCorr()$plot提取plot效应的方差  

    # 计算条件R平方  

    #conditionnal R-square 的计算公式

解释：

hist函数用于绘制直方图，展示lrt.sim（可能是似然比检验的统计量模拟值）的分布情况。直方图的x轴范围设置为从0到lrt.sim和lrt.obs（观察到的似然比检验统计量）中的最大值。同时设置了直方图的颜色、x轴和y轴的标签以及标签的大小。
abline函数在直方图上添加了一条垂直于x轴的线，线的位置为lrt.obs的值，线的颜色为橙色，线宽为3。这通常用于在直方图上标识某个特定的观察值或阈值。
注释部分说明了下一步的分析计划，即先删除NAP变量，然后检查似然比检验是否显著，以及比较先删除NAP还是先删除Exposure变量对LRT统计量的影响。同时提到了其他分析方法，如AIC（赤池信息准则）。
接下来的代码计算了线性混合效应模型mod1_lmer的条件R平方。这包括计算固定效应的方差（VarF），提取模型的方差分量（VarCorr），以及计算条件R平方的值。条件R平方是解释模型中固定效应和随机效应共同解释的变异比例。

R复制代码
    # 从模型中推断  

    # lme 和 glmer 可以获取 p 值，但 lmer 不行  



      

    # 使用 glmer 拟合模型  


    summary(mod1_glmer)  

      

    # 使用 arm 包中的 sim 函数  

    n.sim <- 1000  


    head(simu@fixef)  

    # 95% 可信区间  



    # 绘制 NAP 对丰富度的影响  



    predmat <- matrix(ncol = nsim, nrow = nrow(newdat))  

      


    lines(newdat$NAP, newdat$lower, col = "red", lty = 2, lwd = 1.5)  

      

    # 比较不同项的系数，标准化变量  

    data$stdNAP <- scale(data$NAP)  



    # 95% 可信区间  

    
      

    # 绘图  

    plot(1:3, coeff[, 2], ylim = c(-0.8, 2), xaxt = "n", xlab = "", ylab = "Estimated values")

代码解释

这段代码的主要目的是从广义线性混合效应模型（GLMM）中推断结果，并绘制相关效应的可视化图形。

模型推断：
- summary(mod1) 和 summary(mod1_lmer) 用于展示模型 mod1 和 mod1_lmer 的统计摘要，包括系数估计、标准误差、z值或t值以及对应的p值。
拟合GLMM模型：
- glmer 函数用于拟合广义线性混合效应模型，这里以物种丰富度（Richness）作为响应变量，NAP和Exposure作为固定效应，同时考虑到Beach作为随机效应。
模拟后验分布

2.生态学模拟对广义线性混合模型GLMM进行功率（功效、效能、效力）分析power analysis环境监测数据|附代码数据

假设检验的功效定义为假设原假设为假，检验拒绝原假设的概率。换句话说，如果一个效应是真实的，那么分析判断该效应具有统计显着性的概率是多少？

概括

r 语言允许用户计算 lme 4 包中广义线性混合模型的功效。功率计算基于蒙特卡罗模拟。
它包括用于 (i) 对给定模型和设计进行功效分析的工具；(ii) 计算功效曲线以评估功效和样本量之间的权衡。
本文提供了一个教程，使用具有混合效果的计数数据的简单示例（具有代表环境监测数据的结构）。

介绍

如果一项研究的功效不足，资源可能被浪费，真正的效果可能被遗漏。另一方面，一项大型研究的花费可能过大，因此其费用也会超过必要的范围。因此，在收集数据之前进行功效分析是一个很好的做法，以确保样本具有适当的规模来回答正在考虑的任何研究问题。

广义线性混合模型 (GLMM) 在生态学中很重要，它允许分析计数和比例以及连续数据，并控制空间非独立性.

蒙特卡罗模拟是一种灵活且准确的方法，适用于现实的生态研究设计。在某些情况下，我们可以使用解析公式来计算功效，但这些通常是近似值或需要特殊形式的设计。仿真是一种适用于各种模型和方法的单一方法。即使公式可用于特定模型和设计，定位和应用适当的公式也可能非常困难，因此首选仿真。

对于对 r 不够熟悉的研究人员，设置模拟实验可能太复杂了。在本文中，我们介绍了一个工具来自动化这个过程。

r 包

有一系列的 r 包目前可用于混合模型的功效分析。然而，没有一个可以同时处理非正态因变量和广泛的固定和随机效应规范。

图1

r 旨在与任何可以与 lme 4 中的 lmer 或 glmer 配合的线性混合模型 (LMM) 或 GLMM 一起使用。这允许具有不同固定和随机效应规范的各种模型。还支持在 r 中使用 lm 和 glm 的线性模型和广义线性模型，以允许没有随机效应的模型。

r 中的功效分析从适合 lme 4 的模型开始。

在 r 中，通过重复以下三个步骤来计算功效：(i) 使用提供的模型模拟因变量的新值；(ii) 将模型重新拟合为模拟因变量；(iii) 对模拟拟合应用统计检验。在此设置中，已知存在测试效果，因此每个阳性测试都是真正的阳性，每个阴性测试都是 II 类错误。可以根据步骤 3 的成功和失败次数计算测试的功效。

教程

本教程使用包含的数据集。该数据集代表环境监测数据，在连续固定效应变量\_x \_（例如研究年份）的10 个水平上测量三个组 \_g \_（例如研究地点）的因变量 \_z \_（例如鸟类丰度）。还有一个连续因变量 \_y \_，在本教程中没有使用。

拟合模型

我们首先将 lme 4 中的一个非常简单的泊松混合效应模型拟合到数据集。在这种情况下，我们有一个随机截距模型，其中每个组 ( g ) 都有自己的截距，但这些组共享一个共同的趋势。

glm
summary

本教程重点介绍关于\_x \_趋势的推断。在这种情况下，\_x \_的估计效应大小为 -0.11，使用默认\_z\_检验在 0.01 水平上显着。

请注意，我们特意使用了一个非常简单的模型来使本文易于理解。例如，适当的分析会包含更多的组，并会考虑过度分散等问题。

简单的功率分析

假设我们想重复这项研究。如果效果是真实的，我们是否有足够的功效来期待积极的结果？

指定效应量

在开始功效分析之前，重要的是要考虑您感兴趣的效果大小类型。功效通常随效果大小而增加，较大的效果更容易检测。回顾性“观察功效”计算，其中目标效应大小来自数据，给出误导性结果.

对于此示例，我们将考虑检测 -0.05 斜率的功效。可以使用 lme 4 函数拟合 glmer 模型中的固定效应。然后可以更改固定效应的大小。变量\_x \_的固定效应的大小可以从 -0.11 更改为 -0.05，如下所示：

fixe<‐ ‐0.05

在本教程中，我们只更改变量\_x \_的固定斜率。但是，我们也可以更改随机效应参数或残差方差（适用于合适的模型）。

运行功效分析

一旦指定了模型和效应大小，在 r 中进行功效分析就非常容易了。由于这些计算基于蒙特卡罗模拟，因此您的结果可能略有不同。如果你想得到和本文一样的结果，你可以使用 set.seed(123)。

power

鉴于此特定设置，拒绝\_x \_中零趋势的零假设的能力约为 33%。这几乎总是被认为是不够的；传统上，80% 的功率被认为是足够的.

在实践中， z\_检验可能不适合这样一个小例子。参数引导测试可能是最终分析的首选。但是，更快的 \_z -test 更适合学习使用该包以及在功效分析期间进行初始探索性工作。

增加样本量

在第一个示例中，估计功率很低。小型试点研究通常没有足够的功效来检测微小的影响，但更大的研究可能会。

试点研究对\_x 的 \_10 个值进行了观察，例如代表研究第 1 年到第 10 年。在此步骤中，我们将计算将其增加到 20 年的影响。

modl2 <‐ extend
power(modl2)

沿参数指定要扩展的变量，n 指定要替换它的级别。扩展模型 2 现在将具有从 1 到 20 的\_x \_值，与以前一样分为三组，总共 60 行（与模型 1 中的 30 行相比）。

通过观察\_x 的 \_20 个值，我们将有足够的能力来检测大小为 -0.05 的效应。

各种样本量的功效分析

当数据收集成本高昂时，用户可能只想收集达到一定统计能力所需的数据量。功效曲线函数可用于探索样本大小和功效之间的权衡。

确定所需的最小样本量

在前面的示例中，当对变量\_x 的\_20 个值进行观察时，我们发现了非常高的 \_功效 \_。我们能否减少这个数字，同时保持我们的功效高于通常的 80% 阈值？

 poerCure

print
plot

请注意，我们已将此结果保存到变量 pc2 以匹配模型 2 中的编号。由于模型 1 没有足够的功率，我们没有通过 powerCurve 运行它。绘制的输出如图所示。我们可以看到，检测\_x \_趋势的能力随着采样大小的增加而增加。这里的结果基于将模型拟合到 10 个不同的自动选择的子集。最小的子集仅使用前 3 年（即 9 个观测值），最大的子集使用所有 20 个假设研究年份（即 60 行数据）。该分析表明，该研究必须运行 16 年才能有≥80% 的功效来检测指定大小的影响。