rmcmc

Gelman-Rubin statistic for MCMCglmm model in R


I have a multivariate model with this (approximate) form:

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)

For publication purposes, I've been asked to report the Gelman-Rubin statistic to indicate that the model has converged.

I have been trying to run:

gelman.diag(mod.1) 

But, I get this error:

Error in mcmc.list(x) : Arguments must be mcmc objects

Any suggestions on the proper approach? I assume that the error means I can't pass my mod.1 output through gelman.diag(), but I am not sure what it is I am supposed to put there instead? My knowledge is quite limited here, so I'd appreciate any and all help!

Note that I haven't added the data here, but I suspect the answer is more code syntax and not data related.


Solution

  • The gelman.diag requires a mcmc.list. If we are running models with different set of parameters, extract the 'Sol' and place it in a list (Below, it is the same model)

    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
    

    According to the documentation, it seems to be applicable for more than one mcmc chain

    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.

    The input x here is

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