Taglibro de H

library(MASS)
data(bacteria)

n <- nrow(bacteria)
yInt <- as.integer(bacteria$y == "y")
week <- bacteria$week
trtdrug <- as.integer(bacteria$trt == "drug")
trtdrugplus <- as.integer(bacteria$trt == "drug+")
id <- as.integer(bacteria$ID)
n.id <- length(levels(bacteria$ID))

model <- "
model {
  for (i in 1:n) {
    yInt[i] ~ dbern(p[i])
    logit(p[i]) <- b + b.week * week[i] +
                   b.drug * trtdrug[i] +
                   b.drugplus * trtdrugplus[i] +
                   e[id[i]]
  }
  for (j in 1:n.id) {
    e[j] ~ dnorm(0.0, tau)
  }
  b ~ dnorm(0.0, 1.0E-4)
  b.week ~ dnorm(0.0, 1.0E-4)
  b.drug ~ dnorm(0.0, 1.0E-4)
  b.drugplus ~ dnorm(0.0, 1.0E-4)
  tau ~ dgamma(1.0E-3, 1.0E-3)
  sigma <- 1/sqrt(tau)
}
"
cat(model, file = "model.bug")

vars <- c("b", "b.week", "b.drug", "b.drugplus", "sigma")
data <- c("yInt", "week", "trtdrug", "trtdrugplus",
          "id", "n", "n.id")
init1 <- list(.RNG.seed = 123,
              .RNG.name = "base::Mersenne-Twister",
              b = -10, b.week = -3, b.drug = 0,
              b.drugplus = 10, tau = 1)
init2 <- list(.RNG.seed = 1234,
              .RNG.name = "base::Mersenne-Twister",
              b = -3, b.week = 0, b.drug = 3,
              b.drugplus = 5, tau = 10)
init3 <- list(.RNG.seed = 12345,
              .RNG.name = "base::Mersenne-Twister",
              b = 10, b.week = -5, b.drug = 10,
              b.drugplus = 10, tau = 0.1)
inits <- list(init1, init2, init3)

library(R2jags)
fit <- jags(data = data, inits = inits,
            parameters.to.save = vars,
            model.file = "model.bug",
            n.chain = 3, n.iter = 2000)
results <- update(fit, n.iter = 30000, thin = 30)

> print(results, digits = 3)
Inference for Bugs model at "model.bug", fit using jags,
 3 chains, each with 30000 iterations (first 2000 discarded)
 n.sims = 90000 iterations saved
              mean     sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff
b            3.411  0.753   2.156   2.882   3.337   3.858   5.093 1.002  1500
b.drug      -1.439  0.774  -3.099  -1.898  -1.396  -0.934  -0.021 1.002  2000
b.drugplus  -0.886  0.765  -2.489  -1.360  -0.857  -0.382   0.550 1.002  2400
b.week      -0.154  0.055  -0.265  -0.190  -0.153  -0.116  -0.051 1.001  5900
deviance   171.779 13.197 150.079 162.463 170.280 179.291 204.829 1.004  1200
sigma        1.402  0.546   0.224   1.059   1.378   1.731   2.560 1.052   330

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 86.9 and DIC = 258.7
DIC is an estimate of expected predictive error (lower deviance is better).