2

I'm unsure of how to save a coda (mcmc.list) object in R. Others have asked similar questions, but I found that the answers given were not particularly clear. Ideally I'd like to save the coda object as an R.data file or a text file (e.g., csv), so that I can reimport it and analyze the JAGS chains without having to rerun the model (which takes about 30 minutes on my computer). Right now my coda object "coda.samples" looks like this:

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"

As you can see, it's a list of three matrices that contain 3334 estimates of 1094 parameters (i.e., 3 chains of length 3334). I'd like to store this coda object so that I can call it back into R without having to rerun the model over again each time. I also want to preserve the fact that there are three unique chains.

David Johnson
  • 439
  • 1
  • 3
  • 13
  • 2
    Take a look at `?save` or `?saveRDS`. Any R object can be saved with, e.g., `save(coda.samples, file='my.coda.samples.rda')`. This can then be loaded back into R with `load('my.coda.samples.rda')`. – jbaums Feb 17 '15 at 04:27
  • 1
    You're absolutely right. I was overthinking this. Using the command `save(coda.samples, file = "coda.samples.RData")` saves the coda object. Then loading it with `load(coda.samples.RData)` works perfectly. – David Johnson Feb 17 '15 at 04:42

1 Answers1

1

This is the script I use when saving chains from MCMCglmm. Replace the OBJECT with the name of your coda object (e.g. the name of the model as made in MCMCglmm), and replace FILEPATH with a suitable save destination.

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

Handy Tip: I also wrap my models in switches using if and else - this can be combined to make a useful system that either runs and saves the script, or loads a previous run (setting runmod to y or n, and choosing loaddate to get the right file. For example:

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"))
}
rg255
  • 4,119
  • 3
  • 22
  • 40