###########################
# Approximate residual analysis
# "Stacks" example (see WinBUGS Examples Vol I)
# In Carlin and Louis (2008), see Example 2.16
###########################
# Approximate "leave-one-out" CPOs and standardized residuals:
model
{
## first, standardise x's and coefficients##
for (j in 1 : p) {
b[j] <- beta[j] / sd(x[ , j ])
for (i in 1 : N) {z[i, j] <- (x[i, j] - mean(x[, j])) / sd(x[ , j])}
}
b0 <- beta0 - b[1] * mean(x[,1]) - b[2] * mean(x[,2]) - b[3] * mean(x[,3])
# statistical model
for (i in 1 : N) {
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 +beta[1] * z[i, 1] + beta[2] * z[i, 2] +beta[3] * z[i, 3]
# standardized residuals without the "leave-one-out", and taking the
# the *ratio* before posterior averaging. Results tend to be too close to 0.
sresid_apr[i] <- (Y[i] - mu[i]) / sigma
# flag observations with sresid_apr values bigger than 1.5 as outliers:
outlier_apr[i] <- step(sresid_apr[i] - 1.5) + step(-(sresid_apr[i] +1.5) )
# approximate CPO:
CPO_apr[i] <-sqrt(tau)* exp(-tau/2* (Y[i] - mu[i])* (Y[i] -mu[i]))
}
# priors
beta0 ~ dflat()
for (j in 1 : p) {beta[j] ~ dflat()}
sigma ~ dunif(0.01,100) # vague Gelman prior for sigma
tau<- 1/(sigma*sigma)
} # end of stacks model
#Inits:
list(beta0 = 10, beta=c(0,0, 0), sigma = 0.1)
list(beta0 = -10, beta=c(10,10, 10), sigma = 10)
#Data:
list(p = 3, N = 21,
Y = c(42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15),
x = structure(.Data = c( 80, 27, 89,
80, 27, 88,
75, 25, 90,
62, 24, 87,
62, 22, 87,
62, 23, 87,
62, 24, 93,
62, 24, 93,
58, 23, 87,
58, 18, 80,
58, 18, 89,
58, 17, 88,
58, 18, 82,
58, 19, 93,
50, 18, 89,
50, 18, 86,
50, 19, 72,
50, 19, 79,
50, 20, 80,
56, 20, 82,
70, 20, 91), .Dim = c(21, 3)))