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。