[精彩WinBUGS答问]How to convert winbugs code to R

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Since I am new to winbugs, I need to conduct a bayesian analysis in R. To do so, I need to convert the winbugs code below into R code:

model{
# model’s likelihood
for (i in 1:n){
time ~ dnorm( mu, tau ) # stochastic componenent
# link and linear predictor
mu <- beta0 + beta1 * cases + beta2 * distance
}
# prior distributions
tau ~ dgamma( 0.01, 0.01 )
beta0 ~ dnorm( 0.0, 1.0E-4)
beta1 ~ dnorm( 0.0, 1.0E-4)
beta2 ~ dnorm( 0.0, 1.0E-4)
# definition of sigma
s2<-1/tau
s <-sqrt(s2)
# calculation of the sample variance
for (i in 1:n){ c.time<-time-mean(time[]) }
sy2 <- inprod( c.time[], c.time[] )/(n-1)
# calculation of Bayesian version R squared
R2B <- 1 - s2/sy2
# Expected y for a typical delivery time
typical.y <- beta0 + beta1 * mean(cases[]) + beta2 * mean(distance[])
}
INITS
list( tau=1, beta0=1, beta1=0, beta2=0 )
DATA (LIST)
list( n=25,
time = c(16.68, 11.5, 12.03, 14.88, 13.75, 18.11, 8, 17.83,
79.24, 21.5, 40.33, 21, 13.5, 19.75, 24, 29, 15.35,
19, 9.5, 35.1, 17.9, 52.32, 18.75, 19.83, 10.75),
distance = c(560, 220, 340, 80, 150, 330, 110, 210, 1460,
605, 688, 215, 255, 462, 448, 776, 200, 132,
36, 770, 140, 810, 450, 635, 150),
cases = c( 7, 3, 3, 4, 6, 7, 2, 7, 30, 5, 16, 10, 4, 6, 9,
10, 6, 7, 3, 17, 10, 26, 9, 8, 4) )
and this is the R code I used

library(R2WinBUGS)

time <-      c(16.68,11.5,12.03,14.88,13.75,18.11,8,17.83,79.24,21.5,40.33,21,13.5,19.75,24,29,15.35,19,9.5,35.1,17.9,52.32,18.75,19.83,10.75)
cases <- c(7,3,3,4,6,7,2,7,30,5,16,10,4,6,9,10,6,7,3,17,10,26,9,8,4)
distance <-   c(560,220,340,80,150,330,110,210,1260,605,688,215,255,462,448,776,200,132,36,770,140,810,450,635,150)
n <- length(time)

data <- list("time","cases","distance")
inits <- function(){
list(tau=1,n=25,beta1=rep(0,25),beta2=rep(0.25),beta3=rep(0,25))
}

sim <- bugs(data, inits, model.file = "C:/Users/Gunal/Desktop/dummy/reg.txt",
parameters = c("beta1", "beta2","beta3"),
n.chains = 3, n.iter = 1000,codaPkg = FALSE,
bugs.directory = "D:/PROGRAMLAR/WinBUGS14/",debug=TRUE)

print(sim)
and here is the "reg.txt" file aused above

model{
# model’s likelihood
for (i in 1:n){
time ~ dnorm( mu, tau ) # stochastic componenent
# link and linear predictor
mu <- beta0 + beta1 * cases + beta2 * distance
}
# prior distributions
tau ~ dgamma( 0.01, 0.01 )
beta0 ~ dnorm( 0.0, 1.0E-4)
beta1 ~ dnorm( 0.0, 1.0E-4)
beta2 ~ dnorm( 0.0, 1.0E-4)
# definition of sigma
s2<-1/tau
s <-sqrt(s2)
# calculation of the sample variance
for (i in 1:n){ c.time<-time-mean(time[]) }
sy2 <- inprod( c.time[], c.time[] )/(n-1)
# calculation of Bayesian version R squared
R2B <- 1 - s2/sy2
# Expected y for a typical delivery time
typical.y <- beta0 + beta1 * mean(cases[]) + beta2 * mean(distance[])
}
finally this is the error launched by winbugs:

display(log)
check(C:/Users/Gunal/Desktop/dummy/reg.txt)
model is syntactically correct
data(C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/data.txt)
data loaded
compile(3)
variable n is not defined
inits(1,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/inits1.txt)
command #Bugs:inits cannot be executed (is greyed out)
inits(2,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/inits2.txt)
command #Bugs:inits cannot be executed (is greyed out)
inits(3,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/inits3.txt)
command #Bugs:inits cannot be executed (is greyed out)
gen.inits()
command #Bugs:gen.inits cannot be executed (is greyed out)
thin.updater(1)
update(500)
command #Bugs:update cannot be executed (is greyed out)
set(beta1)
command #Bugs:set cannot be executed (is greyed out)
set(beta2)
command #Bugs:set cannot be executed (is greyed out)
set(beta3)
command #Bugs:set cannot be executed (is greyed out)
set(deviance)
command #Bugs:set cannot be executed (is greyed out)
dic.set()
command #Bugs:dic.set cannot be executed (is greyed out)
update(500)
command #Bugs:update cannot be executed (is greyed out)
coda(*,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/coda)
command #Bugs:coda cannot be executed (is greyed out)
stats(*)
command #Bugs:stats cannot be executed (is greyed out)
dic.stats()

DIC
history(*,C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/history.odc)
command #Bugs:history cannot be executed (is greyed out)
save(C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/log.odc)
save(C:/Users/Gunal/AppData/Local/Temp/RtmpUbKAAJ/log.txt)
Obviously, I failed to define variable n. Does anyone know what is wrong and how to fix it. Any help would be greatly appreciated.

Many thanks in advance

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：winbugs Convert WINBUG BUGS code convert conduct

Something like...

time <-c(16.68,11.5,12.03,14.88,13.75,18.11,8,17.83,79.24,21.5,40.33,21,13.5,19.75,24,29,15.35,19,9.5,35.1,17.9,52.32,18.75,19.83,10.75)
cases <- c(7,3,3,4,6,7,2,7,30,5,16,10,4,6,9,10,6,7,3,17,10,26,9,8,4)
distance <-c(560,220,340,80,150,330,110,210,1260,605,688,215,255,462,448,776,200,132,36,770,140,810,450,635,150)
n <- length(time)
data <- list(n=n,time=time,cases=cases,distance=distance)

sim <- bugs(data, inits=NULL, model.file = "C:/Users/Gunal/Desktop/dummy/reg.txt",
        parameters = c("beta0","beta1", "beta2"),
        n.chains = 3, n.iter = 1000, codaPkg = FALSE,
        bugs.directory = "D:/PROGRAMLAR/WinBUGS14/")
should work. A problem was in the beta's. Your beta's in the BUGS code run from beta0 to beta2. Your beta's in your R code run from beta1 to beta3.

If you want to set initial values too (the initial values in your R code for beta's are vectors of length 25, when they should only be a single number, as in your BUGS file) then this should work...

inits <- function(){
  list(tau=1,beta0=0,beta1=0,beta2=0)
}
sim <- bugs(data, inits=inits, model.file = "C:/Users/Gunal/Desktop/dummy/reg.txt",
        parameters = c("beta0","beta1", "beta2"),
        n.chains = 3, n.iter = 1000, codaPkg = FALSE,
        bugs.directory = "D:/PROGRAMLAR/WinBUGS14/")
some other initial values (for time[i] nodes) will still be generated by WinBUGS.

并没有将winbugs转成R，只是用R调用winbugs。