rm(list = ls())

#Read posterior samples of beta
# Assumes they are in a p x t matrix
setwd("c:\\E-drive\\Books\\LOGLIN3\\BAYES\\")
load("trauma_samp.Rda")
# establish t (number of iterates)
iterates=length(beta[1,])

# Matrix for plots
par(mfrow=c(2,2))

# determine where predictions are made
ISS=seq(0,75,1)
# create a column of 1s
J=1+ISS-ISS
AGE= 60
RTS=3.34


# TI=0, blunt linear predictor
# the command outer computes xy'
w<- outer(J,beta[1,]) + outer(ISS,beta[2,]) +
outer(RTS*J, beta[3,]) + outer(AGE*J, beta[4,])
#TI=1, linear predictor
ww=  w+  outer(J, beta[5,])  + outer(AGE*J,beta[6,])

# turn linear predictors into probabilities
w=exp(w)/(1+exp(w))
ww=exp(ww)/(1+exp(ww))

# compute posterior means of predictive probabilities
y=rowMeans(w)
yy=rowMeans(ww)

# Plot for Age=60 and RTS=3.34
plot(ISS,y,type="l",xlim=c(0,80),ylab="Pr(death)",
    xlab="ISS",main="a. AGE=60 and RTS=3.34",lty=1)
lines(ISS,yy,type="l",lty=2)
legend("bottomright",c("Blunt","Penetrating"),lty=c(1,2))


AGE=10
RTS=3.34

# TI=0, blunt linear predictor
w<- outer(J,beta[1,]) + outer(ISS,beta[2,]) +
outer(RTS*J, beta[3,]) + outer(AGE*J, beta[4,])
#TI=1, linear predictor
ww=  w+  outer(J, beta[5,])  + outer(AGE*J,beta[6,])

w=exp(w)/(1+exp(w))
ww=exp(ww)/(1+exp(ww))
# compute posterior means of predictive probabilities
y=rowMeans(w)
yy=rowMeans(ww)
plot(ISS,y,type="l",xlim=c(0,80),ylab="Pr(death)",
    xlab="ISS",main="b. AGE=10 and RTS=3.34",lty=1)
lines(ISS,yy,type="l",lty=2)
legend("bottomright",c("Blunt","Penetrating"),lty=c(1,2))


AGE= 60
RTS=5.47

# TI=0, blunt linear predictor
w<- outer(J,beta[1,]) + outer(ISS,beta[2,]) +
outer(RTS*J, beta[3,]) + outer(AGE*J, beta[4,])
#TI=1, linear predictor
ww=  w+  outer(J, beta[5,])  + outer(AGE*J,beta[6,])

w=exp(w)/(1+exp(w))
ww=exp(ww)/(1+exp(ww))
# compute posterior means of predictive probabilities
y=rowMeans(w)
yy=rowMeans(ww)
plot(ISS,y,type="l",xlim=c(0,80),ylab="Pr(death)",
    xlab="ISS",main="c. AGE=60 and RTS=5.47",lty=1)
lines(ISS,yy,type="l",lty=2)
legend("bottomright",c("Blunt","Penetrating"),lty=c(1,2))


AGE= 10
RTS=5.47

# TI=0, blunt linear predictor
w<- outer(J,beta[1,]) + outer(ISS,beta[2,]) +
outer(RTS*J, beta[3,]) + outer(AGE*J, beta[4,])
#TI=1, linear predictor
ww=  w+  outer(J, beta[5,])  + outer(AGE*J,beta[6,])

w=exp(w)/(1+exp(w))
ww=exp(ww)/(1+exp(ww))
# compute posterior means of predictive probabilities
y=rowMeans(w)
yy=rowMeans(ww)
plot(ISS,y,type="l",xlim=c(0,80),ylab="Pr(death)",
    xlab="ISS",main="d. AGE=10 and RTS=5.47",lty=1)
lines(ISS,yy,type="l",lty=2)
legend("topleft",c("Blunt","Penetrating"),lty=c(1,2))