##############################################

# Section 10.4.3 Hierarchial model for tables

##############################################

model

{

sigma2 <- 0.65*0.65 # Assume it is known.

for(i in 1:8)

{

for(j in 1:5)

{

y[i,j] <- gpa[(i-1)*5+j]

m[i,j] <- n[(i-1)*5+j]

x1[i,j] <- ACT[(i-1)*5+j]

x2[i,j] <- HSR[(i-1)*5+j]

y[i,j] ~ dnorm(mu[i,j], tau2[i,j])

tau2[i,j] <- m[i,j]/sigma2

mu[i,j] ~dnorm(xb[i,j], sigma2pi.inv)

xb[i,j] <- b[1] + b[2]*x1[i,j] + b[3]*x2[i,j]

}

}

sigma2pi.inv ~ dgamma(8, 0.01)

b[1] ~ dnorm( 0, 0.0001)

b[2] ~ dnorm( 0, 0.0001)

b[3] ~ dnorm( 0, 0.0001)

sigma2pi <- 1/sigma2pi.inv

# Posterior prediction of Pr(GPA > 2.5|y)

for(i in 1:8)

{

for(j in 1:5)

{

prob[i,j] <- 1-phi( (2.5-mu[i,j])/sqrt(sigma2) )

}

}

}

# Initial values. Initialize other values using "gen inits".

list( sigma2pi.inv = 1, b = c(0, 0, 0) )

# Data

gpa[] n[] HSR[] ACT[]

2.64 8 95 17

3.10 15 95 20

3.01 78 95 23

3.07 182 95 26

3.34 166 95 29

2.24 20 85 17

2.63 71 85 20

2.74 168 85 23

2.76 178 85 26

2.91 91 85 29

2.43 40 75 17

2.47 116 75 20

2.64 180 75 23

2.73 133 75 26

2.47 46 75 29

2.31 34 65 17

2.37 93 65 20

2.32 124 65 23

2.24 101 65 26

2.31 19 65 29

2.04 41 55 17

2.20 73 55 20

2.01 62 55 23

2.43 58 55 26

2.38 9 55 29

1.88 19 45 17

1.82 25 45 20

1.84 36 45 23

2.12 49 45 26

2.05 16 45 29

1.86 8 35 17

2.28 9 35 20

1.67 15 35 23

1.89 29 35 26

1.79 9 35 29

1.70 4 25 17

1.65 5 25 20

1.51 9 25 23

1.67 11 25 26

2.33 1 25 29

END