The example I'm attempting to replicate is the varying intercept model
taken from the ARM book (http://www.stat.columbia.edu/~gelman/arm/).
The example and data can be found here:
http://www.stat.columbia.edu/~gelman/arm/examples/radon/
I'm attempting to replicate the bug file "radon.1.bug" in which
intercepts for each county take their prior from the mean of the
global dataset (mu.a).
The concept that I'm struggling with in pymc is how to define the
linkage between the a[i]'s and mu.a.
Any help is appreciated.
-Whit
warmstrong@research:~/scripts/gelman.examples/radon$ cat radon.1.bug
# Bugs code for multilevel model for radon
# with bsmt as an individual predictor
# varying-intercept model
model {
for (i in 1:n){
y[i] ~ dnorm (y.hat[i], tau.y)
y.hat[i] <- a[county[i]] + b*x[i]
}
b ~ dnorm (0, .0001)
tau.y <- pow(sigma.y, -2)
sigma.y ~ dunif (0, 100)
for (j in 1:J){
a[j] ~ dnorm (mu.a, tau.a)
}
mu.a ~ dnorm (0, .0001)
tau.a <- pow(sigma.a, -2)
sigma.a ~ dunif (0, 100)