## Diasorin as a diagnostic test for low bone turnover in kidney dialysis patients ## Comparing two means ## Ch 5, Section 2 n1 <- 19 n2 <- 15 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) diasorin <- c(low,normal) group <- c(rep(1,n1),rep(2,n2)) #################################### # Normality. Fig 5.2 # #################################### group1 <- factor(group, 1:2, c('Low','Normal')) par(mfrow=c(2,2), pty="s") boxplot(diasorin~group1,range=0, ylab="Diasorin", lwd=1.5) boxplot(log(diasorin)~group1,range=0, ylab="log(Diasorin)", lwd=1.5) qqnorm(log(diasorin[group==1]),main = "Low Group \n Log transformed data", lwd=2) qqline(log(diasorin[group==1]), lwd=2) qqnorm(log(diasorin[group==2]),main = "Normal Group \n Log transformed data", lwd=2) qqline(log(diasorin[group==2]), lwd=2) #################################### # Descriptive Statistics # #################################### library(Hmisc) g1 <- function(x)c(MEAN=mean(x,na.rm=TRUE),SD=sd(x,na.rm=TRUE),MEDIAN=median(x,na.rm=TRUE),MIN=min(x,na.rm=TRUE),MAX=max(x,na.rm=TRUE),N=length(x)) summarize(log(diasorin),by=group1,FUN=g1) Low Normal Log(Diasorin) N 19 15 Mean (SD) 4.71 (0.94) 5.49 (0.97) Median (Min, Max) 4.62 (3.22, 7.00) 5.59 (3.87, 7.16) ############################################################ # Example 5.2.2. # # Use R2WinBUGS for the 2 sample log normal Diasorin data # ############################################################ library(R2WinBUGS) n1 <- 19 n2 <- 15 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) diasorindata <- list("n1","n2","low","normal") parameters <- c("muL","muN","tauL","tauN","medianL","medianN","relmedian", "test", "lowf", "normalf","relprec") inits <- list( list(muL=4.87, muN=5.39, tauL=1, tauN=1, lowf=50, normalf=50)) diasorin.fit <- bugs(diasorindata, inits, parameters,"DiasorinModel.txt", working.directory = "F:\\MyDocuments\\BAYES\\OneandTwoSample\\Code", n.chains=1, n.iter=60000, n.thin=1, n.burnin=10000) parameters1 <- c("muL1","muN1","tauL1","tauN1","medianL1","medianN1","relmedian1", "test1", "lowf1", "normalf1") inits1 <- list( list(muL1=4.87, muN1=5.39, tauL1=1, tauN1=1, lowf1=50, normalf1=50)) diasorin.fit1 <- bugs(diasorindata, inits1, parameters1,"DiasorinModel1.txt", working.directory = "F:\\MyDocuments\\BAYES\\OneandTwoSample\\Code", n.chains=1, n.iter=60000, n.thin=1, n.burnin=10000) print(diasorin.fit,digits=3) attach.bugs(diasorin.fit) round(print(diasorin.fit1),3) attach.bugs(diasorin.fit1) ### Figure 5.4. Predictive distributions on the log scale for the low and normal groups attach.bugs(diasorin.fit) plot(density(log(lowf),n=1500), type='l', xlab="", main="", ylab="", xlim=c(2,8), ylim=c(0,0.5), lwd=5) lines(density(log(normalf),n=1500),lty=2,lwd=5) mtext("log (Diasorin Score)",line=3,side=1,cex=1.5) mtext("Predictive Density",line=2.5,side=2,cex=1.5) text(3.2,0.3,"Low",lwd=2,cex=1.2) text(7,.3,"Normal",lwd=2,cex=1.2) ### Figure 5.5. Predictive distributions for the low and normal groups attach.bugs(diasorin.fit) ## Fully informative priors plot(density(lowf,n=2500,from=0), type='l', xlab="", main="", ylab="", xlim=c(0,1000), ylim=c(0,0.006), lwd=5) lines(density(normalf,n=2500,from=0),lty=2,lwd=5) lines(c(0,0),c(0,0.000593),lwd=5,lty=1) mtext("Diasorin Score",line=3,side=1,cex=1.5) mtext("Predictive Density",line=2.5,side=2,cex=1.5) text(50,0.0055,"Low",lwd=2,cex=1.2) text(380,0.0018,"Normal",lwd=2,cex=1.2) attach.bugs(diasorin.fit1) ## Priors centered at elicited values, but otherwise fat. plot(density(lowf1,n=2500), type='l', xlab="", main="", ylab="", xlim=c(0,600), ylim=c(0,0.006), lwd=5) lines(density(normalf1,n=2500),lty=2,lwd=5) mtext("Diasorin Score",line=3,side=1,cex=1.5) mtext("Predictive Density",line=2.5,side=2,cex=1.5) text(50,0.0055,"Low",lwd=2,cex=1.2) text(350,0.0018,"Normal",lwd=2,cex=1.2) c(mean(log(lowf)), sd(log(lowf))) c(mean(log(normalf)), sd(log(normalf))) Mean and SD of the predictive densities for low and normal Low: mean=4.85 sd=0.93 Normal: mean= 5.39 sd=0.93 SElow <- sd(low)/sqrt(n1) SEnormal <- sd(normal)/sqrt(n2) SElow SEnormal SEloglow <- sd(log(low))/sqrt(n1) SElognormal <- sd(log(normal))/sqrt(n2) SEloglow SElognormal