R 中 MCMCglmm 模型的 Gelman-Rubin 统计

Question

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

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", "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

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

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

Answer 1

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, random=~FSfamily, data=PlodiaPO, verbose=FALSE,
  nitt=1300, burnin=300, 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 and Rubin (1992) propose a general approach to monitoring convergence of MCMC output in which m > 1 parallel chains are run with starting values that are overdispersed relative to the posterior distribution.

这里的输入x是

x - An mcmc.list object with more than one chain, and with starting values that are overdispersed with respect to the posterior distribution.

R 中 MCMCglmm 模型的 Gelman-Rubin 统计

Gelman-Rubin statistic for MCMCglmm model in R

r

mcmc