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 100s) in area i
0DX 0DX # 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] } # (J1)th source (overall intercept; represents background rate) lambda[i,J1] <- beta0 # (J2)th source (effect of observed NO2 covariate in area i) lambda[i,J2] <- beta1 * NO2[i] # Note: if additional covariates have been measured, these # should be included in the same way as NO2, giving terms # lambda[i, J3], lambda[i, J4],.....etc. for (j in 1 : J2) { 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
0DX 0DX # prior pointsover 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 }