#***************************************************************#
# R code for the ANOVA example in Section 9.5.                  #
# Diasorin data							#
# Adam Branscum                                                 #
# 5-12-2009                                                     #
#***************************************************************#






##### DATA
low = c(91, 46,  95,  60, 33, 410, 105, 43, 189,1097, 54,178, 114, 137, 233, 101, 25,70,357)
normal = c(370, 267, 99,157, 75,1281, 48, 298, 268, 62,804,430,171,694,404)
high = c(75,   52, 1378,  555,  331,  231,  472,  263,  120,   46,  650,  349,  251,  492,  759,
                96,  627,  171, 1584,   69,  368,  509,  486,  354,  351,  839,   88,  162, 1041,  383,
                234, 1130,  503,  244,  606,  457,  460,  283,  767,  576,  628,  239,  583,  428,  452,
                723,  201,  406,  422,  243)
Diasorin <- c(low,normal,high)
logD <- log(Diasorin)
loghigh = log(high)
loglow = log(low)
lognormal = log(normal)





##### LIBRARIES
library(nortest)






##### ANALYSIS
g1 <- c(rep(1,length(low)),rep(2,length(normal)),rep(3,length(high)))
g <- factor(g1, 1:3, c("Low","Normal","High"))


## Boxplots and qq plot from least squares fit to ANOVA
par(mfrow=c(1,2),pty="s")

boxplot(logD~g,range=0, ylab="Diasorin Score",lwd=1.5)


## Test for equal variance
bartlett.test(logD,g) 


## ANOVA
fit <- aov(logD~g)
summary(fit)

qqnorm(resid(fit), main="Residuals from ANOVA \n Anderson-Darling P=0.22")
qqline(resid(fit))

ad.test(resid(fit))   






## Section 9.5.1.1. Allocation and Diagnosis
library(R2WinBUGS)

n1 <- 19
n2 <- 15
n3 <- 50
x <- c(50,150,250,350,450,550,650,750,850,950,1050)
piL <- 0.40
piN <- 0.20
piH <- 0.40
diasorinanovadata <- list("n1","n2","n3","low","normal","high","x","piL","piN","piH")
parameters <- c("mu1","mu2","mu3","tau","med1","med2","med3","relmed21","relmed31","relmed32","MAD","prob", "lowf",                                 "normalf","highf","P","diff21","diff31","diff32","prob21","prob32","pL","pN","pH") 
inits <- list( list(mu1=4.87, mu2=5.39, mu3=6.4, tau=1, lowf=50, normalf=50, highf=50))
diasorinfit <- bugs(diasorinanovadata, inits, parameters,"DiasorinANOVAModel.txt", 
                 working.directory = "H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", 
                 n.chains=1, n.iter=105000, n.thin=1, n.burnin=5000)
print(diasorinfit,digits=3)
attach.bugs(diasorinfit)




### Figure 9.2.  Predictive distributions on the log scale for the low, normal, and high groups
plot(density(log(lowf),n=1500), type='l', xlab="", main="", ylab="", xlim=c(1,10), ylim=c(0,0.5), lwd=3)
lines(density(log(normalf),n=1500),lty=2,lwd=3)
lines(density(log(highf),n=1500),lty=3,lwd=3)
mtext("log (Diasorin Score)",line=3,side=1,cex=1.2)
mtext("Predictive Density",line=2.5,side=2,cex=1.2)
legend("topleft",c("Low","Normal","High"),lty=1:3,lwd=c(2,2,2))

g1 <- apply(pL,2,mean)
g2 <- apply(pN,2,mean)
g3 <- apply(pH,2,mean)


### Figure 9.3.  Plot classification probabilities (posterior means) versus log diasorin score for the low, normal, and high groups
plot(log(x),g1, type='l', xlab="", main="", ylab="", xlim=c(3,7.5), ylim=c(0,1), lwd=3)
lines(log(x),g2,lty=2,lwd=3)
lines(log(x),g3,lty=3,lwd=3)
mtext("log (Diasorin Score)",line=3,side=1,cex=1.2)
mtext("Probability",line=2.5,side=2,cex=1.2)
legend("topleft",c("Low","Normal","High"),lty=1:3,lwd=c(2,2,2))