2

R で coda (mcmc.list) オブジェクトを保存する方法がわかりません。他の人から同様の質問がありましたが、回答が特に明確ではないことがわかりました。理想的には、coda オブジェクトを R.data ファイルまたはテキスト ファイル (csv など) として保存し、モデルを再実行することなく再インポートして JAGS チェーンを分析できるようにしたいと考えています (これには約 30 分かかります)。私のコンピューター）。現在、私の coda オブジェクト "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 の 3 つのチェーン)。この coda オブジェクトを保存して、毎回モデルを再実行することなく R にコールバックできるようにしたいと考えています。また、独自のチェーンが 3 つあるという事実も維持したいと考えています。

4

1 に答える 1

1

これは、 からチェーンを保存するときに使用するスクリプトですMCMCglmm。OBJECTを coda オブジェクトの名前 (例: で作成したモデルの名前MCMCglmm) に置き換え、FILEPATH適切な保存先に置き換えます。

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

if便利なヒント: また、 andを使用してモデルをスイッチでラップします。elseこれを組み合わせて、スクリプトを実行して保存するか、以前の実行をロードする (またはに設定runmodし、適切なファイルを取得することを選択する) 便利なシステムを作成できます。たとえば、 :ynloaddate

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 に答える