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