# Extracting values from R output of bugs command [closed]

I asked this question on https://stackoverflow.com/, but I couldn't get what I want. So, I am asking it here.

I am running winbugs from R and I need to use some variables in R output. When I type schools.sim\$mean\$theta[1] in R, I get 10.2. However, when I type schools.sim\$2.5%\$theta[1]an error message come up. Any one what I am doing wrong or any other way to get the bayesian intervals?

here is an example:

This is the R code

library(R2WinBUGS)
data(schools)
J <- nrow(schools)
y <- schools$estimate sigma.y <- schools$sd
data <- list ("J", "y", "sigma.y")

inits <- function(){
list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),
sigma.theta = runif(1, 0, 100))
}

schools.sim <- bugs(data, inits, model.file = "D:/model.txt",
parameters = c("theta", "mu.theta", "sigma.theta"),
n.chains = 3, n.iter = 1000,
bugs.directory = "D:/PROGRAMLAR/WinBUGS14/")

schools.sim


and this is the winbugs code which must be stored as model.txt in D.

 model {
for (j in 1:J)
{
y[j] ~ dnorm (theta[j], tau.y[j])
theta[j] ~ dnorm (mu.theta, tau.theta)
tau.y[j] <- pow(sigma.y[j], -2)
}
mu.theta ~ dnorm (0.0, 1.0E-6)
tau.theta <- pow(sigma.theta, -2)
sigma.theta ~ dunif (0, 1000)
}

• You said: "I asked this question on stackoverflow.com, but I couldn't get what I want"; I imagine that your failure to include the 'r' tag on that post really didn't help. When I am there I filter by the r tag, so I didn't even see your post. Jan 13, 2013 at 13:31
• but the output is a R output, bot winbugs even though the model is written is winbugs Jan 13, 2013 at 13:35

I get the following when I run your code on OpenBUGS. You can pick the interval out of the table.

> schools.sim$summary mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff theta[1] 12.163174 7.894141 -1.2702000 7.489750 11.170000 16.422500 32.14200 1.039443 62 theta[2] 9.143817 6.461938 -4.0285250 5.101000 9.396500 13.250000 21.42525 1.015981 150 theta[3] 7.754204 7.665143 -9.3537750 3.583250 8.474500 12.585000 21.12100 1.015841 360 theta[4] 8.812925 6.602441 -4.5149500 4.506500 9.231000 13.270000 20.40100 1.027968 110 theta[5] 6.754827 6.859304 -8.1867000 2.302000 7.487500 11.380000 17.67000 1.017059 410 theta[6] 7.265947 7.232761 -8.6278000 2.745500 8.155000 11.787500 18.89025 1.022475 190 theta[7] 11.501118 6.360316 -0.2822100 7.494500 11.165000 15.692500 25.01250 1.054659 42 theta[8] 9.684910 7.605208 -4.7286250 5.092000 9.574000 14.362500 25.14250 1.019421 140 mu.theta 9.182670 5.203911 -1.1587500 5.822500 9.265000 12.512500 18.18775 1.029443 88 sigma.theta 5.929714 5.622338 0.2344007 1.685247 4.395495 8.472499 20.17768 1.065369 50 deviance 60.671273 2.242092 57.2300000 59.230000 60.075000 61.930000 65.59050 1.034880 120 > schools.sim$summary[23]
[1] -1.2702


Here is the code I ran.

library(R2OpenBUGS)

data(schools)
J <- nrow(schools)
y <- schools$estimate sigma.y <- schools$sd
data <- list ("J", "y", "sigma.y")

model <- function() {
for (j in 1:J)
{
y[j] ~ dnorm (theta[j], tau.y[j])
theta[j] ~ dnorm (mu.theta, tau.theta)
tau.y[j] <- pow(sigma.y[j], -2)
}
mu.theta ~ dnorm (0.0, 1.0E-6)
tau.theta <- pow(sigma.theta, -2)
sigma.theta ~ dunif (0, 1000)
}
write.model(model,"model.txt")

inits <- function(){
list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),
sigma.theta = runif(1, 0, 100))
}

schools.sim <- bugs(data, inits, model.file = "model.txt",
parameters = c("theta", "mu.theta", "sigma.theta"),
n.chains = 3, n.iter = 1000)

schools.sim$summary schools.sim$summary[23]

• I could not be more grateful. schools.sim\$summary[,3] gives exactly what I want. Thank you very much Jan 13, 2013 at 19:02