Analysis_Stage.txt

##########################################
### Model code for the Analysis Stage  ###
##########################################
# In this model, we assume a linear relationship between Y and Z
# numareas = 625 (total number of electoral wards in London)

# eps.mean is the posterior mean of the generalized EPS estimated in the Design stage
# eps.prec is the precision of the generalized EPS estimated in the Design stage
# EPS.mis is a vector of lenght = 625, which includes for the in-sample areas the posterior mean of the generalized EPS 
# and for the out-of-sdample areas the value NA. 
# EPS.tot takes value "EPS" for the in-sample areas and value "EPS.mis" imputed for the out-of-sample areas. In fact, the estimated value of 
# the generalized EPS to be assigned to areal unit i for confounding adjustment is evaluated by using indicator functions implemented via step functions.

# NOTE: for eps.mean and eps.prec we used a trick, as the loop goes from 1 to the total number of areas=625.
# In particular, for the eps.mean we set the values of the out-of-sample areas equal to 0, and for eps.prec we set the 
# values of the out-of-sample areas equal to inverse of the mean of the posterior variance. 
# Importantly, this does not affect the computation, as for the out-of-sample areas the
# value used at each iteration for the generalized EPS is predicted by the imputation model. 

model {

# Uncertainty model  #
for(i in 1:numareas) { 
EPS[i] ~ dnorm(eps.mean[i], eps.prec[i]) 
}

# Imputation model #
for(i in 1:numareas) {
EPS.mis[i] ~ dnorm(mu.eps[i], tau.eps)
mu.eps[i] <- eta + gamma1*X[i] + gamma2*C[i] + theta[i]


# Health outcome model #
Y[i] ~ dpois(mu[i])	

EPS.tot[i] <- EPS[i]*(1-step(ind.mis[i]-1))+EPS.mis[i]*step(ind.mis[i]-1) 

log(mu[i]) <- log(E[i])+beta0+beta1*X[i]+beta2*C[i]+beta3*EPS.tot[i]+V[i]

V[i] ~ dnorm(0.0, tau.V) 

}

# Priors #
eta ~ dflat()  
gamma1 ~ dnorm(0, 0.1)
gamma2 ~ dnorm(0, 0.1)

theta[1:numareas] ~ car.normal(adjSp[],  weights[], numSp[], tau.car)
tau.car ~ dgamma(2,1)
sigma2.car<-1/tau.car

tau.eps ~ dgamma(1,0.1)
sigma2.taueps <- 1/tau.eps

beta0 ~ dnorm(0, 0.01)
beta1 ~ dnorm(0, 0.01)
beta2 ~ dnorm(0, 0.01)
beta3 ~ dnorm(0, 0.01) 

sigma ~ dunif(0,1)
tau.V <- 1/(sigma*sigma) 
}