model { # Likelihood for (i in 1 : Nareas) { for (k in 1 : Ndiseases) { Y[i, k] ~ dpois(mu[i, k]) log(mu[i, k]) <- log(E[i, k]) + alpha[k] + eta[i, k] } } for(i in 1:Nareas) { # Define log relative risk in terms of disease-specific (psi) and shared (phi) # random effects # changed order of k and i index for psi (needed because car.normal assumes # right hand index is areas) eta[i, 1] <- phi[i] * delta + psi[1, i] eta[i, 2] <- phi[i] / delta + psi[2, i] } # Spatial priors (BYM) for the disease-specific random effects for (k in 1 : Ndiseases) { for (i in 1 : Nareas) { # convolution prior = sum of unstructured and spatial effects psi[k, i] <- U.sp[k, i] + S.sp[k, i] # unstructured disease-specific random effects U.sp[k, i] ~ dnorm(0, tau.unstr[k]) } # spatial disease-specific effects S.sp[k,1 : Nareas] ~ car.normal(adj[], weights[], num[], tau.spatial[k]) } # Spatial priors (BYM) for the shared random effects for (i in 1:Nareas) { # convolution prior = sum of unstructured and spatial effects phi[i] <- U.sh[i] + S.sh[i] # unstructured shared random effects U.sh[i] ~ dnorm(0, omega.unstr) } # spatial shared random effects S.sh[1:Nareas] ~ car.normal(adj[], weights[], num[], omega.spatial) for (k in 1:sumNumNeigh) { weights[k] <- 1 } # Other priors for (k in 1:Ndiseases) { alpha[k] ~ dflat() tau.unstr[k] ~ dgamma(0.5, 0.0005) tau.spatial[k] ~ dgamma(0.5, 0.0005) } omega.unstr ~ dgamma(0.5, 0.0005) omega.spatial ~ dgamma(0.5, 0.0005) # scaling factor for relative strength of shared component for each disease logdelta ~ dnorm(0, 5.9) # (prior assumes probability that delta^2 is between 1/5 and 5; delta <- exp(logdelta) # lognormal assumption is invariant to which disease is labelled 1 # and which is labelled 2) # ratio (relative risk of disease 1 associated with shared component) to # (relative risk of disease 2 associated with shared component) # - see Knorr-Held and Best (2001) for further details RR.ratio <- pow(delta, 2) # Relative risks and other summary quantities # The GeoBUGS map tool can only map vectors, so need to create separate vector # of quantities to be mapped, rather than an array (i.e. totalRR[i,k] wont work!) for (i in 1 : Nareas) { SMR1[i] <- Y[i,1] / E[i,1] # SMR for disease 1 (oral) SMR2[i] <- Y[i,2] / E[i,2] # SMR for disease 2 (lung) totalRR1[i] <- exp(eta[i,1]) # overall RR of disease 1 (oral) in area i totalRR2[i] <- exp(eta[i,2]) # overall RR of disease 2 (lung) in area i # residulal RR specific to disease 1 (oral cancer) specificRR1[i]<- exp(psi[1,i]) # residulal RR specific to disease 2 (lung cancer) specificRR2[i]<- exp(psi[2,i]) # shared component of risk common to both diseases sharedRR[i] <- exp(phi[i]) # Note that this needs to be scaled by delta or 1/delta if the # absolute magnitude of shared RR for each disease is of interest logsharedRR1[i] <- phi[i] * delta logsharedRR2[i] <- phi[i] /delta } # empirical variance of shared effects (scaled for disease 1) var.shared[1] <- sd(logsharedRR1[])*sd(logsharedRR1[]) # empirical variance of shared effects (scaled for disease 2) var.shared[2] <- sd(logsharedRR2[])*sd(logsharedRR2[]) # empirical variance of disease 1 specific effects var.specific[1] <- sd(psi[1,])*sd(psi[1,]) # empirical variance of disease 2 specific effects var.specific[2] <- sd(psi[2,])*sd(psi[2,]) # fraction of total variation in relative risks for each disease that is explained # by the shared component frac.shared[1] <- var.shared[1] / (var.shared[1] + var.specific[1]) frac.shared[2] <- var.shared[2] / (var.shared[2] + var.specific[2]) }