我有一个与 R 代码相关的问题,该代码称为 BUGS。我已经在 WinBUGS 中运行了该模型,它运行良好,给出了预期的结果。下面是当我有 Y 的单一结果或单变量数据时使用的自动化代码。现在我想用它来实现多种结果。我尝试了一种不同的方式来读取数据。有 2 个用于测试的模拟是从 csv 文件中读取的。不确定在代码中的何处指定,以便可以对 2 个结果而不是一个结果重复相同的过程。 setwd(“C://Tina/USB_Backup_042213/Testing/CSV”)
matrix=NULL
csvs <- paste("MVN", 1:2, ".csv", sep="")
for(i in 1:length(csvs)){
matrix[[i]] <- read.csv(file=csvs[i], header=T)
print(matrix[[i]])
}
Y1 Y2
1 11 6
2 8 5
3 25 13
4 1 13
5 8 22
Y1 Y2
1 9 1
2 7 9
3 25 13
4 1 18
5 9 12
library("R2WinBUGS")
bugs.output <- list()
for(sim in 1:2){
Y <-(matrix[[sim]])
bugs.output[sim] <- bugs(
data=list(Y=as.matrix(Y), Nf=5, n=60, mn=c(-1.59, -2.44), prec=matrix(c(.0001,0,0,.0001),nrow=2,ncol=2), R=matrix(c(.001,0,0,.001),nrow=2,ncol=2)),
inits=list(list(gamma=c(0,0), T=matrix(c(0.9,0,0,0.9),nrow=2,ncol=2))),
model.file="M-LN_model_trial.txt",
parameters.to.save = c("p","rho","sigma2"),
n.chains=1, n.iter=12000, n.burnin=5000, debug=TRUE,
bugs.directory="C://Tina/USB_Backup_042213/winbugs14/WinBUGS14",
working.directory=NULL)
}
警告消息: 1:在 bugs.output[sim] <- bugs(data = list(Y = as.matrix(Y), Nf = 5, : 要替换的项目数量不是替换长度的倍数 2:在 bugs.output[sim] <- bugs(data = list(Y = as.matrix(Y), Nf = 5, : 要替换的项目数量不是替换长度的倍数
最佳答案
当您从 R 运行 BUGS 模型时遇到错误,一种选择是尝试在 OpenBUGS 或 WinBUGS 本身中模拟运行该模型。它可以帮助您(通过点击检查模型按钮后的光标位置)找到有问题的行。
我用你的 BUGS 模型做了这个。我发现BUGS模型中mn
、prec
和R
的定义有问题。您可以删除它们,因为它们已经在数据中定义了(看起来是定义它们的适当位置)。一旦我从你的 BUGS 模型中删除这些,一切都运行良好。
注意,要在 OpenBUGS 中运行模型,您必须编辑数据格式,例如我运行的脚本是:
model{
#likelihood
for(j in 1 : Nf){
p1[j, 1:2 ] ~ dmnorm(gamma[1:2], T[1:2 ,1:2])
for (i in 1:2){
logit(p[j,i]) <- p1[j,i]
Y[j,i] ~ dbin(p[j,i],n)
}
}
#priors
gamma[1:2] ~ dmnorm(mn[1:2],prec[1:2 ,1:2])
expit[1] <- exp(gamma[1])/(1+exp(gamma[1]))
expit[2] <- exp(gamma[2])/(1+exp(gamma[2]))
T[1:2 ,1:2] ~ dwish(R[1:2 ,1:2], 2)
sigma2[1:2, 1:2] <- inverse(T[,])
rho <- sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2])
}
#data
list(Y=structure(.Data=c(1,11,6,1,8,5,1,25,13,1,1,13,1,8,22),.Dim=c(5,3)),
Nf=5, n=60, mn=c(-1.59,-2.44),
prec=structure(.Data=c(0.0001,0,0,0.0001),.Dim=c(2,2)),
R=structure(.Data=c(0.001,0,0,0.001),.Dim=c(2,2)))
#inits
list(gamma=c(0,0), T=structure(.Data=c(0.9,0,0,0.9),.Dim=c(2,2)))
其中数据和初始化需要一些工作才能从 R 脚本转换。
其他几点:1)我不确定 Y 的格式是否正确,因为它有 3 列,您的分布仅考虑前两列(X 和 Y1)。 2)你可能有一组不必要的大括号。
要通过 R 运行 BUGS 中的代码,您可以使用以下 R 语法...
#BUGS code as a character string
bugs1<-
"model{
#likelihood
for(j in 1 : Nf){
p1[j, 1:2 ] ~ dmnorm(gamma[1:2], T[1:2 ,1:2])
for (i in 1:2){
logit(p[j,i]) <- p1[j,i]
Y[j,i] ~ dbin(p[j,i],n)
}
}
#priors
gamma[1:2] ~ dmnorm(mn[1:2],prec[1:2 ,1:2])
expit[1] <- exp(gamma[1])/(1+exp(gamma[1]))
expit[2] <- exp(gamma[2])/(1+exp(gamma[2]))
T[1:2 ,1:2] ~ dwish(R[1:2 ,1:2], 2)
sigma2[1:2, 1:2] <- inverse(T[,])
rho <- sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2])
}"
#write the BUGS code to a txt file in current working directory
writeLines(bugs1, "bugs1.txt")
#create data
Y<-data.frame(X=1,Y1=c(11,8,25,1,8),Y2=c(6,5,13,13,22))
#run BUGS from R
library("R2OpenBUGS")
mcmc1 <- bugs(data = list(Y=as.matrix(Y), Nf=5, n=60, mn=c(-1.59, -2.44),
prec=matrix(c(.0001,0,0,.0001),nrow=2,ncol=2),
R=matrix(c(.001,0,0,.001),nrow=2,ncol=2)),
inits = list(list(gamma=c(0,0), T=matrix(c(0.9,0,0,0.9),nrow=2,ncol=2))),
param = c("gamma", "sigma2"),
model = "bugs1.txt",
n.iter = 11000, n.burnin = 1000, n.chains = 1)
这里有几点需要注意。 1) 这使用 OpenBUGS 而不是 WinBUGS。 2) 如果您使用 R2WinBUGS,如果您没有以管理员身份运行 R(或 Rstudio,或您正在使用的任何内容),则可能会遇到陷阱。
要运行上面的代码 1000 次,您可以将其放入循环中,例如......
#create and write the BUGS code to a txt file in current working directory (outside the loop)
bugs1<-...
#loop
for(i in 1:1000){
Y <- read.csv(file=paste0("MVN",i,".csv"))
#run BUGS from R
library("R2OpenBUGS")
mcmc1 <- bugs(data = list(Y=as.matrix(Y), Nf=5, n=60, mn=c(-1.59, -2.44),
prec=matrix(c(.0001,0,0,.0001),nrow=2,ncol=2),
R=matrix(c(.001,0,0,.001),nrow=2,ncol=2)),
inits = list(list(gamma=c(0,0), T=matrix(c(0.9,0,0,0.9),nrow=2,ncol=2))),
param = c("gamma", "sigma2"),
model = "bugs1.txt",
n.iter = 11000, n.burnin = 1000, n.chains = 1)
#save mcmc
write.csv(mcmc1$sims.matrix,paste0("mcmc",i,".csv"))
}
关于r - R 中的 WinBUGS 自动化,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/19068493/