Half-Cauchy population distribution for standard deviations...
model {
for (i in 1:10) {
for (j in offset[i]:(offset[i+1]-1)) {
y[j] ~ dnorm(psi[j], inv.sigma.squared[i])
psi[j] <- D*exp(-CL[i]*time[j]/V[i])/V[i]
}
CL[i] <- exp(theta[i, 1])
V[i] <- exp(theta[i, 2])
theta[i, 1:2] ~ dmnorm(mu.theta[], inv.Omega[,])
sigma[i] <- abs(z[i])/sqrt(gamma[i])
z[i] ~ dnorm(0, inv.B.squared)
gamma[i] ~ dgamma(0.5, 0.5)
inv.sigma.squared[i] <- 1/pow(sigma[i], 2)
}
inv.B.squared <- 1/pow(B, 2)
B ~ dunif(0, 100)
mu.theta[1:2] ~ dmnorm(m[], T[,])
inv.Omega[1:2, 1:2] ~ dwish(R[,], k)
Omega[1:2, 1:2] <- inverse(inv.Omega[,])
}
Inits:
list(
B = 0.1,
z = c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1),
gamma = c(1,1,1,1,1,1,1,1,1,1),
theta = structure(
.Data = c(
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201),
.Dim = c(10, 2)),
mu.theta = c(1.064710737, 2.708050201),
inv.Omega = structure(
.Data = c(
25.33071788, 0.0,
0.0, 25.33071788),
.Dim = c(2, 2)))
Data:
list(
D = 30,
y = c(
1.09, 0.75, 0.53, 0.34, 0.23, 0.02,
2.03, 1.28, 1.2, 1.02, 0.83, 0.28,
1.44, 1.3, 0.95, 0.68, 0.52, 0.06,
1.55, 0.96, 0.8, 0.62, 0.46, 0.08,
1.35, 0.78, 0.5, 0.33, 0.18, 0.02,
1.08, 0.59, 0.37, 0.23, 0.17, 0.0,
1.32, 0.74, 0.46, 0.28, 0.27, 0.03,
0.02, 0.0, 1.63, 1.01, 0.73, 0.55,
0.41, 0.01, 0.06, 0.02, 1.26, 0.73,
0.4, 0.3, 0.21, 0.0, 1.3, 0.7,
0.4, 0.25, 0.14, 0.0),
offset = c(1, 7, 13, 19, 25, 31, 37, 45, 53, 59, 65),
time = c(
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
28.0, 32.0, 2.0, 4.0, 6.0, 8.0,
10.0, 24.0, 28.0, 32.0, 2.0, 4.0,
6.0, 8.0, 10.0, 24.0, 2.0, 4.0,
6.0, 8.0, 10.0, 24.0),
m = c(1.064710737, 2.708050201),
T = structure(
.Data = c(
1.0E-4, 0.0,
0.0, 1.0E-4),
.Dim = c(2, 2)),
R = structure(
.Data = c(
0.08, 0.0,
0.0, 0.08),
.Dim = c(2, 2)),
k = 2.0)
node mean sd MC error 2.5% median 97.5% start sample
B 0.05963 0.0261 0.001209 0.02445 0.05482 0.1223 4001 10000
Omega[1,1] 0.1456 0.08493 0.001125 0.05295 0.1257 0.3672 4001 10000
Omega[1,2] 0.01737 0.02722 4.092E-4 -0.02687 0.01392 0.08125 4001 10000
Omega[2,1] 0.01737 0.02722 4.092E-4 -0.02687 0.01392 0.08125 4001 10000
Omega[2,2] 0.02803 0.01713 2.426E-4 0.01013 0.0237 0.07106 4001 10000
mu.theta[1] 1.063 0.121 0.001284 0.8211 1.064 1.299 4001 10000
mu.theta[2] 2.681 0.05737 8.368E-4 2.564 2.681 2.794 4001 10000
sigma[1] 0.01746 0.009934 3.385E-4 0.007621 0.01474 0.04441 4001 10000
sigma[2] 0.1729 0.07763 0.00184 0.08827 0.1538 0.3749 4001 10000
sigma[3] 0.09103 0.03785 8.047E-4 0.04667 0.08188 0.188 4001 10000
sigma[4] 0.08641 0.03409 5.681E-4 0.04597 0.07877 0.1723 4001 10000
sigma[5] 0.03198 0.01504 3.929E-4 0.01542 0.02844 0.07016 4001 10000
sigma[6] 0.03808 0.01724 4.756E-4 0.01847 0.03375 0.08201 4001 10000
sigma[7] 0.05226 0.01705 3.636E-4 0.02997 0.04875 0.09471 4001 10000
sigma[8] 0.058 0.01915 4.001E-4 0.03344 0.05394 0.1066 4001 10000
sigma[9] 0.0477 0.02047 4.759E-4 0.02388 0.04285 0.1017 4001 10000
sigma[10] 0.02581 0.01365 4.568E-4 0.0118 0.02225 0.06013 4001 10000
Log-normal population distribution for standard deviations...
model {
for (i in 1:10) {
for (j in offset[i]:(offset[i+1]-1)) {
y[j] ~ dnorm(psi[j], inv.sigma.squared[i])
psi[j] <- D*exp(-CL[i]*time[j]/V[i])/V[i]
}
CL[i] <- exp(theta[i, 1])
V[i] <- exp(theta[i, 2])
theta[i, 1:2] ~ dmnorm(mu.theta[], inv.Omega[,])
log.sigma[i] ~ dnorm(mu.sigma, inv.omega.sigma.squared)
log(sigma[i]) <- log.sigma[i]
inv.sigma.squared[i] <- 1/pow(sigma[i], 2)
}
mu.sigma ~ dnorm(0, 0.0001)
med.sigma <- exp(mu.sigma)
omega.sigma ~ dunif(0, 100)
inv.omega.sigma.squared <- 1 / pow(omega.sigma, 2)
mu.theta[1:2] ~ dmnorm(m[], T[,])
inv.Omega[1:2, 1:2] ~ dwish(R[,], k)
Omega[1:2, 1:2] <- inverse(inv.Omega[,])
}
Inits:
list(
log.sigma = c(-2,-2,-2,-2,-2,-2,-2,-2,-2,-2),
mu.sigma = -2,
omega.sigma = 0.1,
theta = structure(
.Data = c(
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201,
1.064710737, 2.708050201),
.Dim = c(10, 2)),
mu.theta = c(1.064710737, 2.708050201),
inv.Omega = structure(
.Data = c(
25.33071788, 0.0,
0.0, 25.33071788),
.Dim = c(2, 2)))
Data:
list(
D = 30,
y = c(
1.09, 0.75, 0.53, 0.34, 0.23, 0.02,
2.03, 1.28, 1.2, 1.02, 0.83, 0.28,
1.44, 1.3, 0.95, 0.68, 0.52, 0.06,
1.55, 0.96, 0.8, 0.62, 0.46, 0.08,
1.35, 0.78, 0.5, 0.33, 0.18, 0.02,
1.08, 0.59, 0.37, 0.23, 0.17, 0.0,
1.32, 0.74, 0.46, 0.28, 0.27, 0.03,
0.02, 0.0, 1.63, 1.01, 0.73, 0.55,
0.41, 0.01, 0.06, 0.02, 1.26, 0.73,
0.4, 0.3, 0.21, 0.0, 1.3, 0.7,
0.4, 0.25, 0.14, 0.0),
offset = c(1, 7, 13, 19, 25, 31, 37, 45, 53, 59, 65),
time = c(
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
2.0, 4.0, 6.0, 8.0, 10.0, 24.0,
28.0, 32.0, 2.0, 4.0, 6.0, 8.0,
10.0, 24.0, 28.0, 32.0, 2.0, 4.0,
6.0, 8.0, 10.0, 24.0, 2.0, 4.0,
6.0, 8.0, 10.0, 24.0),
m = c(1.064710737, 2.708050201),
T = structure(
.Data = c(
1.0E-4, 0.0,
0.0, 1.0E-4),
.Dim = c(2, 2)),
R = structure(
.Data = c(
0.08, 0.0,
0.0, 0.08),
.Dim = c(2, 2)),
k = 2.0)
node mean sd MC error 2.5% median 97.5% start sample
Omega[1,1] 0.1511 0.08737 0.001159 0.05688 0.1285 0.3882 4001 10000
Omega[1,2] 0.01709 0.02689 4.377E-4 -0.02746 0.01413 0.07985 4001 10000
Omega[2,1] 0.01709 0.02689 4.377E-4 -0.02746 0.01413 0.07985 4001 10000
Omega[2,2] 0.02723 0.01609 2.11E-4 0.009814 0.02312 0.06801 4001 10000
med.sigma 0.04825 0.01344 2.822E-4 0.02629 0.04683 0.07868 4001 10000
mu.sigma -3.069 0.2751 0.005989 -3.639 -3.061 -2.542 4001 10000
mu.theta[1] 1.061 0.1251 0.001298 0.8044 1.061 1.31 4001 10000
mu.theta[2] 2.684 0.05565 8.193E-4 2.572 2.684 2.793 4001 10000
omega.sigma 0.7455 0.2769 0.008735 0.3354 0.7023 1.386 4001 10000
sigma[1] 0.02228 0.01237 6.062E-4 0.008849 0.01862 0.05339 4001 10000
sigma[2] 0.1359 0.04667 0.001781 0.07818 0.1262 0.2557 4001 10000
sigma[3] 0.07937 0.02785 8.88E-4 0.04446 0.07327 0.1496 4001 10000
sigma[4] 0.07706 0.0258 7.437E-4 0.04368 0.07139 0.1435 4001 10000
sigma[5] 0.03483 0.0152 5.929E-4 0.01652 0.0311 0.07385 4001 10000
sigma[6] 0.038 0.01494 4.896E-4 0.01941 0.03522 0.07479 4001 10000
sigma[7] 0.05098 0.01556 4.476E-4 0.03083 0.04819 0.08736 4001 10000
sigma[8] 0.0552 0.01558 4.655E-4 0.03328 0.05214 0.09341 4001 10000
sigma[9] 0.04546 0.01596 4.755E-4 0.02441 0.0427 0.08433 4001 10000
sigma[10] 0.02899 0.01374 5.145E-4 0.01243 0.02611 0.06437 4001 10000