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

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

1.广义线性混合模型（Generalized Linear Mixed Models，GLMMs）在生态学中的应用

热门文章

最新文章

相关课程

相关电子书

相关实验场景