mcmcglmm 循环创建许多链

Question

从this question 跟进（请参阅可重现的数据框）我想运行 MCMCGLMM n 次，其中 n 是随机化的次数。我试图构建一个运行所有链并保存它们的循环（以便稍后检索随机变量的后验分布），但我遇到了问题。

这是数据框的样子（当n = 5，因此R1-R5），A =响应变量，L和V是随机效应变量，B是固定效应，R1 -R5 是 L 的随机分配，保持 V 的结构：

  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

我可以创建我想分配给我的链的名称，以及随着 MCMC 链的每次运行而变化的变量的名称 (R1-Rn)：

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

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

输出：

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

然后我尝试了这个循环来产生 5 个链：

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)
}

• 感谢 Roland 帮助获得正确调用的随机效果，之前我收到错误 Error in buildZ(rmodel.terms[r] ... object Rs[i] not found- 由 as.formula 修复

但这将所有数据存储在chains 中，并且似乎只有$Sol 组件，但我需要能够访问 VCV 中的值，特别是 R的后验分布> 变量（例如summary(chainR1$VCV)）

总结：似乎我在分配链名称时犯了一个错误，有没有人建议如何做到这一点，并保存后验分布甚至整个链？

Answer 1

使用 assign 是一个关键点：

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]
                }

然后可以构建我感兴趣的 VCV 组件的矩阵（即 R1-Rn 列中的随机 L 分配）

Answer 2

您似乎想在一个循环中运行许多不同的 MCMCglmm 公式。 @Roland 帮助您找到了解决方案（尽管我个人会在循环之前创建公式）。 @Roland 还指出，为了保存每个模型的结果，您应该将它们保存在一个列表中 - 而不是像您目前所做的那样保存在一个链中。您还可以将每个模型保存为 .RData 文件，如问题末尾所示。为了正式回答这个问题，我将通过以下方式执行此操作：

Rs = paste0("~R", 1:5, " + V") ## Create all model formulae
chainNames = paste0("chainR", 1:5) ## Names for each model 
chains = list() ## Initialize list
## Loop over models
for(i in 1:length(Rs)){
chains[[i]] =   MCMCglmm(A ~1 + B,
                random = formula(Rs[i]),
                rcov = ~units,
                nitt = 500,
                thin = 2,
                burnin = 50,
                prior = prior2,
                family = "gaussian",
                start = list(QUASI = FALSE),
                data = df)
}
names(chains) = chainNames ## Name each model
save(chains, "chainsR1-R5.Rdata") ## Save all model output 

附注，paste0 与 paste 相同，但默认使用参数sep=""