rm(list = ls())

# Enter data
Trauma <- read.table(
url("http://stat.unm.edu/~fletcher/LLM/DATA/TRAUMAa.DAT"),
#"C:\\E-drive\\Books\\LOGLIN3\\DATA\\TRAUMAa.DAT",
		  sep="",col.names=c("ID","Death","ISS","TI","RTS","AGE"))
attach(Trauma)
n=length(ISS)
y=1-Death
N=1+ISS-ISS

setwd("c:\\E-drive\\Books\\LOGLIN3\\BAYES\\")
load("trauma_samp.Rda")

iterates=length(beta[1,])

# Enter prior information
a=c(1.1,3,5.9,1.3,1.1,1.5)
b=c(8.5,11,1.7,12,4.9,5.5)

#Enter Xtilde
Xtp=c(
 1, 25, 7.84, 60, 0,  0,
 1, 25, 3.34, 10, 0,  0,
 1, 41, 3.34, 60, 1, 60,
 1, 41, 7.84, 10, 1, 10,
 1, 33, 5.74, 35, 0,  0,
 1, 33, 5.74, 35, 1, 35)
Xtilde=t(matrix(Xtp,6,6))
Xtilde

# define likelihood contribution function
LC <- function(beta,x,N,y)
{
LP = sum(beta*x)
lc=((exp(LP)/(1+exp(LP)))**y)*((1-(exp(LP)/(1+exp(LP))))**(N-y))
return(lc)}


# Find n vector Pr(y=1|Y,x_i)
# X is model matrix
X=matrix(c(rep(1,n),ISS,RTS,AGE,TI,AGE*TI),nrow=n)
JJ=exp(X %*% beta)
JJ=JJ/(1+JJ)
PP=rowMeans(JJ)

# Define dummy vectors for diagnostics
# to be redefined below
Pdiff=rep(1,n)
lc=rep(1,iterates)
qti=rep(1,iterates)

# Figure 13.14
par(mfrow=c(3,2))#,mar=c(5, 5, 4, 2) + 0.1)
i=1
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],Xtilde[i,],b[i]+a[i],a[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
PPi=JJ%*% qti

plot(seq(1:n),PP-PPi,xlab="index",ylab=
expression(paste(Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y),")",-Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y)[(1)],")" )),main="Omit location 1")

i=2
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],Xtilde[i,],b[i]+a[i],a[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
PPi=JJ%*% qti

plot(seq(1:n),PP-PPi,xlab="index",ylab=
expression(paste(Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y),")",-Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y)[(2)],")" )),main="Omit location 2")


i=3
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],Xtilde[i,],b[i]+a[i],a[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
PPi=JJ%*% qti

plot(seq(1:n),PP-PPi,xlab="index",ylab=
expression(paste(Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y),")",-Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y)[(3)],")" )),main="Omit location 3")

i=4
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],Xtilde[i,],b[i]+a[i],a[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
PPi=JJ%*% qti

plot(seq(1:n),PP-PPi,xlab="index",ylab=
expression(paste(Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y),")",-Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y)[(4)],")" )),main="Omit location 4")

i=5
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],Xtilde[i,],b[i]+a[i],a[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
PPi=JJ%*% qti

plot(seq(1:n),PP-PPi,xlab="index",ylab=
expression(paste(Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y),")",-Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y)[(5)],")" )),main="Omit location 5")

i=6
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],Xtilde[i,],b[i]+a[i],a[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
PPi=JJ%*% qti

plot(seq(1:n),PP-PPi,xlab="index",ylab=
expression(paste(Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y),")",-Pr,"(","y=1","|",
Y,",",x[j],",",tilde(Y)[(6)],")" )),main="Omit location 6")