mcmcglmm 循环创建许多链

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[=48 修复=]

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

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

Answer 1

似乎您想运行循环中的许多不同的 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=""

Answer 2

使用 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-R 列中的随机 L 分配n）

mcmcglmm 循环创建许多链

mcmcglmm loop to create many chains

loops

r

mcmc