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,])



# 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
D1=rep(1,n)
D2=rep(1,n)
D=rep(1,n)
Dp=rep(1,n)
Pdiff=rep(1,n)
lc=rep(1,iterates)
qti=rep(1,iterates)

# Diagnostics
# outer "for" loop finds diagnostics for each i
for(i in 1:n)
{
# inner "for" loop creates "iterates" length vectors of
# L(\beta^r|y_i) and unstandardized \qtilde_{r(i)} values, i fixed
# Note \qtilde_{r} =  1/t
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],X[i,],N[i],y[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
#D1[i]=log(sum(lc*qti))-sum(log(lc)*qti)
#D2[i]=mean(log(lc))-log(sum(lc*qti))
#D[i]=D1[i]+D2[i]
PPi=JJ%*% qti
Dp[i]=sum( (PP-PPi)*(log(PP/(1-PP)) - log(PPi/(1-PPi))) )
#Pdiff[i]=sum( PP-PPi)
}
matrix(c(seq(1:n),Dp),ncol=2)
# cases with the highest Dp values
seq(1:n)[Dp>.3]

# This plot does not appear in the book
plot(seq(1:n),Dp,xlab="index",ylab=expression(italic(D[i]^p)),cex.lab=.85)


# Figure 13.12
i=52
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],X[i,],N[i],y[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
PPi=JJ%*% qti
par(mfrow=c(1,1),mar=c(5, 5, 4, 2) + 0.1)
plot(seq(1:n),PP-PPi,xlab="index",ylab=
expression(paste(Pr,"(","y=1","|",
Y,",",x[j],")",-Pr,"(","y=1","|",
Y[(52)],",",x[j],")" )))