2

我不确定如何在 R 中保存 coda (mcmc.list) 对象。其他人也问过类似的问题，但我发现给出的答案并不是特别清楚。理想情况下，我想将 coda 对象保存为 R.data 文件或文本文件（例如 csv），这样我就可以重新导入它并分析 JAGS 链，而无需重新运行模型（大约需要 30 分钟）我的电脑）。现在我的尾声对象“coda.samples”看起来像这样：

str(coda.samples)
List of 3
 $ : mcmc [1:3334, 1:1094] 0.904 0.977 0.927 0.945 0.905 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3334] "2001" "2011" "2021" "2031" ...
  .. ..$ : chr [1:1094] "alpha[1,1]" "alpha[2,1]" "alpha[1,2]" "alpha[2,2]" ...
  ..- attr(*, "mcpar")= num [1:3] 2001 35331 10
 $ : mcmc [1:3334, 1:1094] 0.824 0.866 0.839 0.832 1.032 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3334] "2001" "2011" "2021" "2031" ...
  .. ..$ : chr [1:1094] "alpha[1,1]" "alpha[2,1]" "alpha[1,2]" "alpha[2,2]" ...
  ..- attr(*, "mcpar")= num [1:3] 2001 35331 10
 $ : mcmc [1:3334, 1:1094] 0.956 0.944 0.895 0.809 1.064 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3334] "2001" "2011" "2021" "2031" ...
  .. ..$ : chr [1:1094] "alpha[1,1]" "alpha[2,1]" "alpha[1,2]" "alpha[2,2]" ...
  ..- attr(*, "mcpar")= num [1:3] 2001 35331 10
 - attr(*, "class")= chr "mcmc.list"

如您所见，它是三个矩阵的列表，其中包含 1094 个参数的 3334 个估计值（即 3 个长度为 3334 的链）。我想存储这个 coda 对象，这样我就可以将它调用回 R 中，而不必每次都重新运行模型。我还想保留三个独特链的事实。

4

1 回答 1

1

这是我在保存链时使用的脚本MCMCglmm。将 替换OBJECT为您的 coda 对象的名称（例如，在 中制作的模型名称MCMCglmm），并替换FILEPATH为合适的保存目标。

save(OBJECT, file= "FILEPATH")
model = load("FILEPATH")

实用提示：我还使用开关将我的模型if和else- 这可以组合成一个有用的系统来运行和保存脚本，或加载以前的运行（设置runmod为yor n，并选择loaddate获取正确的文件。例如：

runmod = "y"
loaddate = "2017-01-12"

NITT = 130000
BURN =  30000
THIN =    100

# Model
if(runmod=="n"){
model4.1 = MCMCglmm(cbind(BWT_F, BWT_M, TAR_F, TAR_M) ~ trait-1,
        random = ~us(trait):animal + us(trait):BYEAR + us(trait):MOTHER,
        rcov = ~us(trait):units,
        family = rep("gaussian", times = 4),
        pedigree = Ped,
        data = Data,
        burnin = BURN,
        nitt = NITT,
        thin = THIN,
        prior = prior4.1)
        save(model4.1, file=paste0("FILEPATH......",NITT,"_",Sys.Date(),".rda"))
}else{                       
model41 = load(paste0("FILEPATH......",NITT,"_",loaddate,".rda"))
}
于 2017-01-12T11:08:27.350 回答