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

          }