			model {

			# Likelihood
				for (i in 1 : I) {
			# convolution likelihood: 
			#          Y[i] ~ Poisson( sum_ j  { mu[i, j] } )
			# where mu[i, j] = pop[i] * lambda[i, j] and
			# pop[i] = (standardised) population (in 100's) in area i
			# lambda[i, j] = disease rate in area i attributed to jth 'source', where
			#                     the sources include both observed and latent covariates.
					Y[i] ~ dpois.conv(mu[i, ])	
                                            
				for (j in 1 : J) { 
			# rescale kernel matrix
					k1[i, j] <- ker[i, j] / prec
			# jth 'source' (jth latent spatial grid cell)
					lambda[i, j] <- beta2 * gamma[j] * k1[i, j]
				}
			# (J+1)th 'source' (overall intercept; represents background rate)
				lambda[i,J+1] <- beta0
			# (J+2)th 'source' (effect of observed NO2 covariate in area i)
				lambda[i,J+2] <- beta1 * NO2[i]
			# Note: if additional covariates have been measured, these 
			# should be included in the same way as NO2, giving terms
			# lambda[i, J+3], lambda[i, J+4],.....etc. 
				for (j in 1 : J+2) {                        
					mu[i, j] <- lambda[i, j] * pop[i]
				}
			}
					
			# Priors
			# See Table 22.1 in Best et al. (2000b) for further details.
			
			# Priors for latent spatial grid cell (gamma) parameters:
			#  
			# Assume gamma_j ~ gamma(a, t)  
			# 
			# Prior shape parameter for gamma distn on gamma[j]'s will change in proportion to
			# the size of the latent grid cells, to ensure aggregation consistency. The prior inverse
			# scale parameter for  gamma distn on gamma[j]'s stays constant across different
			# partitions, but depends on the  units used to measure the area of the latent grid cells:
			# here we are assuming # 'area'  is in sq kilometres, and we take tauG = 1 (per km^2)
			# which corresponds to assuming # that our prior guess at the value of the gamma[j]'s
			# is based on 'observing' about 1 * (area of study region = 30*20km^2) = 600
			# 'prior points' over the study region. 
				
			# Note: area of latent grid cells is read in from data file
				alphaG <-  tauG * area
				tauG <- 1          
       
				for (j in 1 : J) {
					gamma[j] ~ dgamma(alphaG, tauG)
				}

			# Priors for beta coefficients:
			#
			# Assume priors on each beta are of the form beta ~ gamma(a, t).
			#
			# Here, parameter of the prior are chosen by setting the prior mean to give an equal 
			# number of cases  allocated to each of the 3 'sources' (baseline, NO2 and latent), i.e.
			# prior mean   =  a / t  =  a * 3 * Xbar / Ybar,
			# where Ybar is the overall disease rate in number of cases per unit population,
			# and Xbar is the population-weighted mean of covariate X. 
			# 
			# We also take shape parameter a = 0.575, since this gives  the ratio of the 90th : 10th
			# percentile of the prior distn = 100, which is suitably diffuse.
			
				YBAR <- sum(Y[])/sum(pop[])	# mean number of cases per unit population
				NO2BAR <- inprod(NO2[], pop[]) / sum(pop[])	# population-weighted mean NO2
				
			# Shape parameter for Intercept (beta0)
				alpha0<- 0.575
			# Inverse scale parameter for Intercept (beta0)
				tau0 <- 3 * alpha0 / YBAR
			# Shape parameter for effect of NO2 (beta1)
				alpha1<- 0.575
			# Inverse scale parameter for effect of NO2 (beta1)
				tau1 <- 3 * alpha1 * NO2BAR / YBAR
			# Shape parameter for Latent coefficient (beta2)
				alpha2 <- 0.575
			# Inverse scale parameter for Latent coefficient (beta2)
				tau2 <- 3 * alpha2	 / YBAR

				beta0 ~ dgamma(alpha0, tau0)
				beta1 ~ dgamma(alpha1, tau1)
				beta2 ~ dgamma(alpha2, tau2)


			# Summary quantities for posterior inference 
				for(i in 1:I) {
			# Overall disease rate in area i  
					RATE[i] <- sum(mu[i, ])/pop[i]
			# Rate associated with latent spatial factor in area i 
					LATENT[i] <- beta2 * inprod(k1[i,], gamma[])
				}

			# Expected rate (number of cases per unit population) attributed to each source 
			# (see Table 22.2 in Best et al (2000b) )

			# Background (overall baseline) rate per 100 population
			rate.base <- beta0

			# Number of cases per 100 population attributed
			# to NO2 = beta1 * population-weighted mean value
			# of NO2 across the study region. 
				rate.no2 <- beta1 * inprod(pop[], NO2[]) / sum(pop[])
				
			# Number of cases per 100 population attributed to 
			# latent sources = beta2 * population-weighted mean
			# value of the latent spatial factor (sum_j k_ij gamma_j)
				rate.latent <- inprod(pop[], LATENT[]) / sum(pop[])
     
			# Percentage of cases attributed to each source:
				Total <- rate.base + rate.no2 + rate.latent   
			# % cases attributed to baseline risk
				PC.base <- rate.base / Total * 100
			# % cases attributed to NO2 exposure
				PC.no2 <- rate.no2 / Total * 100
			# % cases attributed to latent spatial factors
				PC.latent <- rate.latent / Total * 100

			}
