POISSON GAMMA MODEL WITH GA(1,1) PRIOR FOR ALPHA model { for (i in 1 : N) { Obs[i] ~ dpois(mu[i]) mu[i] <- Exp[i]*exp(beta0)*theta[i] RR[i] <- exp(beta0)*theta[i] theta[i] ~ dgamma(alpha,alpha) thresh[i] <- step(RR[i]-1.2) } # alpha ~ dgamma(1,1) beta0 ~ dflat() # Functions of interest: sigma.theta <- sqrt(1/alpha) # standard deviation of non-spatial base <- exp(beta0) } POISSON GAMMA MODEL WITH LOGNORMAL PRIOR FOR ALPHA model { for (i in 1 : N) { Obs[i] ~ dpois(mu[i]) mu[i] <- Exp[i]*exp(beta0)*theta[i] RR[i] <- exp(beta0)*theta[i] theta[i] ~ dgamma(alpha,alpha) thresh[i] <- step(RR[i]-1.2) } # lalpha ~ dnorm(2.30,1.40) alpha <- exp(lalpha) beta0 ~ dflat() # Functions of interest: sigma.theta <- sqrt(1/alpha) # standard deviation of non-spatial base <- exp(beta0) } DATA list(N=88,Obs = c(14, 56, 26, 57, 21, 22, 67, 18, 149, 9, 20, 96, 72, 24, 57, 19, 26, 922, 31, 22, 20, 47, 43, 16, 518, 9, 26, 28, 50, 24, 537, 29, 20, 14, 8, 22, 11, 8, 14, 17, 72, 35, 137, 46, 79, 24, 144, 283, 15, 176, 52, 38, 14, 12, 40, 8, 320, 11, 10, 59, 4, 23, 8, 14, 22, 21, 50, 17, 10, 62, 39, 31, 65, 40, 19, 214, 316, 138, 48, 6, 9, 4, 58, 39, 31, 16, 30, 14), Exp = c(15.6783572590858, 62.7864814735788, 26.9533831342569, 59.4483984831333, 25.7109433798694, 24.7643185072109, 52.4373940442046, 19.0822783997527, 129.573778991062, 14.7673346492981, 19.8714574471399, 88.813440086054, 58.7458370764702, 19.7219617366772, 64.5352773524708, 22.2022020547338, 28.5481364466407, 991.173894894725, 31.8040584117781, 19.2049169363729, 26.9133169021863, 45.4659930457387, 50.607787907297, 16.5725520524927, 436.336651157096, 20.1397822888456, 17.2288588864002, 36.2719881107168, 59.0035669254572, 23.9083484859730, 515.2401890566, 34.4611904078476, 17.3406234012350, 10.6547788929387, 15.8400894819940, 23.1232551697127, 14.9160851958418, 13.7311229474182, 29.2640064650174, 18.475080071196, 56.5694958616261, 27.8595192888571, 111.035761795778, 35.5092071123234, 63.9874805541996, 26.1334327270453, 136.298893446484, 262.572434045609, 16.8853349505498, 189.590494143472, 33.8340715399702, 48.5487421987463, 14.5851723617589, 20.9717958517554, 49.7476931663953, 9.8378440745325, 325.426228629749, 9.01477682138723, 13.9236704583501, 46.8912245269686, 7.05604611248674, 26.4089073473897, 10.4540187756399, 16.3808734375746, 23.2946363002087, 13.9867748189397, 58.7640507057366, 22.0490730860987, 17.3890102191324, 71.537511359622, 40.2347204094593, 33.734192243803, 48.7277300305991, 33.9443040768346, 22.2410416503075, 218.642971120352, 303.805530968110, 137.561395276257, 51.7864997608315, 16.6699093850902, 17.8752820104834, 6.9024420887517, 45.7154858000085, 35.7508074400044, 49.058246331259, 20.264821102368, 48.1549342081627, 13.0663927624113)) INITIAL ESTIMATES FOR POISSON GAMMA MODEL WTH GAMMA PRIOR FOR ALPHA list(beta0=0,alpha=1,theta=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)) INITIAL ESTIMATES FOR POISSON GAMMA MODEL WITH LOGNORMAL PRIOR FOR ALPHA list(beta0=0,lalpha=0,theta=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)) POISSON-LOGNORMAL NON-SPATIAL MODEL model { for (i in 1 : N) { Obs[i] ~ dpois(mu[i]) log(mu[i]) <- log(Exp[i]) + beta0 + V[i] RR[i] <- exp(beta0 + V[i]) V[i] ~ dnorm(0,tau.V) thresh[i] <- step(RR[i]-1.2) } # The gamma prior corresponds to df=2, q=0.95, R=log 2. tau.V ~ dgamma(1,0.0260) beta0 ~ dflat() # Functions of interest: sigma.V <- sqrt(1/tau.V) # standard deviation of non-spatial RRRlo <- exp(-1.96*sigma.V) RRRhi <- exp(1.96*sigma.V) } INITIAL ESTIMATES FOR POISSON LOGNORMAL MODEL list(beta0=0,tau.V=1,V=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)) JOINT SPATIAL MODEL + LOGNORMAL NON-SPATIAL MODEL model { for (i in 1 : N) { Obs[i] ~ dpois(mu[i]) log(mu[i]) <- log(Exp[i]) + beta0 + V[i]+ U[i] RR[i] <- exp(beta0 + V[i] + U[i]) V[i] ~ dnorm(0,tau.V) thresh[i] <- step(RR[i]-1.2) mean[i] <- 0 } # Multivariate prior distribution for spatial random effects: U[1:N] ~ spatial.exp(mean[], xm[], ym[], tau.U, phi, 1) tau.T ~ dgamma(1,0.0260) p ~ dbeta(1,1) sigma.U <- sqrt(p/tau.T) sigma.V <- sqrt((1-p)/tau.T) tau.V <- 1/(sigma.V*sigma.V) tau.U <- 1/(sigma.U*sigma.U) # # This prior is derived by assuming that there is a 5% chance that the # correlations die to 0.5 in less that 50km, and a 95% chance that they diE # to 0.5 in less than 200km. # dhalf ~ dlnorm(4.6,0.42) phi <- 0.6931/dhalf beta0 ~ dflat() # Parameters of interest base <- exp(beta0) vratio <- sigma.U*sigma.U/(sigma.U*sigma.U+sigma.V*sigma.V) } DATA list(N=88,Obs = c(14, 56, 26, 57, 21, 22, 67, 18, 149, 9, 20, 96, 72, 24, 57, 19, 26, 922, 31, 22, 20, 47, 43, 16, 518, 9, 26, 28, 50, 24, 537, 29, 20, 14, 8, 22, 11, 8, 14, 17, 72, 35, 137, 46, 79, 24, 144, 283, 15, 176, 52, 38, 14, 12, 40, 8, 320, 11, 10, 59, 4, 23, 8, 14, 22, 21, 50, 17, 10, 62, 39, 31, 65, 40, 19, 214, 316, 138, 48, 6, 9, 4, 58, 39, 31, 16, 30, 14), Exp = c(15.6783572590858, 62.7864814735788, 26.9533831342569, 59.4483984831333, 25.7109433798694, 24.7643185072109, 52.4373940442046, 19.0822783997527, 129.573778991062, 14.7673346492981, 19.8714574471399, 88.813440086054, 58.7458370764702, 19.7219617366772, 64.5352773524708, 22.2022020547338, 28.5481364466407, 991.173894894725, 31.8040584117781, 19.2049169363729, 26.9133169021863, 45.4659930457387, 50.607787907297, 16.5725520524927, 436.336651157096, 20.1397822888456, 17.2288588864002, 36.2719881107168, 59.0035669254572, 23.9083484859730, 515.2401890566, 34.4611904078476, 17.3406234012350, 10.6547788929387, 15.8400894819940, 23.1232551697127, 14.9160851958418, 13.7311229474182, 29.2640064650174, 18.475080071196, 56.5694958616261, 27.8595192888571, 111.035761795778, 35.5092071123234, 63.9874805541996, 26.1334327270453, 136.298893446484, 262.572434045609, 16.8853349505498, 189.590494143472, 33.8340715399702, 48.5487421987463, 14.5851723617589, 20.9717958517554, 49.7476931663953, 9.8378440745325, 325.426228629749, 9.01477682138723, 13.9236704583501, 46.8912245269686, 7.05604611248674, 26.4089073473897, 10.4540187756399, 16.3808734375746, 23.2946363002087, 13.9867748189397, 58.7640507057366, 22.0490730860987, 17.3890102191324, 71.537511359622, 40.2347204094593, 33.734192243803, 48.7277300305991, 33.9443040768346, 22.2410416503075, 218.642971120352, 303.805530968110, 137.561395276257, 51.7864997608315, 16.6699093850902, 17.8752820104834, 6.9024420887517, 45.7154858000085, 35.7508074400044, 49.058246331259, 20.264821102368, 48.1549342081627, 13.0663927624113),xm = c(8, 4, 16, 26, 17, 3, 24, 6, 1, 23, 6, 6, 4, 6, 26, 18, 12, 19, 1, 1, 11, 13, 13, 8, 11, 4, 16, 22, 6, 21, 1, 7, 7, 23, 4, 7, 15, 18, 14, 14, 25, 15, 23, 14, 14, 6, 17, 6, 9, 26, 10, 18, 17, 1, 3, 24, 3, 18, 12, 18, 21, 10, 1, 16, 11, 11, 23, 1, 4, 14, 11, 10, 11, 10, 3, 22, 20, 26, 21, 9, 1, 15, 4, 21, 18, 1, 7, 9), ym = c(3, 19, 20, 27, 7, 18, 13, 3, 8, 17, 14, 12, 5, 8, 19, 15, 20, 25, 14, 24, 15, 25, 10, 9, 12, 26, 2, 25, 10, 13, 6, 21, 18, 15, 24, 5, 8, 17, 23, 4, 17, 16, 28, 1, 14, 16, 24, 27, 12, 21, 18, 22, 5, 18, 13, 10, 10, 9, 17, 12, 10, 26, 22, 10, 10, 4, 23, 10, 21, 19, 7, 24, 13, 22, 15, 20, 23, 24, 17, 15, 20, 6, 8, 8, 20, 26, 24, 20)) INITIAL ESTIMATES FOR THE JOINT MODEL list(beta0=0,tau.T=1,p=0.5,dhalf=10,V=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), U=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)) ICAR SPATIAL MODEL + LOGNORMAL NON-SPATIAL MODEL model { for (i in 1 : N) { Obs[i] ~ dpois(mu[i]) log(mu[i]) <- log(Exp[i]) + beta0 + V[i] + U[i] RR[i] <- exp(beta0 + V[i] + U[i]) V[i] ~ dnorm(0,tau.V) thresh[i] <- step(RR[i]-1.2) } # ICAR prior distribution for spatial random effects: U[1:N] ~ car.normal(adj[], weights[], num[], tau.U) for(k in 1:sumNumNeigh) { weights[k] <- 1 } tau.V ~ dgamma(1,0.0260) tau.U ~ dgamma(1,0.0260) sigma.V <- sqrt(1/tau.V) beta0 ~ dflat() sd.U <- sd(U[1:N]) vratio <- sd.U*sd.U/(sd.U*sd.U+sigma.V*sigma.V) } DATA FOR ICAR MODEL list(N=88,Obs = c(14, 56, 26, 57, 21, 22, 67, 18, 149, 9, 20, 96, 72, 24, 57, 19, 26, 922, 31, 22, 20, 47, 43, 16, 518, 9, 26, 28, 50, 24, 537, 29, 20, 14, 8, 22, 11, 8, 14, 17, 72, 35, 137, 46, 79, 24, 144, 283, 15, 176, 52, 38, 14, 12, 40, 8, 320, 11, 10, 59, 4, 23, 8, 14, 22, 21, 50, 17, 10, 62, 39, 31, 65, 40, 19, 214, 316, 138, 48, 6, 9, 4, 58, 39, 31, 16, 30, 14), Exp = c(15.6783572590858, 62.7864814735788, 26.9533831342569, 59.4483984831333, 25.7109433798694, 24.7643185072109, 52.4373940442046, 19.0822783997527, 129.573778991062, 14.7673346492981, 19.8714574471399, 88.813440086054, 58.7458370764702, 19.7219617366772, 64.5352773524708, 22.2022020547338, 28.5481364466407, 991.173894894725, 31.8040584117781, 19.2049169363729, 26.9133169021863, 45.4659930457387, 50.607787907297, 16.5725520524927, 436.336651157096, 20.1397822888456, 17.2288588864002, 36.2719881107168, 59.0035669254572, 23.9083484859730, 515.2401890566, 34.4611904078476, 17.3406234012350, 10.6547788929387, 15.8400894819940, 23.1232551697127, 14.9160851958418, 13.7311229474182, 29.2640064650174, 18.475080071196, 56.5694958616261, 27.8595192888571, 111.035761795778, 35.5092071123234, 63.9874805541996, 26.1334327270453, 136.298893446484, 262.572434045609, 16.8853349505498, 189.590494143472, 33.8340715399702, 48.5487421987463, 14.5851723617589, 20.9717958517554, 49.7476931663953, 9.8378440745325, 325.426228629749, 9.01477682138723, 13.9236704583501, 46.8912245269686, 7.05604611248674, 26.4089073473897, 10.4540187756399, 16.3808734375746, 23.2946363002087, 13.9867748189397, 58.7640507057366, 22.0490730860987, 17.3890102191324, 71.537511359622, 40.2347204094593, 33.734192243803, 48.7277300305991, 33.9443040768346, 22.2410416503075, 218.642971120352, 303.805530968110, 137.561395276257, 51.7864997608315, 16.6699093850902, 17.8752820104834, 6.9024420887517, 45.7154858000085, 35.7508074400044, 49.058246331259, 20.264821102368, 48.1549342081627, 13.0663927624113), list( num = c(4, 5, 7, 3, 6, 7, 5, 4, 4, 5, 6, 5, 4, 6, 4, 6, 6, 6, 6, 4, 6, 3, 5, 6, 6, 3, 4, 6, 6, 6, 3, 7, 7, 5, 7, 6, 6, 6, 7, 6, 4, 7, 3, 3, 7, 5, 5, 4, 7, 4, 6, 5, 3, 3, 5, 3, 7, 5, 5, 6, 6, 3, 3, 6, 6, 5, 6, 3, 7, 5, 7, 5, 4, 6, 5, 8, 6, 4, 6, 7, 5, 6, 6, 4, 5, 3, 7, 5 ), adj = c( 73, 66, 36, 8, 81, 69, 33, 32, 6, 85, 70, 52, 47, 42, 39, 38, 78, 43, 28, 84, 82, 64, 58, 53, 37, 81, 75, 54, 46, 33, 19, 2, 61, 56, 41, 34, 30, 36, 14, 13, 1, 83, 68, 57, 31, 79, 76, 41, 34, 15, 80, 75, 55, 49, 46, 12, 57, 55, 49, 29, 11, 83, 31, 14, 8, 83, 36, 29, 24, 13, 8, 76, 50, 41, 10, 79, 60, 45, 42, 38, 30, 88, 74, 70, 59, 51, 39, 77, 67, 52, 47, 43, 28, 75, 68, 57, 55, 54, 6, 86, 69, 63, 35, 80, 59, 51, 45, 42, 25, 72, 47, 39, 65, 64, 45, 37, 25, 71, 65, 49, 36, 29, 14, 80, 65, 49, 45, 23, 21, 86, 48, 35, 82, 53, 44, 40, 78, 77, 67, 43, 18, 4, 83, 57, 49, 24, 14, 12, 79, 61, 60, 34, 16, 7, 83, 13, 9, 88, 87, 74, 69, 35, 33, 2, 88, 80, 51, 46, 32, 6, 2, 79, 41, 30, 10, 7, 87, 86, 69, 48, 32, 26, 20, 71, 66, 24, 14, 8, 1, 82, 71, 65, 64, 23, 5, 85, 79, 76, 42, 16, 3, 74, 72, 70, 47, 22, 17, 3, 82, 73, 71, 66, 44, 27, 34, 15, 10, 7, 70, 59, 45, 38, 21, 16, 3, 28, 18, 4, 73, 40, 27, 64, 60, 42, 25, 23, 21, 16, 80, 75, 33, 11, 6, 52, 39, 22, 18, 3, 87, 62, 35, 26, 80, 65, 29, 25, 24, 12, 11, 78, 76, 67, 15, 88, 80, 59, 33, 21, 17, 85, 77, 47, 18, 3, 82, 27, 5, 81, 19, 6, 75, 57, 19, 12, 11, 84, 61, 7, 83, 68, 55, 29, 19, 12, 9, 84, 64, 61, 60, 5, 70, 51, 42, 21, 17, 64, 61, 58, 45, 30, 16, 84, 60, 58, 56, 30, 7, 87, 72, 48, 81, 69, 20, 60, 58, 45, 37, 23, 5, 71, 49, 37, 25, 24, 23, 73, 71, 40, 36, 1, 78, 77, 76, 50, 28, 18, 57, 19, 9, 87, 81, 63, 35, 32, 20, 2, 59, 42, 39, 17, 3, 82, 66, 65, 40, 37, 36, 24, 87, 74, 62, 39, 22, 66, 44, 40, 1, 88, 87, 72, 39, 32, 17, 55, 46, 19, 11, 6, 85, 79, 77, 67, 50, 38, 15, 10, 85, 76, 67, 52, 28, 18, 67, 50, 28, 4, 76, 38, 34, 30, 16, 10, 51, 49, 46, 33, 25, 21, 11, 69, 63, 54, 6, 2, 71, 53, 40, 37, 27, 5, 57, 31, 29, 14, 13, 9, 61, 58, 56, 5, 77, 76, 52, 38, 3, 35, 26, 20, 74, 72, 69, 62, 48, 35, 32, 74, 51, 33, 32, 17 ), sumNumNeigh = 460)) INITIAL ESTIMATES FOR ICAR MODEL list(beta0=0,tau.U=1.0,tau.V=1,V=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), U=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))