Intrinsic CAR model with spatially structured random effects...
model {
for (i in 1:n) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + alpha0 + alpha1*X[i]/10 + S[i]
RR[i] <- exp(alpha0 + alpha1*X[i]/10 + S[i])
}
S[1:n] ~ car.normal(adj[], weights[], num[], inv.omega.squared)
for(k in 1:sumNumNeigh) { # sumNumNeigh = length of adj[],
weights[k] <- 1 # obtained from Map->Adjacency Tool
}
alpha0 ~ dflat()
alpha1 ~ dnorm(0.0, 1.0E-5)
rr.x <- exp(alpha1)
inv.omega.squared ~ dgamma(0.5, 0.0005)
}
Inits:
list(inv.omega.squared = 1, alpha0 = 0, alpha1 = 0,
S=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))
Data:
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(16,16,10,24,10,24,10, 7, 7,16,
7,16,10,24, 7,16,10, 7, 7,10,
7,16,10, 7, 1, 1, 7, 7,10,10,
7,24,10, 7, 7, 0,10, 1,16, 0,
1,16,16, 0, 1, 7, 1, 1, 0, 1,
1, 0, 1, 1,16,10),
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)
Model with spatially structured and unstructured random effects...
model {
for (i in 1:n) {
Y[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + alpha0 + alpha1*X[i]/10
+ S[i] + H[i]
RR[i] <- exp(alpha0 + alpha1*X[i]/10 + S[i] + H[i])
residRR[i] <- exp(S[i] + H[i]) # residual RR
H[i] ~ dnorm(0, inv.omega.sq.h)
X.pred[i] <- alpha1*X[i]/10
lE[i] <- log(E[i])
}
S[1:n] ~ car.normal(adj[], weights[], num[], tau)
for(k in 1:sumNumNeigh) { # sumNumNeigh = length of adj[],
weights[k] <- 1 # obtained from Map->Adjacency Tool
}
alpha0 ~ dflat()
alpha1 ~ dnorm(0.0, 1.0E-5)
rr.x <- exp(alpha1)
tau ~ dgamma(0.5, 0.0005)
inv.omega.sq.h ~ dgamma(0.5, 0.0005)
sdS <- sd(S[])
sdH <- sd(H[])
sdX <- sd(X.pred[])
sdE <- sd(lE[])
sumvar <- sdS*sdS + sdH*sdH + sdX*sdX + sdE*sdE
pS <- sdS*sdS/sumvar
pH <- sdH*sdH/sumvar
pX <- sdX*sdX/sumvar
pE <- sdE*sdE/sumvar
QR80 <- ranked(RR[], 51)/ranked(RR[], 6)
resQR80 <- ranked(residRR[], 51)/ranked(residRR[], 6)
}
Inits:
list(tau = 1,
inv.omega.sq.h = 1,
alpha0 = 0, alpha1 = 0,
S=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),
H = 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)
)
Data:
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(16,16,10,24,10,24,10, 7, 7,16,
7,16,10,24, 7,16,10, 7, 7,10,
7,16,10, 7, 1, 1, 7, 7,10,10,
7,24,10, 7, 7, 0,10, 1,16, 0,
1,16,16, 0, 1, 7, 1, 1, 0, 1,
1, 0, 1, 1,16,10),
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)
node mean sd MC error 2.5% median 97.5% start sample
QR80 7.555 1.297 0.01185 5.412 7.424 10.46 1001 40000
pE 0.6238 0.03673 5.298E-4 0.549 0.6246 0.6931 1001 40000
pH 0.0182 0.03086 0.001449 1.703E-4 0.003441 0.1125 1001 40000
pS 0.2718 0.05497 0.001634 0.1559 0.2733 0.3773 1001 40000
pX 0.08619 0.04071 9.25E-4 0.01754 0.0829 0.1749 1001 40000
resQR80 4.934 0.9109 0.01352 3.481 4.818 7.016 1001 40000
rr.x 1.587 0.1907 0.004348 1.236 1.58 1.986 1001 40000