###########################
# Exact residual analysis
# ("tricking" WinBUGS into computing 22 posteriors simultaneously)
# "Stacks" example (see WinBUGS Examples Vol I)
# In Carlin and Louis (2008), see second part of Example 2.16
#########################
# Preliminary R code to format dataset:
#
#Data: replicate original Y and "leave-one-out" for each replicate
#Y<-c(42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15)
#N<-length(Y)
#Ymat<-matrix(rep(Y,(N+1)),byrow=T, nrow=(N+1),ncol=N)
#for(i in 1:N){Ymat[i,i]<-NA}
# Need to change row and col as WinBUGS and R use different order for arrays
#dput(t(Ymat),"f:\\Ymat.txt")
#
# HERE IS THE WINBUGS CODE:
#
model
{
# Standardise x's and coefficients
for(i in 1:(N+1)){for (j in 1:p) {b[i,j] <- beta[i,j]/sd(x[,j]) }}
for (j in 1:p) {for (i in 1:N) {
z[i,j]<-(x[i,j]-mean(x[,j]))/sd(x[,j])}}
for(i in 1:(N+1)){ b0[i]<-beta0[i]-b[i,1]*mean(x[,1])
-b[i,2]*mean(x[,2])-b[i,3]*mean(x[,3])}
# statistical model:
for (i in 1:N) { # replicated dataset with NA in ith location
for(j in 1:N){ # datapoint in each replicate set
Y[i,j] ~ dnorm(mu[i,j],tau[i])
mu[i,j] <- beta0[i]+beta[i,1]*z[j,1]
+beta[i,2]*z[j,2]+beta[i,3]*z[j,3]}
# exact CPO:
CPO[i]<-sqrt(tau[i])* exp(-tau[i]/2*(Y[(N+1),i]-mu[i,i])
*(Y[(N+1),i]-mu[i,i]))
# exact standardized residuals:
sresid[i] <- (Y[(N+1),i]-mu[i,i])/sigma[i]
outlier[i] <- step(sresid[i]-1.5) + step(-(sresid[i]+1.5))
} # end of i loop
# full data (approximate) model:
for(j in 1:N){
Y[(N+1),j] ~ dnorm(mu[(N+1),j], tau[(N+1)])
mu[(N+1),j]<-beta0[(N+1)] + beta[(N+1),1]*z[j,1]
+beta[(N+1),2]*z[j,2] + beta[(N+1),3]*z[j,3]
# approximate CPO:
CPO_apr[j]<-sqrt(tau[(N+1)])*exp(-tau[(N+1)]/2
*(Y[(N+1),j]-mu[(N+1),j])*(Y[(N+1),j]-mu[(N+1),j]))
# approximate standardized residuals:
sresid_apr[j] <- (Y[(N+1),j]-mu[(N+1),j])/sigma[(N+1)]
outlier_apr[j]<-step(sresid_apr[j]-1.5)
+ step(-(sresid_apr[j]+1.5))
} # end of j loop
# Priors
for(i in 1:(N+1)) { beta0[i] ~ dnorm(0, 0.00001)
sigma[i] ~ dunif(0.01,100) # vague Gelman prior for sigma
tau[i]<- 1/(sigma[i]*sigma[i])
for (j in 1:p) {beta[i,j] ~ dnorm(0, 0.00001)} # coeffs indep
} # end of i loop
} # end of WinBUGS code
# HERE ARE THE INITS (TWO DIFFERENT CHOICES):
list(beta0 =c(10,10,10,10,10,10,10,10,10,10,10,
10,10,10,10,10,10,10,10,10,10,10),
beta=structure(.Data=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), .Dim = c(22,3)),
sigma =c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,
0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1))
list(beta0 =c(-10,-10,-10,-10,-10,-10,-10,-10,-10,-10,-10,
-10,-10,-10,-10,-10,-10,-10,-10,-10,-10,-10),
beta=structure(.Data=c(10,10,10,10,10,10,10,10,10,10,10,10,
10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,
10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,
10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10),
.Dim = c(22, 3)),
sigma =c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
# HERE IS THE DATA (NOTE DIFFERENT ROW FOR EACH LEAVE-ONE-OUT CASE:
list(p = 3, N = 21,
Y = structure(.Data =
c(NA, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, NA, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, NA, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, NA, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, NA, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, NA, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, NA, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, NA, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, NA, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, NA, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, NA, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, NA, 11, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, NA, 12, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, NA, 8, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, NA, 7, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, NA, 8, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, NA, 8, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, NA, 9, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, NA, 15, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, NA, 15,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, NA,
42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15),
.Dim = c(22, 21)),
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)))