Когда вы получаете ошибку при запуске вашей модели BUGS из R, одним из вариантов является попытка макетного запуска модели в OpenBUGS или WinBUGS сам. Он может помочь вам (через размещение курсора после нажатия кнопки модели проверки), чтобы найти проблемные линии.
Я сделал это с вашей моделью BUGS. Я нашел проблемы в определении mn
, prec
и R
в модели BUGS. Вы можете отказаться от них, поскольку они уже определены в данных (что, похоже, подходит для их определения). Как только я сбросил их с вашей модели 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)))
, где данные и inits нужно немного работы, чтобы преобразовать из вашего R.
Несколько других моментов: 1) Я не уверен, что у вас есть правильный формат для Y, так как он имеет 3 столбца, в вашем распределении учитываются только первые два (X и Y1). 2) Вероятно, у вас был ненужный набор фигурных скобок.
Для запуска кода в BUGS с помощью R вы можете использовать следующий синтаксис 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"))
}
это официально не запрещено, но, пожалуйста, пересмотреть кросс-постинг на SO и г-помощь HTTP : //article.gmane.org/gmane.comp.lang.r.general/300245 ... он тратит силы. –
Можем ли мы увидеть модель BUGS? Я предполагаю, что его нужно будет настроить для обработки многомерных данных (и это наиболее вероятно, где ваша проблема). – gjabel
@ Бен Болкер, передумает в следующий раз. Я не был уверен, что в listservs и stackexchange есть одна и та же аудитория. – user1560215