Classical measurement error model...
model {
for (i in 1:N) {
for (j in 1:3) {
y[i,j] ~ dnorm(psi[i,j], inv.sigma.squared)
psi[i,j] <- alpha[i] + beta[i]*(t[i,j] - tbar)
+ gamma*(
mu0
[i] - y0bar)
}
alpha[i] ~ dnorm(mu.alpha, inv.omega.alpha.squared)
beta[i] ~ dnorm(mu.beta, inv.omega.beta.squared)
mu0[i] ~ dnorm(mu.eps, inv.omega.eps.squared)
y0[i] ~ dnorm(mu0[i], inv.sigma.squared)
}
inv.sigma.squared <- 1/sigma.squared
inv.omega.alpha.squared <- 1/omega.alpha.squared
inv.omega.beta.squared <- 1/omega.beta.squared
inv.omega.eps.squared <- 1/omega.eps.squared
sigma.squared <- pow(sigma, 2)
omega.alpha.squared <- pow(omega.alpha, 2)
omega.beta.squared <- pow(omega.beta, 2)
omega.eps.squared <- pow(omega.eps, 2)
log(sigma) <- log.sigma
log.sigma ~ dunif(-10, 10)
omega.alpha ~ dunif(0, 100)
omega.beta ~ dunif(0, 100)
omega.eps ~ dunif(0, 100)
mu.alpha ~ dnorm(0, 0.0001)
mu.beta ~ dnorm(0, 0.0001)
mu.eps ~ dnorm(0, 0.0001)
gamma ~ dnorm(0, 0.0001)
y0bar <- mean(y0[])
}
Inits:
list(log.sigma = 0, omega.alpha = 0.1, omega.beta = 0.1,
mu.alpha = 0, mu.beta = 0, gamma = 0
,
mu.eps = 0, omega.eps = 0.1,
alpha = 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,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0),
beta = 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,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0),
mu0 = 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,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0)
)
Data...
node mean sd MC error 2.5% median 97.5% start sample
deviance 1119.0 36.68 1.029 1048.0 1119.0 1193.0 10001 50000
gamma 1.004 0.1597 0.004629 0.7242 0.9923 1.353 10001 50000
mu.alpha 6.134 0.1619 0.00212 5.819 6.135 6.453 10001 50000
mu.beta -1.064 0.1406 0.005251 -1.336 -1.065 -0.789 10001 50000
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
Dbar Dhat pD DIC
y 815.571 719.083 96.488 912.059
y0 303.517 250.644 52.873 356.389
total 1119.090 969.727 149.361 1268.450
Baseline measurement as a direct covariate...
model {
for (i in 1:N) {
for (j in 1:3) {
y[i,j] ~ dnorm(psi[i,j], inv.sigma.squared)
psi[i,j] <- alpha[i] + beta[i]*(t[i,j] - tbar)
+ gamma*(
y0
[i] - y0bar)
}
alpha[i] ~ dnorm(mu.alpha, inv.omega.alpha.squared)
beta[i] ~ dnorm(mu.beta, inv.omega.beta.squared)
mu0[i] ~ dnorm(mu.eps, inv.omega.eps.squared)
y0[i] ~ dnorm(mu0[i], inv.sigma.squared)
}
inv.sigma.squared <- 1/sigma.squared
inv.omega.alpha.squared <- 1/omega.alpha.squared
inv.omega.beta.squared <- 1/omega.beta.squared
inv.omega.eps.squared <- 1/omega.eps.squared
sigma.squared <- pow(sigma, 2)
omega.alpha.squared <- pow(omega.alpha, 2)
omega.beta.squared <- pow(omega.beta, 2)
omega.eps.squared <- pow(omega.eps, 2)
log(sigma) <- log.sigma
log.sigma ~ dunif(-10, 10)
omega.alpha ~ dunif(0, 100)
omega.beta ~ dunif(0, 100)
omega.eps ~ dunif(0, 100)
mu.alpha ~ dnorm(0, 0.0001)
mu.beta ~ dnorm(0, 0.0001)
mu.eps ~ dnorm(0, 0.0001)
gamma ~ dnorm(0, 0.0001)
y0bar <- mean(y0[])
}
Inits:
list(log.sigma = 0, omega.alpha = 0.1, omega.beta = 0.1,
mu.alpha = 0, mu.beta = 0, gamma = 0,
mu.eps = 0, omega.eps = 0.1,
alpha = 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,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0),
beta = 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,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0),
mu0 = 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,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0))
Data...
node mean sd MC error 2.5% median 97.5% start sample
deviance 1115.0 34.94 0.6869 1046.0 1115.0 1184.0 10001 50000
gamma 0.675 0.08526 0.001505 0.5055 0.6752 0.8403 10001 50000
mu.alpha 6.137 0.1518 7.699E-4 5.839 6.136 6.434 10001 50000
mu.beta -1.056 0.1466 0.006404 -1.339 -1.056 -0.7769 10001 50000
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
Dbar Dhat pD DIC
y 814.684 716.752 97.932 912.616
y0 300.487 226.607 73.879 374.366
total 1115.170 943.359 171.811 1286.980