R 中 MCMCglmm 模型的 Gelman-Rubin 统计量

问题描述

我有一个具有这种（近似）形式的多元模型：

library(MCMCglmm)

mod.1 <- MCMCglmm(
    cbind(OFT1,MIS1,PC1,PC2) ~ 
    trait-1 +
    trait:sex +
    trait:date,random = ~us(trait):squirrel_id + us(trait):year,rcov = ~us(trait):units,family = c("gaussian","gaussian","gaussian"),data= final_MCMC,prior = prior.invgamma,verbose = FALSE,pr=TRUE,#this saves the BLUPs 
    nitt=103000,#number of iterations
    thin=100,#interval at which the Markov chain is stored
    burnin=3000)

出于发表目的，我被要求报告 Gelman-Rubin 统计数据以表明模型已经收敛。

我一直在尝试跑步：

gelman.diag(mod.1)

但是，我收到此错误：

Error in mcmc.list(x) : Arguments must be mcmc objects

对正确方法有什么建议吗？我认为错误意味着我无法通过 mod.1 传递我的 gelman.diag() 输出，但我不确定我应该把它放在那里吗？我的知识在这里非常有限，所以我很感激任何帮助！

请注意，我没有在此处添加数据，但我怀疑答案更多是代码语法而不是数据相关。

解决方法

gelman.diag 需要 mcmc.list。如果我们正在运行具有不同参数集的模型，则提取“Sol”并将其放置在 list 中（下面，它是相同的模型）

library(MCMCglmm)
model1 <- MCMCglmm(PO~1,random=~FSfamily,data=PlodiaPO,verbose=FALSE,nitt=1300,burnin=300,thin=1)
model2 <- MCMCglmm(PO~1,thin=1 )

mclist <- mcmc.list(model1$Sol,model2$Sol)
gelman.diag(mclist)
# gelman.diag(mclist)
#Potential scale reduction factors:
#            Point est. Upper C.I.
#(Intercept)          1          1

根据文档，它似乎适用于一个以上的 mcmc 链

Gelman 和 Rubin (1992) 提出了一种监测 MCMC 输出收敛的通用方法，其中 m > 1 个平行链以相对于后验分布过度分散的起始值运行。

这里的输入 x 是

x - 一个 mcmc.list 对象，具有多个链，并且具有相对于后验分布过度分散的起始值。

mcmc r r