Search code examples
rloopsmcmc

mcmcglmm loop to create many chains


Following up from this question (see for reproducible data frame) I want to run MCMCGLMM n times, where n is the number of randomisations. I have tried to construct a loop which runs all the chains, and saves them (to retrieve the posterior distributions of the randomised variable later) but I am encountering problems.

This is what the data frame looks like (when n = 5, hence R1-R5), A = response variable, L and V are random effect variables, B is a fixed effect, R1-R5 are random assignments of L with structure of V maintained:

  ID    L B V    A R1 R2 R3 R4 R5
1 1_1_1 1 1 1 11.1  6 19 21  1 31
2 1_1_1 1 1 1  6.9  6 19 21  1 31
3 1_1_4 1 1 4  7.7  2 24  8 22 22
4 1_1_4 1 1 4 10.5  2 24  8 22 22
5 1_1_5 1 1 5  8.5 11 27 14 17 22
6 1_1_7 1 1 7 11.2  5 24 13 18 25

I can create the names I want to assign to my chains, and the names of the variable that changes with each run of the MCMC chain (R1-Rn):

n = 5
Rs = as.vector(rep(NA,n))

for(i in 1:n){
 Rs[i] =     paste("R",i, sep = "")
 }
Rs

Output:

> Rs
[1] "R1" "R2" "R3" "R4" "R5"

I then tried this loop to produce 5 chains:

for(i in 1:n){
chains[i] =     MCMCglmm(A ~1 + B,
                random = as.formula(paste0("~" ,Rs[i], " + Vial")),
                rcov = ~units,
                nitt = 500,
                thin = 2,
                burnin = 50,
                prior = prior2,
                family = "gaussian",
                start = list(QUASI = FALSE),
                data = df)
}
  • Thanks Roland for helping to get the random effect to call properly, previously I was getting an error Error in buildZ(rmodel.terms[r] ... object Rs[i] not found- fixed by as.formula

But this stores all of the data in chains and seemingly only the $Sol components, but I need to be able to access the values within the VCV, specifically the posterior distributions of the R variables (e.g. summary(chainR1$VCV))

In summary: It seems I am making a mistake in how I assign the chain names, does anyone have a suggestion of how to do this, and save the posterior distributions or even the whole chain?


Solution

  • Using assign was a key point:

    n = 10 #Number of chains to run
    chainVCVdf = matrix(rep(NA, times = ((nitt-burnin)/thin)*n), ncol = n)
    colnames(chainVCVdf)=c(rep("X", times = n))
    
    for(i in 1:n){
    assign("chainX",paste0("chain",Rs[i]))
    chainX =    MCMCglmm(A ~1 + B,
                    random = as.formula(paste0("~" ,Rs[i], " + V")),
                    rcov = ~units,
                    nitt = nitt,
                    thin = thin,
                    burnin = burnin,
                    prior = prior1,
                    family = "gaussian",
                    start = list(QUASI = FALSE),
                    data = df)
    assign("chainVCV",  chainX$VCV[,1]) 
    chainVCVdf[,i]=(chainVCV)   
    colnames(chainVCVdf)[i] = colnames(chainX$VCV)[1]
                    }
    

    It then became possible to build a matrix of the VCV component that I am interested in (namely the randomised L assignment in columns R1-Rn)