时间戳

首页 > 解决方案 > R中MCMCglmm模型的Gelman-Rubin统计

r - 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", "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()，但我不确定我应该把它放在那里吗？我的知识在这里非常有限，所以我将不胜感激任何和所有的帮助！

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

标签： rmcmc

解决方案

gelman.diag需要mcmc.list一个. 如果我们正在运行具有不同参数集的模型，请提取“Sol”并将其放在 a 中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 和 Rubin (1992) 提出了一种监测 MCMC 输出收敛性的通用方法，其中 m > 1 个平行链以相对于后验分布过度分散的起始值运行。

这里的输入x是

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

推荐阅读