va92 <- read.table("../datasets/va92.orig.dat", header=T)
 va93 <- read.table("../datasets/va93.orig.dat", header=T)
 va92$f <- va92$y/va92$n
 va93$f <- va93$y/va93$n


plot(sort(va92$f))

se92_sqrt(va92$f*(1-va92$f)/va92$n)
error.bar( 1:length(va92$f), va92$f[order(va92$f)], lower=1.96*se92[order(va92$f)], bar.end=F,  ylim=c(0,1), ylab="f", xlab="Rank")

postsim_matrix(0,ncol=10000,nrow=length(va92$n))
for(i in 1:10000)
{
# order and rank do the same thing; i.e. what is the rank of the ith obj in the# vector) simulations are done in the order of the observed proportions
  postsim[,i]_rbeta(length(va92$n),
                    va92$y[order(va92$f)]+1,
                    va92$n[order(va92$f)]-va92$y[order(va92$f)]+1)
}

postrank_apply(postsim,2,rank)
# applies the function "rank" to the columns (simulations) to get the rank of each hospital's mean in that simulation
dim(postrank)

round(postsim[1:6,1:10],3)
postrank[1:6,1:10]

# compute expected and median plus quantiles of ranks of each hospital
# (apply function across rows)
exprank <- apply(postrank, 1, mean)
medrank <- apply(postrank, 1, median)
postquant_apply(postrank,1,quantile,c(0.025,0.975))

dim(postquant)


par(mfrow=c(2,1))
error.bar( 1:length(va92$f), va92$f[order(va92$f)], lower=1.96*se92[order(va92$f)], bar.end=F,  ylim=c(0,1), ylab="sample proportion")
plot(exprank, xlab="Ranks of Sample Proportions")
lines(postquant[1,],col=2)
lines(postquant[2,],col=3)
lines(medrank, col=4)

summary(as.factor(postrank[1,]))
tabulate(postrank[1,])

tabulate(postrank[6,])
summary(as.factor(postrank[6,]))

mean(postrank[1,]<1.5)
mean(postrank[6,]<1.5)

apply(postrank,1,function(x,y){mean(x<y)},y=1.5)

apply(postrank,1,function(x,y){mean(x<y)},y=6.5)

apply(postrank,1,function(x){mean(x==6)})

#==================================================================
error.bar <- 
function(x, y = NULL, lower, upper, incr = T, bar.ends = T, gap = T, add = F,
	horizontal = F, ..., xlab = deparse(substitute(x)), xlim, ylim)
{
	draw.null.warn <- function(draw, gap)
	{
		if(!any(draw)) {
			warning("Not enough room for a gap.")
			draw <- !draw
			gap <- 0
		}
		invisible(list(draw = draw, gap = gap))
	}
	if(missing(x))
		stop("no data for x or y")
	if(missing(y)) {
		if(missing(xlab))
			xlab <- "Index"
		y <- x
		x <- time(x)
	}
	n <- length(x)
	if(length(y) != n)
		stop("length of y must equal the length of x")
	center <- if(horizontal) x else y
	if(missing(lower))
		stop("you must provide lower")
	if(length(lower) > 1 && length(lower) != n)
		stop("length of lower must be 1 or equal to the length of x")
	#if incr=T lower is assumed >=0
	if(incr) lower <- center - abs(lower) else lower <- rep(lower, length
			 = n)
	if(any(lower >= center))
		warning(paste(
			"There are values of 'lower' which are greater or equal to ",
			if(horizontal) "x" else "y"))
	if(missing(upper))
		upper <- 2 * center - lower
	else {
		if(length(upper) > 1 && length(upper) != n)
			stop("length of upper must be 1 or equal to the length of x"
				)
		if(incr)
			upper <- center + upper
		else upper <- rep(upper, length = n)
	}
	if(any(upper <= center))
		warning(paste(
			"There are values of 'upper' which are smaller or\nequal to ",
			if(horizontal) "x" else "y"))
	if(!add)
		if(horizontal) {
			if(missing(ylim))
				plot(x, y, xlim = if(missing(xlim)) range(
						c(lower, upper), na.rm = T)
					 else xlim, xlab = xlab, ...)
			else plot(x, y, xlim = if(missing(xlim)) range(c(lower,
						upper), na.rm = T) else xlim,
					ylim = ylim, xlab = xlab, ...)
		}
		else {
			if(missing(xlim))
				plot(x, y, ylim = if(missing(ylim)) range(
						c(lower, upper), na.rm = T)
					 else ylim, xlab = xlab, ...)
			else plot(x, y, ylim = if(missing(ylim)) range(c(lower,
						upper), na.rm = T) else ylim,
					xlim = xlim, xlab = xlab, ...)
		}
	if(horizontal) {
		if(gap)
			gap <- 0.75 * par("cxy")[1]
		draw <- x - lower > gap
		z <- draw.null.warn(draw, gap)
		draw <- z$draw
		gap <- z$gap
		segments(lower[draw], y[draw], x[draw] - gap, y[draw])
		draw <- upper - x > gap
		z <- draw.null.warn(draw, gap)
		draw <- z$draw
		gap <- z$gap
		segments(x[draw] + gap, y[draw], upper[draw], y[draw])
		if(bar.ends) {
			size.bar <- par("cxy")[2]
			segments(lower, y - size.bar, lower, y + size.bar)
			segments(upper, y - size.bar, upper, y + size.bar)
		}
	}
	else {
		if(gap)
			gap <- 0.75 * par("cxy")[2]
		draw <- upper - y > gap
		z <- draw.null.warn(draw, gap)
		draw <- z$draw
		gap <- z$gap
		segments(x[draw], y[draw] + gap, x[draw], upper[draw])
		draw <- y - lower > gap
		z <- draw.null.warn(draw, gap)
		draw <- z$draw
		gap <- z$gap
		segments(x[draw], y[draw] - gap, x[draw], lower[draw])
		if(bar.ends) {
			size.bar <- par("cxy")[1]
			segments(x - size.bar, upper, x + size.bar, upper)
			segments(x - size.bar, lower, x + size.bar, lower)
		}
	}
}