R语言+AI提示词：贝叶斯广义线性混合效应模型GLMM生物学Meta分析

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

本文旨在帮助0基础或只有简单编程基础的研究学者，通过 AI 的提示词工程，使用 R 语言完成元分析，包括数据处理、模型构建、评估以及结果解读等步骤（点击文末“阅读原文”获取完整代码、数据、文档）。

主要介绍了元分析的概念、工作原理、固定效应与随机效应元分析的区别，以及模型的语法、输出结果解读、先验知识的应用等内容，还探讨了如何控制测量误差和处理不同的方差结构。

一、引言

元分析是对多个相似主题的单个研究结果进行的统计分析，它能提供比单个研究更可靠的估计，还能揭示研究间的模式和趋势。在生物研究中，我们常常寻找生物对不同处理或环境响应的预测因子，元分析是实现这一目标的有效方法。
在贝叶斯统计框架下，使用马尔可夫链蒙特卡罗（MCMC）方法拟合广义线性混合效应模型（GLMM）。即使不完全理解贝叶斯统计和混合建模的复杂概念，也可以使用该包进行分析。

二、元分析的概念

元分析是对众多关于相似主题的单个研究结果的统计分析。与单个研究相比，它能提供更稳健的估计，还能在控制单个研究中固有的非独立性和测量误差来源的同时，揭示研究间的模式和趋势。
例如，当比较来自不同位置（如纬度、海拔、半球、气候区）、不同物种（如具有不同行为或生活史特征）或不同时间段（如研究的时间和持续时间）的研究时，会引入非独立性来源，在估计所有研究的平均效应时需要控制这些来源。
然而，这些非独立性来源也可能是我们感兴趣的。比如，在控制纬度时，可能会发现它能解释我们所关注的响应中研究间的很大一部分方差，从而可以说纬度是该响应的一个良好预测因子。
作为生物学家，我们经常寻找生物对不同处理或环境等响应的预测因子（如上述提到的位置差异、物种或时间段），元分析是实现这一目标的好方法。

三、MCMCglmm介绍

MCMCglmm使用马尔可夫链蒙特卡罗（MCMC）方法，在贝叶斯统计框架下拟合广义线性混合效应模型（GLMM）。贝叶斯统计听起来可能很复杂，但实际上比频率统计更直观。
以抛硬币为例，频率统计认为抛硬币出现正面的概率是0.5，这是在大量抛硬币实验中正面出现的频率，即基于一组参数值观察数据的概率。而贝叶斯统计认为，如果硬币已经抛出但未让你看到结果，你说出现正面的概率是0.5，是因为你不知道实际结果，硬币要么是正面要么是反面，概率为0或1，贝叶斯统计的基础是存在一个真实的参数值，但你不知道它是什么，它考虑的是基于观察数据的一组参数值的概率，表征了对真实值的主观不确定性。
在贝叶斯统计中，我们基于对先前情况的了解，在模型中纳入先验概率。此时，数据是固定的，而参数根据我们的先验知识以及我们对某种结果发生可能性的判断而改变。模型的输出是一个后验分布，它是数据、先验知识和似然函数的组合。

四、固定效应与随机效应元分析

在贝叶斯分析中，固定效应和随机效应没有根本区别，关键在于理解每种类型的分析如何处理方差。
固定效应元分析假设任何观察间的方差仅由抽样误差引起（即当精度趋于无穷大时，漏斗图会收敛到一个点）。在组合实验研究时特别有用，例如在分析大量个体处理的元分析中，如果我们想明确回答“雄性或雌性对某种处理的反应更有效吗”这样的问题，可以将性别作为固定效应。当我们将一个效应视为固定效应时，认为其他效应的知识不会提供关于我们关注效应的可能大小的信息。
随机效应元分析则会利用这些信息，认为抽样误差可能导致极端估计，因此这些估计应被给予较少的统计权重。此外，观察间的一些方差被认为是真实的并被估计。在处理来自野外系统的数据时特别有用，因为我们假设存在自然变异，并希望找到有助于解释它的模式。
例如，在研究候鸟到达繁殖地的时间是否与过去几十年大致相同，或者由于气候变化现在是更早还是更晚到达的问题中，我们收集了来自许多不同已发表研究的数据，这些数据包含了不同物种、国家、纬度等信息。我们可以使用包含这些信息的数据集，通过元分析计算所有种群的平均天数差异，同时将物种作为随机效应，以了解物种间的方差，并估计每个物种的平均响应。

五、熟悉MCMCglmm的语法和模型输出

准备工作

