# 
# Seizure example in Section 9.6
# Models 1 and 2 in Table 9.5
#
model
{
for (i in 1:n){
    for (j in 1:k){
        Y[i,j] ~ dpois(mu[i,j])
        log(mu[i,j]) <- log(t[j])+beta0+beta1*x1[i]+beta2*x2[j]+
                                beta3*x1[i]*x2[j]+b[i]
    }
    b[i] ~ dnorm(0,tau)
}
tau ~ dgamma(1,0.26)
# tau ~ dgamma(2,1.376) # Model 2 Prior
sigma <- sqrt(1/tau)
beta0 ~ dflat()
beta1 ~ dflat()
beta2 ~ dflat()
beta3 ~ dflat()
}
#
# Initial estimates     
#          
list(b=c(0.1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2), 
beta0=1, beta1=0.05, beta2=0, beta3=0, tau=1)
#
# Data
#
list(k = 5, n = 59, t = c(8,2,2,2,2), x2 = c(0,1,1,1,1),x1 = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 
     1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),
Y = structure(.Data = c(11,5,3,3,3,11,3,5,3,3,6,2,4,0,5,8,4,4,1,4,66,7,18,9,21,
27,5,2,8,7,12,6,4,0,2,52,40,20,23,12,23,5,6,6,5,10,14,13,6,0,52,26,12,6,22,
33,12,6,8,4,18,4,4,6,2,42,7,9,12,14,87,16,24,10,9,50,11,0,0,5,18,0,0,3,3,
111,37,29,28,29,18,3,5,2,5,20,3,0,6,7,12,3,4,3,4,9,3,4,3,4,17,2,3,3,5,
28,8,12,2,8,55,18,24,76,25,9,2,1,2,1,10,3,1,4,2,47,13,15,13,12,76,11,14,9,8,
38,8,7,9,4,19,0,4,3,0,10,3,6,1,3,19,2,6,7,4,24,4,3,1,3,31,22,17,19,16,
14,5,4,7,4,11,2,4,0,4,67,3,7,7,7,41,4,18,2,5,7,2,1,1,0,22,0,2,4,0,13,5,4,0,3,
46,11,14,25,15,36,10,5,3,8,38,19,7,6,7,7,1,1,2,3,36,6,10,8,8,11,2,1,0,0,
151,102,65,72,63,22,4,3,2,4,41,8,6,5,7,32,1,3,1,5,56,18,11,28,13,24,6,3,4,0,
16,3,5,4,3,22,1,23,19,8,25,2,3,0,1,13,0,0,0,0,12,1,4,3,2),.Dim = c(59,5)))
#
# Model 3 - Poisson lognormal for nugget also
#
model
{
for (i in 1:n){
    for (j in 1:k){
        Y[i,j] ~ dpois(mu[i,j])
        log(mu[i,j]) <- log(t[j])+beta0+beta1*x1[i]+beta2*x2[j]+
                                beta3*x1[i]*x2[j]+b[i]+be[i,j]
        be[i,j] ~ dnorm(0,taue)
    }
    b[i] ~ dnorm(0,tau)
}
taue ~ dgamma(1,0.26)
tau ~ dgamma(1,0.26)
sigma <- sqrt(1/tau)
sigmae <- 1/sqrt(taue)
beta0 ~ dflat()
beta1 ~ dflat()
beta2 ~ dflat()
beta3 ~ dflat()
}