store()
ind <- inddat()
out <- readdat()

which <- 5

ti <- c('Accelerated t-PA vs. Combination Therapy',
	    'SK+SQ Heparin vs. SK+IV Heparin',
	    'Combined SK vs. Accelerated t-PA')[min(which,3)]

##prior <- 'Noninformative prior'
##prior <- 'Skeptical prior'   ## (1/4,4) possible, .025 prob > 2, <1/2
## sd=.3537729
prior   <- 'Very skeptical prior' ## (1/4,4) possible, .025 prob < .75,>1.333
## sd=.146762
fi <- c('or.tpa.combo','or.sk','or.sk.tpa','or.sk.tpa.skeptical',
        'or.sk.tpa.skeptical2')[which]

if(F) {
	ps.slide('priors', type=3)
	dtruncnorm <- 
	function(x, mean = 0, sd = 1, lower = NA, upper = NA) {
#
# density of truncated normal  - taken from BART
#
        k.upper <- if(!is.na(upper)) pnorm((upper - mean)/sd) else 1
        k.lower <- if(!is.na(lower)) pnorm((lower - mean)/sd) else 0
        K <- 1/(k.upper - k.lower)
        y <- K * dnorm(x, mean, sd)
        y[x < lower] <- 0
        y[x > upper] <- 0
        y
}

	x <- seq(.1,3,length=200)
	for(i in 1:2) {
		d <- dtruncnorm(log(x), mean=0, sd=c(.3537729,.146762)[i],
				        lower=-log(4), upper=log(4))
		if(i==1) plot(x, d, xlab='Odds Ratio', ylab='', 
			          ylim=c(0,3), type='l') else
		lines(x, d, lty=3)
	}
	abline(v=1, lty=2, lwd=1)
	dev.off()
}

ps.slide(fi, type=3)

if(F) {
drawdat2 <- xedit(drawdat) #text() use cex=cex, remove cex= from plot(),
#comment out points(), par(), add xlab, posterior mode, remove title,
#get digits from options(), add xlim

#This was done to create drawdat2, which is now in /suserlib/.Data so
#everyone can access it

store(drawdat2)
}

cex <- .75  # was 1.25 for large plot
options(digits=3)
drawdat2(v='or', trace=F, cex=cex, xlab='Odds Ratio', xlim=c(.5,1.5))

or <- out[,'or']
cl <- quantile(or, c(.025,.975))
options(digits=3)
fcl <- format(cl)
xpos <- c(1.01,1.09,.955,.963,.972)[which]
ypos <- c(5.4,5,6.2,6.4,6.4)[which]


text(xpos,ypos,
	 paste('2.5% = ',fcl[1],
		   '\n97.5% = ',fcl[2],
		   '\n\nProb(OR < 1) = ',format(mean(or < 1)),
		   '\nProb(OR < .95) = ',format(mean(or < .95)),
		   '\nProb(OR < .90) = ',format(mean(or < .90)),
		   '\nProb(.90 < OR < 1/.9) = ',format(mean(or > .9 & or < 1/.9)),
		   sep=''),
	 adj=0, cex=cex)

pstamp(paste(ti,prior,sep='   '))
dev.off()