安装软件：确保你已经安装了 R 语言环境。你可以从 R 官方网站（https://www.r-project.org/）下载并安装适合你操作系统的版本。
安装相关包：打开 R 语言环境，运行以下代码安装本教程所需的包：
数据准备：准备好你的数据文件，本教程中使用的示例数据文件名为“metadata.csv”，请确保数据文件的格式正确且包含所需的变量（如 Predictor、Slope、SE、Species、Location、Study 等），并将其放置在你指定的工作目录中。

AI 工具选择

你可以选择使用各种支持自然语言生成代码的 AI 工具，如 deepseek、ChatGPT等。在本教程中，我们将以通用的方式介绍提示词工程，你可以根据自己的喜好选择合适的 AI 工具。

提示词工程步骤

启动 AI 工具并输入基本提示：打开你选择的 AI 工具，输入以下提示词，让 AI 了解我们的任务背景：
“我是一个研究学者，想使用 R 语言进行元分析，数据已经准备好，存储在名为‘metadata.csv’的文件中，数据包含 Predictor、Slope、SE、Species、Location、Study 等变量，你能逐步指导我完成这个元分析过程吗？”

数据导入与初步查看：根据 AI 的回复，你会得到类似以下的代码及解释，按照提示在 R 语言中运行代码：

在进行分析时，首先需要下载数据并导入到R中，加载所需的包。通过设置工作目录来指定数据所在的文件夹，可以使用setwd("your-filepath")代码（将your-filepath替换为实际文件路径），或者通过点击Session/Set working directory/Choose directory来选择文件夹。

# 加载包
library("dplyr") # 用于数据操作
migrata <- read.csv("metadata.csv", header = T) # 导入数据集
View(migrata) # 查看数据集。检查Predictor变量。有两个，时间和温度

查看数据集后，为了专注于研究的第一部分，我们可以筛选出Predictor为“year”的行，即仅关注随时间测量的年度到达日期的数据。

此时，AI 会解释说，library()函数用于加载所需的包，read.csv()函数用于读取 CSV 格式的数据文件，View()函数可以让你在数据查看器中直观地查看数据内容。

数据筛选：继续向 AI 提问，例如：“我只想分析 Predictor 为‘year’的数据，该怎么做？” AI 会给出以下代码：

migrata%>%
 filter(Predictor == "year") -> migrationtime # 这将数据集简化为一个预测变量，时间

在开始建模之前，绘制数据是很有帮助的。漏斗图通常用于可视化元分析的数据，通过将预测变量与每个数据点的1/标准误差绘制在一起，根据精度对每个研究进行加权，标准误差高的研究权重较低。

并解释dplyr包中的filter()函数用于根据条件筛选数据，这里筛选出了 Predictor 变量值为“year”的行，并将结果存储在migrationtime变量中。

数据可视化（漏斗图绘制）：询问 AI 如何绘制漏斗图来可视化数据，AI会回复：

从漏斗图中可以看出，数据似乎在零附近聚集，并且正值和负值都有很好的表示，说明该研究没有受到发表偏倚的影响。进一步放大查看，可以更清楚地看到真实值似乎在零的左侧聚集，并且周围有相当大的变化。

点击标题查阅往期内容

R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）

左右滑动查看更多

解释plot()函数用于绘制图形，这里将Slope变量作为 x 轴，1/SE（精度）作为 y 轴绘制了漏斗图，以根据精度对每个研究进行加权。
如果你想进一步调整图形的范围，如 x 轴和 y 轴的范围，可以继续向 AI 提问：我希望你能扮演一名数据分析师。我有一个名为‘midata’的数据集，包含‘Slope’和‘SE’列。请绘制一个散点图，x 轴为‘Slope’，y 轴为‘1/SE’，并添加图表标题和轴标签，AI 会给出：

接下来，我们运行一个随机效应模型，仅将截距作为固定效应。截距将估计所有数据点的到达日期的平均变化，随机效应将估计截距周围是否存在真实变化，以及有多少变化可以由物种、位置或研究的效应来解释。

模型构建（随机效应模型）：向 AI 询问如何构建一个仅包含截距作为固定效应的随机效应模型，AI 会提供：

运行模型后，我们得到了一个估计参数的分布，因为运行了13000次模型迭代，并采样了1000次以提供后验分布。通过查看summary(randomtest)的摘要统计信息，我们可以看到每个效应的后验均值、分布的95%可信区间（不是置信区间）、有效样本大小以及固定效应的pMCMC值。
有效样本大小应该相当高（通常目标是1000-2000），更复杂的模型通常需要更多的迭代来达到可比的有效样本大小。

六、模型评估

（一）显著性评估

