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