attach('/users/feh/projects/misc/gusto/.Data')
options(digits=3)
store()
fit.lrm <- F
bayes <- T

if(fit.lrm) {
  dead <- day30=='dead'
  print(summary(dead ~ tx))

  accel <- tx=='tpa'
  f.tpa <- lrm(dead ~ accel, subset=tx %in% c('tpa','combo'),
			   linear.predictors=F)

  sq <- tx=='sk+sq'
  f.sk  <- lrm(dead ~ sq, subset=tx %in% c('sk+iv','sk+sq'),
			   linear.predictors=F)

  f.accel <- lrm(dead ~ accel, subset=tx!='combo', linear.predictors=F)
  combo <- tx=='combo'
  f.combo <- lrm(dead ~ combo, subset=tx!='tpa', linear.predictors=F)

  print(f.tpa)
  print(f.sk)
  print(f.accel)
  print(f.combo)
  store(f.tpa); store(f.sk); store(f.accel); store(f.combo)

}

if(bayes) {
  library(bart)

  pbart <- function(fit, main) {
	effect <- fit$coef[2]
	sd     <- sqrt(fit$var[2,2])
	par(mfrow=c(2,1))
	z <- paste('Odds ratio=',format(exp(fit$coef[2])),sep='')
	plot(bart(effect, sd), main=paste(main,'\n',z,sep=''), 
		 which=c(1,3), xlim=c(-.75,.75), xlab='Log Odds Ratio')
	invisible()
  }
  postscript('bart.ps',hor=F)
  pbart(f.tpa, 'Accelerated t-PA vs. t-PA Combo')
  pbart(f.sk,  'SK+IV Heparin vs. SK+SQ Heparin')
  pbart(f.accel, 'Accelerated t-PA vs. Both SK')
  pbart(f.combo, 'Combo t-PA vs. Both SK')
  graphics.off()

  ## Analyze accel t-pa vs. combined SK for a variety of truncated
  ## normal prior distributions

  effect <- f.tpa$coef[2]
  sd <- sqrt(f.tpa$var[2,2])
  prior.sds <- c(25,10,5,2.5,1)
  prior.equiv <- p.benefit <- p.equiv <- p.efficacy <- prior.sds
  y.prior <- y.post <- list()
  i <- 0
  for(sd.prior in prior.sds) {
	i <- i+1
	b <- bart(effect, sd, prior1.n=1/(sd.prior^2), prior1.type='gaussian')
	d <- plot(b, pl=F)
	if(i==1) x <- d$x
	prior.equiv[i] <- diff(d$tails.prior)
	p.benefit[i] <- d$tails.post[1]
	p.equiv[i]   <- diff(d$tails.post)
	y.prior[[i]] <- d$y.prior
	y.post[[i]]  <- d$y.post
	d <- plot(bart(effect, sd, prior1.n=1/(sd.prior^2), prior1.type='gaussian',
			  equiv=c(0,0)), pl=F)
	p.efficacy[i] <- d$tails.post[1]
  }
  ps.slide('many.priors',type=1)
  plot(0,0,xlab='Log Odds Ratio',type='n',xlim=c(-.75,.75),ylim=c(0,10),
	   ylab='Prior Density')
  abline(v=log(c(1,1/1.05,1.05)), lty=1, lwd=.5)
  k <- length(y.prior)
  par(err=-1)
  for(i in 1:k) lines(x, y.prior[[i]], lty=i)
  title('Prior Distributions')
  txt <- function(lab, vals) paste(paste(lab,'\n',sep=''),
								   paste(format(vals),collapse='\n'),sep='')
  cx <- .7
  text(.55,8,txt('P{Equivalence}', prior.equiv), cex=cx)

  plot(0,0,xlab='Log Odds Ratio',type='n',xlim=c(-.75,.75),ylim=c(0,10),
	   ylab='Posterior Density')
  abline(v=log(c(1,1/1.05,1.05)), lty=1, lwd=.5)
  for(i in 1:k) lines(x, y.post[[i]], lty=i)
  title('Accelerated t-PA vs. Combined SK')
  text(.55,8, txt('P{Equivalence}',p.equiv), cex=cx)
  text(-.65,8,txt('P{Clinical Benefit}',p.benefit), cex=cx)
  text(-.35,8,txt('P{Efficacy}',p.efficacy), cex=cx)
  graphics.off()

}