对于固定效应，当可信区间不跨越零时，我们可以认为该固定效应是显著的，因为如果后验分布跨越零，我们就不能确定它不是零。虽然会报告pMCMC值，但更应关注可信区间。理想情况下，后验分布应该很窄，表明该参数值被精确估计。
对于随机效应，我们估计方差。由于方差不能为零或负数，当方差的分布不接近零时，我们认为随机效应是显著的。可以通过绘制每个后验分布的直方图来检查这一点。

解释函数用于构建模型，Slope ~ 1表示只有截距作为固定效应，random = ~Species + Location + Study指定了随机效应为Species、Location和Study，data = migrationtime指定了使用的数据集。summary()函数用于查看模型的摘要统计信息。

模型评估 - 显著性评估：询问 AI 如何评估模型中固定效应和随机效应的显著性，AI 会回复：

从直方图中可以看出，位置和物种的方差分布接近零，对于随机效应要显著，我们希望其尾部远离零。

（二）模型收敛评估

检查模型收敛需要分别对固定效应和随机效应进行。对于固定效应，可以绘制randomtest$Sol来查看截距的轨迹和密度估计。轨迹类似于模型运行时的时间序列，可用于评估混合（或收敛）情况，而密度类似于模型每次迭代产生的后验分布估计的平滑直方图。

解释对于固定效应，当可信区间不跨越零时，认为该固定效应显著；对于随机效应，通过绘制方差的后验分布直方图来评估，当方差的分布不接近零时，认为随机效应显著。上述代码使用par()函数设置绘图布局，hist()函数绘制直方图，mcmc()函数用于处理模型的输出数据。

模型评估 - 收敛评估：AI 会给出：

为了确保模型已经收敛，轨迹图应该看起来像一只模糊的毛毛虫，从图中可以看出截距的混合情况良好。
如果怀疑存在过多的自相关，可以采取以下措施：

增加迭代次数，默认值为13000（例如nitt = 60000，对于更复杂的模型，可能需要使用几十万次迭代）。
增加燃烧期，默认情况下MCMCglmm会忽略前3000次迭代，因为此时模型尚未收敛，可以增加这个值（例如burnin = 5000）。
增加抽样间隔，默认值为10（例如thin = 30）。
考虑使用更强的先验。
对于随机效应的方差，同样绘制randomtest$VCV进行检查。根据电脑和屏幕的情况，可能会收到绘图太大无法显示的错误消息，可以通过向上和向左拖动绘图面板来扩大它，以便绘图有足够的空间显示。

从图中可以看出，一些随机效应的方差混合得不太好，有效样本大小也很小。可能需要增加迭代次数，但由于链似乎停留在零附近，看起来需要使用比默认值更强的先验。

七、先验知识

贝叶斯分析中最困难的部分是如何拟合正确的先验。先验是对我们认为参数的均值和/或方差可能是什么的先验知识的数学量化。我们为每个固定效应、随机效应和残差分别拟合一个先验。
先验可以用来告知模型我们认为后验分布将采取的形状。先验的信息性可以有所不同，弱信息先验用于我们没有太多先验知识并且希望数据自己说话的情况，它不会将后验分布从数据表明可能的参数值上拉开。信息性先验提供对模型估计至关重要的信息，并会在很大程度上影响后验分布的形状。
每个先验都遵循类似的公式，其是强信息还是弱信息取决于所包含的值。需要注意的是，没有完全无信息的先验。
默认先验假设固定效应的后验分布为正态分布，方差非常大，而对于随机效应的方差，实现了逆Wishart先验。逆Wishart先验包含方差矩阵V和置信度参数nu。
随着模型变得更加复杂，更有可能最终收到错误消息，或者模型从一开始就无法混合。在这种情况下，我们应该使用自己的参数扩展先验。参数扩展意味着先验不再是逆Wishart先验，而是缩放F先验（如果不理解也不用担心）。这不一定是坏事，因为参数扩展先验比逆Wishart先验更不容易指定错误。

八、参数扩展先验和测量误差

我们再次运行模型，但这次为随机效应使用参数扩展先验，通过包含prior = prior1。每个随机效应由一个G表示，残差由R表示。参数扩展是指我们包含了先验均值（alpha.mu）和（协）方差矩阵（alpha.V）以及V和nu。

先验知识与参数扩展先验：向 AI 询问，AI 会提供以下代码示例及解释：

