###########################
# 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)))