mcmcglmm循環創建多條鏈路

Question

從this question隨訪（見重複性的數據幀）我想運行MCMCGLMM ñ次，其中ñ是randomisations的數量。我試圖構建一個運行所有鏈的循環，並保存它們（以後檢索隨機變量的後驗分佈），但我遇到了問題。 （n = 5，因此R1-R5），A =響應變量，L和V是隨機效應變量，B是固定效應，R1-R5是隨機的用作爲V的結構L的分配維持：

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-R ñ）：

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

• 感謝羅蘭幫助，以獲得隨機效應正確地調用，以前我得到了一個錯誤Error in buildZ(rmodel.terms[r] ... object Rs[i] not found - 由固定as.formula

但是，這存儲所有t他在chains數據和表面上僅$Sol組件，但我需要能夠在VCV內訪問值，特別是ř變量的後驗分佈（例如summary(chainR1$VCV)）

總結：看來我在做我該怎麼分配鏈的名字搞錯，有沒有人對如何做到這一點的建議，並保存後驗分佈，甚至整個產業鏈？

Answer 1

使用分配是關鍵點：

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

然後它成爲可能建立我對（在列即隨機L處的賦值R1-R ñ）的VCV成分的基質

Answer 2

它好像要在一個循環中運行多個不同MCMCglmm公式。 @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是一樣的貼，但默認參數sep=""