Extracting values from R output of bugs command [closed]

Question

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

Answer 1

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]