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