这里增加了迭代次数到60000，以改善混合和有效样本大小。由于MCMC的随机性，每次重新运行模型时，输出都会略有不同，因此即使在模型中使用相同的效应，结果也会与这里打印的内容略有不同。
检查有效样本大小，发现现在有效样本大小大了很多，这是一个好迹象。再次绘制随机效应的方差图，发现模型的混合情况也更好了。
在进行模型检查之前，我们希望在模型中控制抽样误差，这是使用MCMCglmm进行元分析而不是其他程序或包的关键原因之一。通过拟合idh(SE):units作为随机效应并将相关方差固定为1，可以使用一个计算技巧来实现这一点。

控制测量误差：AI 会给出：

检查摘要统计信息，可以看到包含测量误差后，我们的估计更加保守，标准误差较高的研究被赋予了较低的统计权重。
为了进行模型检查，我们根据相同的参数值（方差/协方差结构（先验））模拟新数据，然后将其与实际数据进行绘图，以确保它们重叠。

模型检查：AI 会回复：

从图中可以看出，模拟数据与实际数据的拟合情况还算合理，尽管模拟数据可能稍微向左倾斜。
这里的问题比较复杂，但我们可以尝试分析一下。问题在于，低精度估计的抽样方差实际上高于SE2SE2（即元分析的假设不成立）。这意味着a) 仍然对低精度研究给予了过多的权重；b) 一些生物变异（单位方差）被高估了；c) 如果发表偏倚在低精度研究中更有可能发生，那么平均效应大小可能会有偏差。
一个快速的解决方案是查看观察间方差是否随着报告的标准误差的增加而比预期更快地增加。为此，我们可以估计方差，而不是假设它为1，并查看估计值是否大于1。
让我们重新运行模型，但这次更改测量误差的先验，使其不再固定为1。

现在我们可以再次模拟新数据，并将其与我们收集的数据进行绘图。

这些参数似乎更适合我们的数据。

九、其他内容

其他分析（固定效应、计算后验均值、非高斯族、协方差结构等）：根据你的具体需求，向 AI 提问关于其他分析的问题，例如：

“如何在模型中添加固定效应？”
“如何计算随机效应的后验均值？”
“如何处理非高斯族数据？”
“如何构建协方差结构？”
AI 会根据你的问题提供相应的代码和解释，你只需按照提示在 R 语言中运行代码并理解其含义即可。

（一）固定效应

除了随机效应，还可以拟合固定效应。MCMCglmm估计随机效应的方式与固定效应类似，但在随机效应元分析中，分析师通常关注的是方差。

需要注意的是，对于固定效应，可能需要根据具体项目研究是否需要更改先验。

（二）计算随机效应的后验均值

MCMCglmm估计随机效应的方差和每个类别内的真实效应大小，但报告随机效应的方差比报告每个效应大小更有信息性。当使用summary()时，R会报告随机效应的方差和可信区间，但不会报告效应大小。然而，可以保存这些效应大小的后验模式并在工作中报告，但绝不能对其进行进一步的统计分析，并始终确保让读者知道预测的来源。

（三）非高斯族

本教程基于使用高斯分布，但MCMCglmm也可以处理非高斯族。通过指定family=来选择正确的分布。

（四）计算95%可信区间

可以使用interval(mcmc())来绘制或报告模型的可信区间，而无需直接从摘要中提取数字。

上述代码的结果应该与查看摘要时截距的后验分布的95%可信区间的数值相似。当想要组合效应时，interval特别有用。例如，若想知道来自欧洲的短距离迁徙者的后验分布的均值以及上下95%可信区间，可以像下面这样操作：

这些数值随后可用于绘图、报告等用途。

（五）（协）方差结构

到目前为止，我们了解到对于模型中的每个随机效应和残差，MCMCglmm会估计该效应内的方差，即方差在一个1x1的矩阵 - [_V_]中。然而，如果我们有需求，可以重新构建随机效应和残差内的方差矩阵。
举个例子，在之前的模型中，我们假设所有到达日期的测量方式都是相同的，但实际上，它有三种不同的测量方式：首先是到达的平均日期和中位数日期，峰值到达也作为一个类别包含在内，不过没有包含该类别的行。

因此，我们的残差方差是异质的，需要在模型中考虑到这一点。我们可以通过在rcov中使用idh():units函数来实现。由于我们希望分别估计每个水平的方差，所以必须更改残差先验的方差结构。在这种情况下，我们使用一个3x3的方差矩阵，因为有三种类型的响应。

如果只是在R控制台中运行prior4，应该能够更轻松地可视化残差先验的矩阵。

现在，当我们打印摘要时可以看到，已经为所有三种到达测量方式估计了残差方差，这是一个成功的结果。
最后，除了使用方差矩阵，还可以使用协方差矩阵，将idh()替换为us()。在这种情况下，需要更新效应的先验。例如，如果有三个水平，则需要将矩阵的大小增加到3。