model {
# Set up 'data' to define spatial dependence structure
# =====================================
for(i in 1 : N) {
m[i] <- 1/E[i] # scaling factor for variance in each cell
}
# The vector
C[]
required as input into the car.proper distribution is a vector
# respresention of the weight matrix with elements Cij. The first J1 elements of the
C[]
# vector contain the weights for the J1 neighbours of area i=1; the (J1+1) to J2
# elements of the
C[]
vector contain the weights for the J2 neighbours of area i=2;
# etc. To set up this vector, we need to define a variable
cumsum
, which gives the
# values of J1, J2, etc.; we then set up an index matrix
pick[,]
with N columns
# corresponding to the i=1,...,N areas, and with the same number of rows as there are
# elements in the
C[]
vector (i.e.
sumNumNeigh
). The elements
#C[ (cumsum[i]+1):cumsum[i+1] ]
correspond to
# the set of weights Cij associated with area i, and so we set up ith column of the
# matrix
pick[,]
to have a 1 in all the rows k for which
#
cumsum[i] < k <= cumsum[i+1]
, and 0's elsewhere.
# For example, let
N
=4 and
cumsum
=c(0,3,5,6,8), so area i=1 has 3 neighbours, area
# i=2 has 2 neighbours, area i=3 has 1 neighbour and area i=4 has 2 neighbours. The
# the matrix
pick[,]
is:
#
pick
# 1, 0, 0, 0,
# 1, 0, 0, 0,
# 1, 0, 0, 0,
# 0, 1, 0, 0,
# 0, 1, 0, 0,
# 0, 0, 1, 0,
# 0, 0, 0, 1,
# 0, 0, 0, 1,
#
# We can then use the inner product (
inprod(,)
) function in WinBUGS and the kth
# row of
pick
to select which area corresponds to the kth element in the vector
C[]
;
# likewise, we can use
inprod(,)
# and the ith column of
pick
to select the elements of
C[]
which correspond to area i.
#
# Note: this way of setting up the C vector is somewhat convoluted!!!! In future
# versions, we hope the GeoBUGS adjacency matrix tool will be able to dump out the
# relevant vectors required. Alternatively, the C vector could be created using another
# package (e.g. Splus) and read into WinBUGS as data.
#
cumsum[1] <- 0
for(i in 2:(N+1)) {
cumsum[i] <- sum(num[1:(i-1)])
}
for(k in 1 : sumNumNeigh) {
for(i in 1:N) {
pick[k,i] <- step(k - cumsum[i] - epsilon) * step(cumsum[i+1] - k)
# pick[k,i] = 1 if cumsum[i] < k <= cumsum[i=1]; otherwise, pick[k,i] = 0
}
C[k] <- sqrt(E[adj[k]] / inprod(E[], pick[k,])) # weight for each pair of neighbours
}
epsilon <- 0.0001
# Model
# =====
# Likelihood
for (i in 1 : N) {
O[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + S[i]
# Area-specific relative risk
RR[i] <- exp(S[i])
theta[i] <- alpha
}
# Proper CAR prior distribution for spatial random effects:
S[1:N] ~ car.proper(theta[], C[], adj[], num[], m[], prec, gamma)
# Other priors:
alpha ~ dnorm(0, 0.0001)
# prior on precision
prec ~ dgamma(0.5, 0.0005)
v <- 1/prec # variance
sigma <- sqrt(1 / prec) # standard deviation
gamma.min <- min.bound(C
[], adj[], num[], m[]
)
gamma.max <- max.bound(C
[], adj[], num[], m[]
)
gamma ~ dunif(gamma.min, gamma.max)
}