MODEL 1 is non-spatial only and no covariate
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + beta0 + V[i]
RR[i] <- exp(beta0 + V[i])
V[i] ~ dnorm(0,tau.V)
}
# 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)
}
DATA 1
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8))
INIT 1
list(tau.V = 100, beta0 = 0,
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))
MODEL 2 is non-spatial but with linear covariate
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + beta0 + beta1*X[i] + V[i]
RR[i] <- exp(beta0 + beta1*X[i]+V[i])
V[i] ~ dnorm(0,tau.V)
}
# The gamma prior corresponds to a log t prior for the residual relative
# risks with df=2, q=0.95, R=log 2.
tau.V ~ dgamma(1,0.0260)
beta0 ~ dflat()
beta1 ~ 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)
}
DATA 2
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10))
INIT 2
list(ltau.V = 100, beta0 = 0, beta1 = 0,
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))
MODEL 3 is non-spatial but with cubic covariate
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
X1c[i] <- X[i]-mean(X[1:N])
X2c[i] <- X1c[i]*X1c[i]
X3c[i] <- X1c[i]*X1c[i]*X1c[i]
log(mu[i]) <- log(E[i]) + beta0 + beta1*X1c[i] + beta2*X2c[i] + beta3*X3c[i] + V[i]
RR[i] <- exp(beta0 + beta1*X1c[i] + beta2*X2c[i]+ beta3*X3c[i] + V[i])
V[i] ~ dnorm(0,tau.V)
}
# The gamma prior corresponds to df=2, q=0.95, R=log 2.
tau.V ~ dgamma(1,0.0260)
beta0 ~ dflat()
beta1 ~ dflat()
beta2 ~ dflat()
beta3 ~ 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)
}
DATA 3
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10))
INIT 3
list(tau.V = 100, beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0,
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))
MODEL 4 is non-spatial and spatial ICAR with linear covariate
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + beta0 + beta1*X[i] + V[i] + U[i]
RR[i] <- exp(beta0 + beta1*X[i] + V[i] + U[i])
V[i] ~ dnorm(0,tau.V)
}
# 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
}
# we assume that r = s_u^2/omega_u^2 so that tau.U = r/s_u^2
r <-
tau.V ~ dgamma(1,0.0260)
tau.U ~ dgamma(1,0.0260*r)
# Other priors:
beta0 ~ dflat()
beta1 ~ dflat()
sigma.V <- sqrt(1 / tau.V) # standard deviation of unstructured
sigma.U <- sd(U[1:N])
vratio <- sigma.U^2/(sigma.U^2+sigma.V^2)
}
DATA 4
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10))
INIT 4
list(tau.V = 0, tau.U=0,beta0 = 0, beta1 = 0,
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),
U=c(0,0,0,0,0,NA,0,NA,0,0,
NA,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0))
MODEL 5: Spatial exp correlation model + non-spatial + linear covariate +
independent variance priors.
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + beta0 + beta1*X[i] + V[i]+ U[i]
RR[i] <- exp(beta0 + beta1*X[i] + V[i] + U[i])
V[i] ~ dnorm(0,tau.V)
mean[i] <- 0
}
# Multivariate prior distribution for spatial random effects:
U[1:N] ~ spatial.exp(mean[], xm[], ym[], tau.U, phi, 1)
tau.V ~ dgamma(1,0.0260)
tau.U ~ dgamma(1,0.0260)
#
# This prior is derived by assuming that there is a 5% chance that the
# correlations die to 0.5 in less that 5km, and a 95% chance that they die
# to 0.5 in less than 100km.
#
dhalf ~ dnorm(3.107,0.9106)
phi <- 0.6931/dhalf
beta0 ~ dflat()
beta1 ~ dflat()
# Parameters of interest
base <- exp(beta0)
RRx <- exp(beta1)
sigma.V <- sqrt(1/tau.V) # standard deviation for non-spatial
sigma.U <- sqrt(1/tau.U) # standard deviation for spatial
vratio <- sigma.U^2/(sigma.U^2+sigma.V^2)
}
DATA 5
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10),
xm = c(
162.1894, 385.7761, 293.9555, 377.9338, 220.6786,
340.1739, 324.9915, 442.2445, 194.5176, 367.6924,
112.8916, 247.7566, 289.5922, 227.9563, 342.3574,
351.3505, 280.4916, 341.6081, 249.6855, 359.5902,
348.7138, 388.7655, 180.4228, 295.4908, 333.1159,
312.0605, 290.1701, 359.4153, 291.3727, 303.4219,
257.4402, 264.9711, 336.4464, 258.0319, 227.1801,
234.5294, 218.3428, 279.1010, 235.0805, 254.1736,
250.8301, 287.1202, 292.3773, 288.0333, 320.5682,
257.8758, 276.9737, 281.9644, 267.8444, 342.226,
274.8713, 257.8069, 265.5934, 267.8921, 321.4991,
322.1780),
ym =c(
834.7496, 852.3782, 946.0722, 650.501,
870.9356, 1015.154, 842.0317, 1168.904, 781.3746,
828.219, 903.1592, 924.9536, 842.3052, 561.1628,
713.0808, 792.1617, 801.0356, 628.6406, 825.8545,
610.6554, 760.2982, 812.7655, 699.6693, 635.7658,
701.8189, 691.102, 586.6673, 669.4746, 746.2605,
670.1395, 605.9585, 568.3428, 658.671, 716.452,
598.2521, 668.0481, 641.4785, 670.285, 697.044,
677.589, 657.4675, 680.7535, 699.3761, 665.2905,
671.6064, 631.046, 640.8285, 654.6629, 666.7073,
736.4561, 678.8585, 683.7104, 646.5754, 682.2943,
640.1429, 589.9408))
INIT 5
list(tau.V = 1, tau.U=1,beta0 = 0, beta1 = 0,dhalf =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),
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))
MODEL 6: Spatial exp correlation model + non-spatial + linear covariate +
dependent variance priors.
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + beta0 + beta1*X[i] + V[i]+ U[i]
RR[i] <- exp(beta0 + beta1*X[i] + V[i] + U[i])
V[i] ~ dnorm(0,tau.V)
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 5km, and a 95% chance that they die
# to 0.5 in less than 100km.
#
dhalf ~ dlnorm(3.107,0.9106)
phi <- 0.6931/dhalf
beta0 ~ dflat()
beta1 ~ dflat()
# Parameters of interest
base <- exp(beta0)
RRx <- exp(beta1)
vratio <- sigma.U*sigma.U/(sigma.U*sigma.U+sigma.V*sigma.V)
}
DATA 6
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10),
xm = c(
162.1894, 385.7761, 293.9555, 377.9338, 220.6786,
340.1739, 324.9915, 442.2445, 194.5176, 367.6924,
112.8916, 247.7566, 289.5922, 227.9563, 342.3574,
351.3505, 280.4916, 341.6081, 249.6855, 359.5902,
348.7138, 388.7655, 180.4228, 295.4908, 333.1159,
312.0605, 290.1701, 359.4153, 291.3727, 303.4219,
257.4402, 264.9711, 336.4464, 258.0319, 227.1801,
234.5294, 218.3428, 279.1010, 235.0805, 254.1736,
250.8301, 287.1202, 292.3773, 288.0333, 320.5682,
257.8758, 276.9737, 281.9644, 267.8444, 342.226,
274.8713, 257.8069, 265.5934, 267.8921, 321.4991,
322.1780),
ym =c(
834.7496, 852.3782, 946.0722, 650.501,
870.9356, 1015.154, 842.0317, 1168.904, 781.3746,
828.219, 903.1592, 924.9536, 842.3052, 561.1628,
713.0808, 792.1617, 801.0356, 628.6406, 825.8545,
610.6554, 760.2982, 812.7655, 699.6693, 635.7658,
701.8189, 691.102, 586.6673, 669.4746, 746.2605,
670.1395, 605.9585, 568.3428, 658.671, 716.452,
598.2521, 668.0481, 641.4785, 670.285, 697.044,
677.589, 657.4675, 680.7535, 699.3761, 665.2905,
671.6064, 631.046, 640.8285, 654.6629, 666.7073,
736.4561, 678.8585, 683.7104, 646.5754, 682.2943,
640.1429, 589.9408))
INIT 6
list(tau.T = 1, p=0.5,beta0 = 0, beta1 = 0,dhalf =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),
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))
MODEL 7: Spatial exp correlation model + non-spatial + cubic covariate +
dependent variance priors.
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
X1c[i] <- X[i]-mean(X[1:N])
X2c[i] <- X1c[i]*X1c[i]
X3c[i] <- X1c[i]*X1c[i]*X1c[i]
log(mu[i]) <- log(E[i]) + beta0 + beta1*X1c[i] + beta2*X2c[i] +
beta3*X3c[i] + V[i] + U[i]
RR[i] <- exp(beta0 + beta1*X1c[i] + beta2*X2c[i]+
beta3*X3c[i] + V[i] + U[i])
V[i] ~ dnorm(0,tau.V)
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)
# tau.T ~ dgamma(1,0.1399)
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 5km, and a 95% chance that they die
# to 0.5 in less than 100km.
#
# dhalf ~ dlnorm(3.107,0.9106)
#
# This prior is derived by assuming that there is a 5% chance that the
# correlations die to 0.5 in less that 5km, and a 95% chance that they die
# to 0.5 in less than 20km.
#
dhalf ~ dlnorm(2.303,0.4214)
phi <- 0.6931/dhalf
beta0 ~ dflat()
beta1 ~ dflat()
beta2 ~ dflat()
beta3 ~ dflat()
# Parameters of interest
vratio <- sigma.U*sigma.U/(sigma.U*sigma.U+sigma.V*sigma.V)
}
DATA 7
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10),
xm = c(
162.1894, 385.7761, 293.9555, 377.9338, 220.6786,
340.1739, 324.9915, 442.2445, 194.5176, 367.6924,
112.8916, 247.7566, 289.5922, 227.9563, 342.3574,
351.3505, 280.4916, 341.6081, 249.6855, 359.5902,
348.7138, 388.7655, 180.4228, 295.4908, 333.1159,
312.0605, 290.1701, 359.4153, 291.3727, 303.4219,
257.4402, 264.9711, 336.4464, 258.0319, 227.1801,
234.5294, 218.3428, 279.1010, 235.0805, 254.1736,
250.8301, 287.1202, 292.3773, 288.0333, 320.5682,
257.8758, 276.9737, 281.9644, 267.8444, 342.226,
274.8713, 257.8069, 265.5934, 267.8921, 321.4991,
322.1780),
ym =c(
834.7496, 852.3782, 946.0722, 650.501,
870.9356, 1015.154, 842.0317, 1168.904, 781.3746,
828.219, 903.1592, 924.9536, 842.3052, 561.1628,
713.0808, 792.1617, 801.0356, 628.6406, 825.8545,
610.6554, 760.2982, 812.7655, 699.6693, 635.7658,
701.8189, 691.102, 586.6673, 669.4746, 746.2605,
670.1395, 605.9585, 568.3428, 658.671, 716.452,
598.2521, 668.0481, 641.4785, 670.285, 697.044,
677.589, 657.4675, 680.7535, 699.3761, 665.2905,
671.6064, 631.046, 640.8285, 654.6629, 666.7073,
736.4561, 678.8585, 683.7104, 646.5754, 682.2943,
640.1429, 589.9408))
INIT 7
list(tau.T = 1, p=0.5,beta0 = 0, beta1 = 0, beta2 = 0, beta3 =0, dhalf =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),
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))
MODEL 8: Spatial exp correlation model + non-spatial + linear covariate +
dependent variance priors.
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
Xc[i] <- X[i]-mean(X[1:N])
log(mu[i]) <- log(E[i]) + beta0 + beta1*Xc[i] + V[i] + U[i]
RR[i] <- exp(beta0 + beta1*Xc[i] + V[i] + U[i])
V[i] ~ dnorm(0,tau.V)
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)
# tau.T ~ dgamma(1,0.1399)
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 5km, and a 95% chance that they die
# to 0.5 in less than 100km.
#
# dhalf ~ dlnorm(3.107,0.9106)
#
# This prior is derived by assuming that there is a 5% chance that the
# correlations die to 0.5 in less that 5km, and a 95% chance that they die
# to 0.5 in less than 20km.
#
dhalf ~ dlnorm(2.303,0.4214)
phi <- 0.6931/dhalf
beta0 ~ dflat()
beta1 ~ dflat()
# Parameters of interest
base <- exp(beta0)
RRx <- exp(beta1)
vratio <- sigma.U*sigma.U/(sigma.U*sigma.U+sigma.V*sigma.V)
}
DATA 8
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10),
xm = c(
162.1894, 385.7761, 293.9555, 377.9338, 220.6786,
340.1739, 324.9915, 442.2445, 194.5176, 367.6924,
112.8916, 247.7566, 289.5922, 227.9563, 342.3574,
351.3505, 280.4916, 341.6081, 249.6855, 359.5902,
348.7138, 388.7655, 180.4228, 295.4908, 333.1159,
312.0605, 290.1701, 359.4153, 291.3727, 303.4219,
257.4402, 264.9711, 336.4464, 258.0319, 227.1801,
234.5294, 218.3428, 279.1010, 235.0805, 254.1736,
250.8301, 287.1202, 292.3773, 288.0333, 320.5682,
257.8758, 276.9737, 281.9644, 267.8444, 342.226,
274.8713, 257.8069, 265.5934, 267.8921, 321.4991,
322.1780),
ym =c(
834.7496, 852.3782, 946.0722, 650.501,
870.9356, 1015.154, 842.0317, 1168.904, 781.3746,
828.219, 903.1592, 924.9536, 842.3052, 561.1628,
713.0808, 792.1617, 801.0356, 628.6406, 825.8545,
610.6554, 760.2982, 812.7655, 699.6693, 635.7658,
701.8189, 691.102, 586.6673, 669.4746, 746.2605,
670.1395, 605.9585, 568.3428, 658.671, 716.452,
598.2521, 668.0481, 641.4785, 670.285, 697.044,
677.589, 657.4675, 680.7535, 699.3761, 665.2905,
671.6064, 631.046, 640.8285, 654.6629, 666.7073,
736.4561, 678.8585, 683.7104, 646.5754, 682.2943,
640.1429, 589.9408))
INIT 8
list(tau.T = 1, p=0.5,beta0 = 0, beta1 = 0, dhalf =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),
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))
MODEL 9: Spatial exp correlation model + non-spatial + indiv + individual model
dependent variance priors.
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + log( (1-X[i])*exp(alpha0) + X[i]*exp(alpha0+alpha1)) + V[i] + U[i]
RR[i] <- ( (1-X[i])*exp(alpha0) + X[i]*exp(alpha0+alpha1) )*exp(V[i] + U[i])
V[i] ~ dnorm(0,tau.V)
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)
# tau.T ~ dgamma(1,0.1399)
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 5km, and a 95% chance that they die
# to 0.5 in less than 100km.
#
dhalf ~ dlnorm(3.107,0.9106)
#
# This prior is derived by assuming that there is a 5% chance that the
# correlations die to 0.5 in less that 5km, and a 95% chance that they die
# to 0.5 in less than 20km.
#
# dhalf ~ dlnorm(2.303,0.4214)
phi <- 0.6931/dhalf
alpha0 ~ dflat()
alpha1 ~ dflat()
# Parameters of interest
base <- exp(alpha0)
RRx <- exp(alpha1)
vratio <- sigma.U*sigma.U/(sigma.U*sigma.U+sigma.V*sigma.V)
}
DATA 9
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10),
xm = c(
162.1894, 385.7761, 293.9555, 377.9338, 220.6786,
340.1739, 324.9915, 442.2445, 194.5176, 367.6924,
112.8916, 247.7566, 289.5922, 227.9563, 342.3574,
351.3505, 280.4916, 341.6081, 249.6855, 359.5902,
348.7138, 388.7655, 180.4228, 295.4908, 333.1159,
312.0605, 290.1701, 359.4153, 291.3727, 303.4219,
257.4402, 264.9711, 336.4464, 258.0319, 227.1801,
234.5294, 218.3428, 279.1010, 235.0805, 254.1736,
250.8301, 287.1202, 292.3773, 288.0333, 320.5682,
257.8758, 276.9737, 281.9644, 267.8444, 342.226,
274.8713, 257.8069, 265.5934, 267.8921, 321.4991,
322.1780),
ym =c(
834.7496, 852.3782, 946.0722, 650.501,
870.9356, 1015.154, 842.0317, 1168.904, 781.3746,
828.219, 903.1592, 924.9536, 842.3052, 561.1628,
713.0808, 792.1617, 801.0356, 628.6406, 825.8545,
610.6554, 760.2982, 812.7655, 699.6693, 635.7658,
701.8189, 691.102, 586.6673, 669.4746, 746.2605,
670.1395, 605.9585, 568.3428, 658.671, 716.452,
598.2521, 668.0481, 641.4785, 670.285, 697.044,
677.589, 657.4675, 680.7535, 699.3761, 665.2905,
671.6064, 631.046, 640.8285, 654.6629, 666.7073,
736.4561, 678.8585, 683.7104, 646.5754, 682.2943,
640.1429, 589.9408))
INIT 9
list(tau.T = 1, p=0.5,alpha0 = 0, alpha1 = 0, dhalf =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),
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))
MODEL 10: Spatial exp correlation model + non-spatial + linear covariate +
dependent variance priors.
model
{
for (i in 1 : N) {
Y[i] ~ dpois(mu[i])
Xc[i] <- X[i]-mean(X[1:N])
log(mu[i]) <- log(E[i]) + beta0 + beta1*Xc[i] + T[i]
RR[i] <- exp(beta0 + beta1*Xc[i] + T[i])
T[i] <- p[i]*U[i] + (1-p[i])*V[i]
p[i] ~ dbeta(a,b)
V[i] ~ dnorm(0,tau.V)
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)
# tau.T ~ dgamma(1,0.1399)
a ~ dgamma(1,1)
b ~ dgamma(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 5km, and a 95% chance that they die
# to 0.5 in less than 100km.
#
# dhalf ~ dlnorm(3.107,0.9106)
#
# This prior is derived by assuming that there is a 5% chance that the
# correlations die to 0.5 in less that 5km, and a 95% chance that they die
# to 0.5 in less than 20km.
#
dhalf ~ dlnorm(2.303,0.4214)
phi <- 0.6931/dhalf
beta0 ~ dflat()
beta1 ~ dflat()
# Parameters of interest
base <- exp(beta0)
RRx <- exp(beta1)
vratio <- sigma.U*sigma.U/(sigma.U*sigma.U+sigma.V*sigma.V)
}
DATA 10
list(N = 56,
Y = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,
13, 5, 3, 8, 17, 9, 2, 7, 9, 7,
16, 31, 11, 7, 19, 15, 7, 10, 16, 11,
5, 3, 7, 8, 11, 9, 11, 8, 6, 4,
10, 8, 2, 6, 19, 3, 2, 3, 28, 6,
1, 1, 1, 1, 0, 0),
E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,
4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,
10.5,22.7, 8.8, 5.6,15.5,12.5, 6.0, 9.0,14.4,10.2,
4.8, 2.9, 7.0, 8.5,12.3,10.1,12.7, 9.4, 7.2, 5.3,
18.8,15.8, 4.3,14.6,50.7, 8.2, 5.6, 9.3,88.7,19.6,
3.4, 3.6, 5.7, 7.0, 4.2, 1.8),
X = c(0.16,0.16,0.10,0.24,0.10,0.24,0.10, 0.07, 0.07,0.16,
0.07,0.16,0.10,0.24, 0.07,0.16,0.10, 0.07, 0.07,0.10,
0.07,0.16,0.10, 0.07, 0.01, 0.01, 0.07, 0.07,0.10,0.10,
0.07,0.24,0.10, 0.07, 0.07, 0,0.10, 0.01,0.16, 0,
0.01,0.16,0.16, 0, 0.01, 0.07, 0.01, 0.01, 0, 0.01,
0.01, 0, 0.01, 0.01,0.16,0.10),
xm = c(
162.1894, 385.7761, 293.9555, 377.9338, 220.6786,
340.1739, 324.9915, 442.2445, 194.5176, 367.6924,
112.8916, 247.7566, 289.5922, 227.9563, 342.3574,
351.3505, 280.4916, 341.6081, 249.6855, 359.5902,
348.7138, 388.7655, 180.4228, 295.4908, 333.1159,
312.0605, 290.1701, 359.4153, 291.3727, 303.4219,
257.4402, 264.9711, 336.4464, 258.0319, 227.1801,
234.5294, 218.3428, 279.1010, 235.0805, 254.1736,
250.8301, 287.1202, 292.3773, 288.0333, 320.5682,
257.8758, 276.9737, 281.9644, 267.8444, 342.226,
274.8713, 257.8069, 265.5934, 267.8921, 321.4991,
322.1780),
ym =c(
834.7496, 852.3782, 946.0722, 650.501,
870.9356, 1015.154, 842.0317, 1168.904, 781.3746,
828.219, 903.1592, 924.9536, 842.3052, 561.1628,
713.0808, 792.1617, 801.0356, 628.6406, 825.8545,
610.6554, 760.2982, 812.7655, 699.6693, 635.7658,
701.8189, 691.102, 586.6673, 669.4746, 746.2605,
670.1395, 605.9585, 568.3428, 658.671, 716.452,
598.2521, 668.0481, 641.4785, 670.285, 697.044,
677.589, 657.4675, 680.7535, 699.3761, 665.2905,
671.6064, 631.046, 640.8285, 654.6629, 666.7073,
736.4561, 678.8585, 683.7104, 646.5754, 682.2943,
640.1429, 589.9408))
INIT 10
list(tau.T = 1, a = 1, b = 1,beta0 = 0, beta1 = 0, dhalf =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),
p=c(.5,.5,.5,.5,.5,.5,.5,.5,.5,.5,
.5,.5,.5,.5,.5,.5,.5,.5,.5,.5,
.5,.5,.5,.5,.5,.5,.5,.5,.5,.5,
.5,.5,.5,.5,.5,.5,.5,.5,.5,.5,
.5,.5,.5,.5,.5,.5,.5,.5,.5,.5,
.5,.5,.5,.5,.5,.5),
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))
MODEL 11: Spatial exp correlation model + non-spatial + indiv + individual model
dependent variance priors. New data!
model
{
for (i in 1 : N) {
for (k in 1:K){
Ymat[i,k] ~ dpois(mu[i,k])
log(mu[i,k]) <- log( Nmat[i,k]*(1-x[i])*exp(alpha0+gamma[k]) + Nmat[i,k]*x[i]*exp(alpha0+alpha1+gamma[k]) ) + V[i] + U[i]
}
V[i] ~ dnorm(0,tau.V)
mean[i] <- 0
}
gamma[1] ~ dnorm(0,1.0E10)
gamma[2] ~ dflat()
gamma[3] ~ dflat()
gamma[4] ~ dflat()
gamma[5] ~ dflat()
gamma[6] ~ dflat()
gamma[7] ~ dflat()
gamma[8] ~ dflat()
gamma[9] ~ dflat()
gamma[10] ~ dflat()
# 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)
# tau.T ~ dgamma(1,0.1399)
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 5km, and a 95% chance that they die
# to 0.5 in less than 100km.
#
dhalf ~ dlnorm(3.107,0.9106)
#
# This prior is derived by assuming that there is a 5% chance that the
# correlations die to 0.5 in less that 5km, and a 95% chance that they die
# to 0.5 in less than 20km.
#
# dhalf ~ dlnorm(2.303,0.4214)
phi <- 0.6931/dhalf
alpha0 ~ dflat()
alpha1 ~ dflat()
# Parameters of interest
base <- exp(alpha0)
RRx <- exp(alpha1)
vratio <- sigma.U*sigma.U/(sigma.U*sigma.U+sigma.V*sigma.V)
}
DATA 11
list(N = 56, K = 10,
Nmat = structure(
.Data = c(17666,1512,1473,1545,1518,1618,1429,942,431,190,
156965,13303,12594,11680,10486,10195,7676,4609,2425,1404,
55384,5493,5053,4487,3711,3704,2688,1501,714,455,
31571,2994,3087,3082,2929,2933,2418,1553,692,451,
91614,6981,6554,6112,5197,4946,3761,2364,1182,560,
33425,3058,3123,3033,2909,2965,2315,1361,656,354,
173251,13414,12857,11476,9997,9324,7415,4541,2226,1012,
43768,3157,3197,2698,2550,2658,2142,1420,664,349,
40991,3631,3318,2829,2596,2398,1695,930,500,295,
112083,9211,8843,7916,7103,7101,5820,4087,2195,1195,
54425,4757,5215,4713,4511,4937,4262,2688,1399,908,
23308,2356,2222,1919,1940,2025,1696,1158,621,276,
20319,1676,1532,1193,1192,1263,1090,637,325,147,
56513,5269,5218,4798,4345,4242,3198,1739,787,335,
121349,10042,10087,9829,9457,9169,7382,4710,2251,1196,
74788,6177,5952,5448,5013,5173,4235,2699,1390,790,
17869,1624,1589,1378,1311,1242,975,659,299,129,
59348,5115,5719,5713,5205,5108,3796,2330,1166,645,
112966,9422,9162,8213,7066,6320,4602,2728,1486,902,
65332,5965,6127,5936,5456,5359,4131,2511,1240,640,
174559,14556,14354,14196,12821,12644,9861,5897,2848,1469,
382834,34545,34642,32745,29023,27798,20898,12366,5490,2986,
123830,10271,10104,9618,8942,9720,8062,5551,3092,1626,
112137,9578,9348,8578,7480,6846,4890,2851,1436,674,
290197,25319,26133,23837,20300,18798,13951,8228,3438,1931,
262344,20099,21301,20542,17585,15853,10985,6147,2682,1408,
108450,10086,10065,8905,7951,7421,5450,3190,1479,706,
151060,13546,13936,13610,11847,10895,8176,4760,2215,1140,
225296,19486,19967,19122,17453,16950,13258,8232,4157,2120,
282017,21015,19398,17406,14261,12175,8496,4721,2126,1087,
93437,8080,8521,8295,6922,6137,4105,2220,921,510,
41605,3859,3850,3725,3276,3328,2878,1719,781,427,
178512,14432,14075,12575,10242,8750,5950,3196,1372,563,
157933,13276,12964,12344,10608,10070,7903,4704,2175,1148,
208221,18920,19747,18648,16223,15139,11448,6575,2955,1440,
198969,17212,18054,17535,15077,12699,8813,4747,2152,980,
275294,22120,20824,19281,16794,14973,10919,6653,3002,1653,
226809,17854,17382,16394,13822,11860,7983,4342,1722,904,
164067,11977,12246,11578,10257,8984,6145,3479,1748,796,
105313,9008,9613,9459,8141,6940,4609,2392,980,469,
431113,37589,35792,32415,27161,23128,16074,8611,3707,1823,
282616,25667,25821,24203,20620,19102,14521,8450,3643,1876,
99014,7990,8036,7534,6235,5331,3793,2043,933,385,
307496,26035,26581,24817,21147,18809,13071,7048,2911,1316,
836636,74629,78320,74158,66654,63778,47989,27213,12279,5905,
160344,14338,13913,13408,11762,10611,7340,4037,1583,834,
181522,16344,15529,12144,8170,5995,3640,2055,941,509,
219819,18507,17545,15997,13279,11557,8297,4426,1861,815,
1497228,130850,143185,141383,128687,118300,85137,45598,18422,7563,
365081,32852,33336,30116,26322,24656,17990,9981,4552,2130,
140722,10223,8478,6543,4599,3797,2617,1357,571,287,
74560,7821,7241,6114,4968,4323,2858,1608,735,479,
94718,9357,9292,8348,7349,6937,4866,2940,1411,894,
175199,16024,14144,12101,9690,8426,5892,3244,1362,662,
24386,2227,2200,2141,2000,2129,1788,1093,500,240,
66573,6319,6293,6061,5651,4889,3683,2301,1099,543),
.Dim=c(56,10)),
Ymat = structure(
.Data = c(0,0,0,1,2,1,1,2,1,1,
1,1,1,0,5,7,6,7,6,5,
0,1,0,2,2,2,2,0,2,0,
0,0,1,1,0,2,3,2,0,0,
1,2,0,1,1,2,2,4,1,1,
0,0,0,0,1,1,4,1,0,1,
1,0,1,1,4,3,7,4,2,3,
0,0,0,0,1,1,2,2,0,1,
0,0,0,1,2,1,0,1,0,1,
1,0,2,0,2,5,2,3,5,0,
0,0,1,0,1,3,4,2,1,1,
0,2,0,0,0,1,1,0,0,1,
0,0,0,0,1,1,0,1,0,0,
0,0,0,1,0,1,2,2,1,1,
0,0,0,0,3,4,6,2,2,0,
0,0,0,1,3,1,1,1,1,1,
0,0,0,0,1,0,1,0,0,0,
0,1,0,1,0,0,1,0,0,4,
2,1,1,0,1,3,0,1,0,0,
0,0,0,0,1,2,1,3,0,0,
0,1,2,1,1,3,4,2,2,0,
0,0,0,3,6,7,5,4,3,3,
1,0,1,4,1,2,1,0,0,1,
0,1,1,1,1,1,0,1,0,1,
2,0,1,2,4,3,2,3,1,1,
0,0,1,2,3,2,3,0,1,3,
0,0,0,2,0,1,2,1,0,1,
0,0,0,1,1,1,0,5,0,2,
0,1,1,0,0,3,4,2,3,2,
1,1,0,0,3,3,0,1,2,0,
0,0,1,0,0,3,0,1,0,0,
0,0,0,0,1,0,1,1,0,0,
1,2,0,0,0,0,0,1,3,0,
0,0,0,4,0,0,1,2,0,1,
1,1,1,1,0,5,1,0,1,0,
1,0,1,0,2,2,1,0,0,2,
0,1,1,2,2,1,2,1,1,0,
0,2,2,0,1,1,1,1,0,0,
0,0,2,0,1,1,0,1,0,1,
0,1,0,2,0,1,0,0,0,0,
0,1,1,1,0,2,1,2,0,2,
0,0,1,1,2,0,1,2,1,0,
0,1,0,1,0,0,0,0,0,0,
0,0,0,0,0,3,1,1,1,0,
2,0,2,1,1,4,5,3,0,1,
0,0,0,1,0,0,0,0,2,0,
0,0,0,1,0,0,1,0,0,0,
0,0,0,0,1,0,1,0,0,1,
1,0,2,1,4,6,5,6,2,1,
0,0,2,2,2,0,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,
0,1,0,0,0,0,0,0,0,0,
0,0,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0),
.Dim=c(56,10)),
x = c(0.13,0.15,0.09,0.25,0.09,0.25,0.1,0.07,0.08,0.13,0.05,0.13,0.09,0.2,
0.08,0.13,0.11,0.08,0.04,0.1,0.08,0.01,0.12,0.06,0.01,0.01,0.08,0.06,0.09,
0.02,0.05,0.23,0.02,0.04,0.04,0,0.02,0.01,0.01,0,0.01,0.01,0.01,0,0.01,0.03,
0.01,0.01,0,0.01,0.01,0,0.01,0.01,0.11,0.14),
xm = c(
162.1894, 385.7761, 293.9555, 377.9338, 220.6786,
340.1739, 324.9915, 442.2445, 194.5176, 367.6924,
112.8916, 247.7566, 289.5922, 227.9563, 342.3574,
351.3505, 280.4916, 341.6081, 249.6855, 359.5902,
348.7138, 388.7655, 180.4228, 295.4908, 333.1159,
312.0605, 290.1701, 359.4153, 291.3727, 303.4219,
257.4402, 264.9711, 336.4464, 258.0319, 227.1801,
234.5294, 218.3428, 279.1010, 235.0805, 254.1736,
250.8301, 287.1202, 292.3773, 288.0333, 320.5682,
257.8758, 276.9737, 281.9644, 267.8444, 342.226,
274.8713, 257.8069, 265.5934, 267.8921, 321.4991,
322.1780),
ym =c(
834.7496, 852.3782, 946.0722, 650.501,
870.9356, 1015.154, 842.0317, 1168.904, 781.3746,
828.219, 903.1592, 924.9536, 842.3052, 561.1628,
713.0808, 792.1617, 801.0356, 628.6406, 825.8545,
610.6554, 760.2982, 812.7655, 699.6693, 635.7658,
701.8189, 691.102, 586.6673, 669.4746, 746.2605,
670.1395, 605.9585, 568.3428, 658.671, 716.452,
598.2521, 668.0481, 641.4785, 670.285, 697.044,
677.589, 657.4675, 680.7535, 699.3761, 665.2905,
671.6064, 631.046, 640.8285, 654.6629, 666.7073,
736.4561, 678.8585, 683.7104, 646.5754, 682.2943,
640.1429, 589.9408))
INIT 11
list(tau.T = 1, p=0.5,alpha0 = -13.5, alpha1 = 2.8, dhalf =100,
gamma=c(0,2.8,3.1,3.5,4.1,4.5,4.7,5.1,5.3,6.0),
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),
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))
MODEL 12 is non-spatial and spatial ICAR with individual model and tau.T for islands
model
{
for (i in 1 : N1) {
for (k in 1:K){
Ymat[i1[i],k] ~ dpois(mu[i1[i],k])
log(mu[i1[i],k]) <- log( Nmat[i1[i],k]*(1-x[i1[i]])*exp(alpha0+gamma[k]) + Nmat[i1[i],k]*x[i1[i]]*exp(alpha0+alpha1+gamma[k]) ) + V[i1[i]] + U[i1[i]]
}
V[i1[i]] ~ dnorm(0,tau.V)
}
for (i in 1 : N2) {
for (k in 1:K){
Ymat[i2[i],k] ~ dpois(mu[i2[i],k])
log(mu[i2[i],k]) <- log( Nmat[i2[i],k]*(1-x[i2[i]])*exp(alpha0+gamma[k]) + Nmat[i2[i],k]*x[i2[i]]*exp(alpha0+alpha1+gamma[k]) ) + V[i2[i]] + U[i2[i]]
}
V[i2[i]] ~ dnorm(0,tau.T)
}
gamma[1] ~ dnorm(0,1.0E10)
gamma[2] ~ dflat()
gamma[3] ~ dflat()
gamma[4] ~ dflat()
gamma[5] ~ dflat()
gamma[6] ~ dflat()
gamma[7] ~ dflat()
gamma[8] ~ dflat()
gamma[9] ~ dflat()
gamma[10] ~ dflat()
# ICAR prior distribution for spatial random effects:
U[1:N] ~ car.normal(adj[], weights[], num[], tauomega.U)
for(k in 1:sumNumNeigh) {
weights[k] <- 1
}
#
# We simulate for the total variance and then transform to marginal
# spatial and non-spatial, and then transform to omega
#
tau.T ~ dgamma(1,0.0260)
#tau.T ~ dgamma(1,0.1399)
p ~ dbeta(1,1)
sigma.Z <- sqrt(p/tau.T)
omega.U <- sigma.Z/sqrt(1.164)
sigma.V <- sqrt((1-p)/tau.T)
tau.V <- 1/(sigma.V*sigma.V)
tauomega.U <- 1/(omega.U*omega.U)
alpha0 ~ dflat()
alpha1 ~ dflat()
# Parameters of interest
sd.U <- sd(U[1:N])
vratio <- sd.U*sd.U/(sd.U*sd.U+sigma.V*sigma.V)
RRx <- exp(alpha1)
}
DATA 12
list(N = 56, N1 = 53, N2 = 3,K = 10,
Nmat = structure(
.Data = c(17666,1512,1473,1545,1518,1618,1429,942,431,190,
156965,13303,12594,11680,10486,10195,7676,4609,2425,1404,
55384,5493,5053,4487,3711,3704,2688,1501,714,455,
31571,2994,3087,3082,2929,2933,2418,1553,692,451,
91614,6981,6554,6112,5197,4946,3761,2364,1182,560,
33425,3058,3123,3033,2909,2965,2315,1361,656,354,
173251,13414,12857,11476,9997,9324,7415,4541,2226,1012,
43768,3157,3197,2698,2550,2658,2142,1420,664,349,
40991,3631,3318,2829,2596,2398,1695,930,500,295,
112083,9211,8843,7916,7103,7101,5820,4087,2195,1195,
54425,4757,5215,4713,4511,4937,4262,2688,1399,908,
23308,2356,2222,1919,1940,2025,1696,1158,621,276,
20319,1676,1532,1193,1192,1263,1090,637,325,147,
56513,5269,5218,4798,4345,4242,3198,1739,787,335,
121349,10042,10087,9829,9457,9169,7382,4710,2251,1196,
74788,6177,5952,5448,5013,5173,4235,2699,1390,790,
17869,1624,1589,1378,1311,1242,975,659,299,129,
59348,5115,5719,5713,5205,5108,3796,2330,1166,645,
112966,9422,9162,8213,7066,6320,4602,2728,1486,902,
65332,5965,6127,5936,5456,5359,4131,2511,1240,640,
174559,14556,14354,14196,12821,12644,9861,5897,2848,1469,
382834,34545,34642,32745,29023,27798,20898,12366,5490,2986,
123830,10271,10104,9618,8942,9720,8062,5551,3092,1626,
112137,9578,9348,8578,7480,6846,4890,2851,1436,674,
290197,25319,26133,23837,20300,18798,13951,8228,3438,1931,
262344,20099,21301,20542,17585,15853,10985,6147,2682,1408,
108450,10086,10065,8905,7951,7421,5450,3190,1479,706,
151060,13546,13936,13610,11847,10895,8176,4760,2215,1140,
225296,19486,19967,19122,17453,16950,13258,8232,4157,2120,
282017,21015,19398,17406,14261,12175,8496,4721,2126,1087,
93437,8080,8521,8295,6922,6137,4105,2220,921,510,
41605,3859,3850,3725,3276,3328,2878,1719,781,427,
178512,14432,14075,12575,10242,8750,5950,3196,1372,563,
157933,13276,12964,12344,10608,10070,7903,4704,2175,1148,
208221,18920,19747,18648,16223,15139,11448,6575,2955,1440,
198969,17212,18054,17535,15077,12699,8813,4747,2152,980,
275294,22120,20824,19281,16794,14973,10919,6653,3002,1653,
226809,17854,17382,16394,13822,11860,7983,4342,1722,904,
164067,11977,12246,11578,10257,8984,6145,3479,1748,796,
105313,9008,9613,9459,8141,6940,4609,2392,980,469,
431113,37589,35792,32415,27161,23128,16074,8611,3707,1823,
282616,25667,25821,24203,20620,19102,14521,8450,3643,1876,
99014,7990,8036,7534,6235,5331,3793,2043,933,385,
307496,26035,26581,24817,21147,18809,13071,7048,2911,1316,
836636,74629,78320,74158,66654,63778,47989,27213,12279,5905,
160344,14338,13913,13408,11762,10611,7340,4037,1583,834,
181522,16344,15529,12144,8170,5995,3640,2055,941,509,
219819,18507,17545,15997,13279,11557,8297,4426,1861,815,
1497228,130850,143185,141383,128687,118300,85137,45598,18422,7563,
365081,32852,33336,30116,26322,24656,17990,9981,4552,2130,
140722,10223,8478,6543,4599,3797,2617,1357,571,287,
74560,7821,7241,6114,4968,4323,2858,1608,735,479,
94718,9357,9292,8348,7349,6937,4866,2940,1411,894,
175199,16024,14144,12101,9690,8426,5892,3244,1362,662,
24386,2227,2200,2141,2000,2129,1788,1093,500,240,
66573,6319,6293,6061,5651,4889,3683,2301,1099,543),
.Dim=c(56,10)),
Ymat = structure(
.Data = c(0,0,0,1,2,1,1,2,1,1,
1,1,1,0,5,7,6,7,6,5,
0,1,0,2,2,2,2,0,2,0,
0,0,1,1,0,2,3,2,0,0,
1,2,0,1,1,2,2,4,1,1,
0,0,0,0,1,1,4,1,0,1,
1,0,1,1,4,3,7,4,2,3,
0,0,0,0,1,1,2,2,0,1,
0,0,0,1,2,1,0,1,0,1,
1,0,2,0,2,5,2,3,5,0,
0,0,1,0,1,3,4,2,1,1,
0,2,0,0,0,1,1,0,0,1,
0,0,0,0,1,1,0,1,0,0,
0,0,0,1,0,1,2,2,1,1,
0,0,0,0,3,4,6,2,2,0,
0,0,0,1,3,1,1,1,1,1,
0,0,0,0,1,0,1,0,0,0,
0,1,0,1,0,0,1,0,0,4,
2,1,1,0,1,3,0,1,0,0,
0,0,0,0,1,2,1,3,0,0,
0,1,2,1,1,3,4,2,2,0,
0,0,0,3,6,7,5,4,3,3,
1,0,1,4,1,2,1,0,0,1,
0,1,1,1,1,1,0,1,0,1,
2,0,1,2,4,3,2,3,1,1,
0,0,1,2,3,2,3,0,1,3,
0,0,0,2,0,1,2,1,0,1,
0,0,0,1,1,1,0,5,0,2,
0,1,1,0,0,3,4,2,3,2,
1,1,0,0,3,3,0,1,2,0,
0,0,1,0,0,3,0,1,0,0,
0,0,0,0,1,0,1,1,0,0,
1,2,0,0,0,0,0,1,3,0,
0,0,0,4,0,0,1,2,0,1,
1,1,1,1,0,5,1,0,1,0,
1,0,1,0,2,2,1,0,0,2,
0,1,1,2,2,1,2,1,1,0,
0,2,2,0,1,1,1,1,0,0,
0,0,2,0,1,1,0,1,0,1,
0,1,0,2,0,1,0,0,0,0,
0,1,1,1,0,2,1,2,0,2,
0,0,1,1,2,0,1,2,1,0,
0,1,0,1,0,0,0,0,0,0,
0,0,0,0,0,3,1,1,1,0,
2,0,2,1,1,4,5,3,0,1,
0,0,0,1,0,0,0,0,2,0,
0,0,0,1,0,0,1,0,0,0,
0,0,0,0,1,0,1,0,0,1,
1,0,2,1,4,6,5,6,2,1,
0,0,2,2,2,0,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,
0,1,0,0,0,0,0,0,0,0,
0,0,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0),
.Dim=c(56,10)),
x = c(0.13,0.15,0.09,0.25,0.09,0.25,0.1,0.07,0.08,0.13,0.05,0.13,0.09,0.2,
0.08,0.13,0.11,0.08,0.04,0.1,0.08,0.01,0.12,0.06,0.01,0.01,0.08,0.06,0.09,
0.02,0.05,0.23,0.02,0.04,0.04,0,0.02,0.01,0.01,0,0.01,0.01,0.01,0,0.01,0.03,
0.01,0.01,0,0.01,0.01,0,0.01,0.01,0.11,0.14),
i1 = c(1,2,3,4,5,7,9,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,
46,47,48,49,50,51,52,53,54,55,56), i2 = c(6,8,11),
num = c(3, 2, 1, 3, 3, 0, 5, 0, 5, 4,
0, 2, 3, 3, 2, 6, 6, 6, 5, 3,
3, 2, 4, 8, 3, 3, 4, 4, 11, 6,
7, 3, 4, 9, 4, 2, 4, 6, 3, 4,
5, 5, 4, 5, 4, 6, 6, 4, 9, 2,
4, 4, 4, 5, 6, 5
),
adj = c(
19, 9, 5,
10, 7,
12,
28, 20, 18,
19, 12, 1,
17, 16, 13, 10, 2,
29, 23, 19, 17, 1,
22, 16, 7, 2,
5, 3,
19, 17, 7,
35, 32, 31,
29, 25,
29, 22, 21, 17, 10, 7,
29, 19, 16, 13, 9, 7,
56, 55, 33, 28, 20, 4,
17, 13, 9, 5, 1,
56, 18, 4,
50, 29, 16,
16, 10,
39, 34, 29, 9,
56, 55, 48, 47, 44, 31, 30, 27,
29, 26, 15,
43, 29, 25,
56, 32, 31, 24,
45, 33, 18, 4,
50, 43, 34, 26, 25, 23, 21, 17, 16, 15, 9,
55, 45, 44, 42, 38, 24,
47, 46, 35, 32, 27, 24, 14,
31, 27, 14,
55, 45, 28, 18,
54, 52, 51, 43, 42, 40, 39, 29, 23,
46, 37, 31, 14,
41, 37,
46, 41, 36, 35,
54, 51, 49, 44, 42, 30,
40, 34, 23,
52, 49, 39, 34,
53, 49, 46, 37, 36,
51, 43, 38, 34, 30,
42, 34, 29, 26,
49, 48, 38, 30, 24,
55, 33, 30, 28,
53, 47, 41, 37, 35, 31,
53, 49, 48, 46, 31, 24,
49, 47, 44, 24,
54, 53, 52, 48, 47, 44, 41, 40, 38,
29, 21,
54, 42, 38, 34,
54, 49, 40, 34,
49, 47, 46, 41,
52, 51, 49, 38, 34,
56, 45, 33, 30, 24, 18,
55, 27, 24, 20, 18
),
sumNumNeigh = 234))
INIT 12
list(tau.T = 1, p=0.5,alpha0 = -13.5, alpha1 = 2.8,
gamma=c(0,2.8,3.1,3.5,4.1,4.5,4.7,5.1,5.3,6.0))