## sens.spec

sens.spec <- function(ttt){
	# Given a 2x2 table, or a 2x2 matrix 
	# This will calculate sensitivity, specificity, 
	# Positive Predicted Value, Negative Predicted Value
	# The entries of the table are
	#   
	#   			Disease Yes  Disease No
	#   Diagnostic Yes         11           21
	#   Diagnostic No          12           22
	#   
	
	x11 <- ttt[1,1] ; x12 <- ttt[1,2] ; x21 <- ttt[2,1] ; x22 <- ttt[2,2]
	rowt1 <- x11 + x12 ; rowt2 <- x21 + x22
	colt1 <- x11 + x21 ; colt2 <- x12 + x22 ; grat <- colt1 + colt2

	Sensitivity <- round(x11/colt1, 4)
	Specificity <- round(x22/colt2, 4)
	PosPredValu <- round(x11/rowt1, 4)
	NegPredValu <- round(x22/rowt2, 4)

	cbind(Sensitivity, Specificity, PosPredValu, NegPredValu)}

## roc0

roc0 <- function(disease,diag.numeric,cut.point,BigIsBad=T,EqualIsBad=T){
	# disease = T/F or 1/0 ; diag.numeric is numeric
	# if BigIsBad, then diagnosis=T when diag.numeric > cut.point
	# if EqualIsBad, then diagnosis=T when diag.numeric == cut.point
	
		if(length(unique(disease)) != 2){ stop(' disease should be 1/0 or TRUE/FALSE. ') }
	dis <- factor( as.logical(disease), levels=c(T,F) )
		fun1 <- function(x,y){x>y} ; fun2 <- function(x,y){x>=y}
		fun <- list(fun1,fun2)[[1+EqualIsBad]]
	if(!BigIsBad){ cut.point <- -cut.point ; diag.numeric <- -diag.numeric }
	dia <- factor( fun(diag.numeric,cut.point) , levels=c(T,F))
	sens.spec(table(dia,dis))
	}

## roc0v

roc0v <- function(disease, diag.numeric, BigIsBad=T, EqualIsBad=T, pl=F, ...){
	cut.pointV <- c(min(diag.numeric)-1,unique(sort(diag.numeric)),max(diag.numeric)+1)
		
	L <- length(cut.pointV) ; if(BigIsBad){ cut.pointV <- rev(cut.pointV) }
	out <- matrix(0, ncol=3, nrow=L)
		for(i in 1:L){
			temp <- roc0(disease, diag.numeric, cut.pointV[i], BigIsBad, EqualIsBad)
			temp[,2] <- 1-temp[,2]
			out[i,] <- c(cut.pointV[i], temp[,1:2])}
	unik <- !duplicated( out[,2]*100 + out[,3] )
	dd <- data.frame(out[unik,]) ; names(dd) <- c('cut point', 'Sensitivity', '1-Specificity')
	
	x <- dd[,3] ; y <- dd[,2]
	area <- rep(0,ii<-nrow(dd)-1)
		for(i in 1:ii){area[i] <- (y[i]+y[i+1])*(x[i+1]-x[i])/2}
	auc <- sum(area)
	
	if(pl){ plot(x,y,xlim=c(0,1),ylim=c(0,1),las=1,xlab='1-Specificity',ylab='Sensitivity',type='b',pch=19, ...) }
	
	list(details=dd,auc=auc)
	}

## show.colors 

show.colors <- function(x=22,y=30){
par(mfrow=c(1,1))
   par(mai=c(.4,.4,.4,.4), oma=c(.2,0,0,.2))

   plot(c(-1,x),c(-1,y), xlab='', ylab='', type='n', xaxt='n', yaxt='n', bty='n')
   for(i in 1:x){for(j in 1:y){
         k <- y*(i-1)+j ; co <- colors()[k]
         rect(i-1, j-1, i, j, col=co, border=1)}}

   text(rep(-.5,y),(1:y)-.5, 1:y, cex=1.2-.016*y)
   text((1:x)-.5, rep(-.5,x), y*(0:(x-1)), cex=1.2-.022*x)
   title('col=colors()[n]')
}

## within.group.id

within.group.id <- wgi <- function(v,g=NULL){

	one.group <- function(v){
		nl <- length(v)
		new <- which(diff(v)!=0)
		gs <- c(new[1],diff(c(new,nl)))
		rep(1:length(gs), gs)
	}
	
	out1 <- one.group(v)
	if(length(g)>0){
		out1 <- NULL
		g <- factor(g)
		n.group <- length(unique(g))
		for(i in 1:n.group){
			vv <- v[g==unique(g)[i]]
			out1 <- c(out1,one.group(vv))
			}
	}
	out2 <- NULL
		a <- c(which(diff(out1)!=0),length(v))
		b <- c(a[1],diff(a))
	for(i in 1:length(b)){out2 <- c(out2,1:b[i])}
	cbind(out1,out2)
	}

## group.sum

group.sum <- gs <- function(vx,gx, r1=10){
    v <- vx[complete.cases(vx,gx)] ; g <- gx[complete.cases(vx,gx)]

	g <- factor(g) ; lg <- levels(g) ; ng <- length(lg)

	n <- g.mean <- g.med <- g.min <- g.max <- g.sd <-  g.25 <- g.75 <- rep(0,ng+1)

	r2 <- r1 ; if( is.integer(v) ){ r2 <- 0 }	
		rr1 <- function(x){round(x,r1)}
		rr2 <- function(x){round(x,r2)}		

		
	for(i in 1:ng){
		vg <- v[g==lg[i]]
		n[i] <- length(vg) ; g.mean[i] <- mean(vg) ; g.sd[i] <- sd(vg)
		g.med[i] <- median(vg) ; g.max[i] <- max(vg) ; g.min[i] <- min(vg)
		g.25[i] <- quantile(vg, .25) ; g.75[i] <- quantile(vg, .75)
	}
	n[ng+1] <- length(v) ; g.mean[ng+1] <- mean(v) ; g.sd[ng+1] <- sd(v)
	g.med[ng+1] <- median(v) ; g.max[ng+1] <- max(v) ; g.min[ng+1] <- min(v)
	g.25[ng+1] <- quantile(v, .25) ; g.75[ng+1] <- quantile(v, .75)
	
	
	
	out <-data.frame(n, rr2(g.min), rr1(g.25), rr1(g.med), rr1(g.75), rr2(g.max), rr1(g.mean), rr1(g.sd), 
		row.names=c(lg, 'All'))
	names(out) <- c('n', 'min', '1Q', 'med', '3Q', 'max', 'mean', 'sd')
	out
	}

odds.ratio2 <- function(mat, alpha=.05, con.cor =T){
	# mat is a 2x2 matirix
	# con.cor = continuity correction
	
	if(sum(dim(mat) == c(2,2))!=2){stop("\n", "Input must be a 2x2 matrix.", "\n")}
	
	cc <- c(0, .5)[con.cor+1]
	
	n11 <- mat[1,1] ; n12 <- mat[1,2] ; n21 <- mat[2,1] ; n22 <- mat[2,2]
	m11 <- n11 + cc ; m12 <- n12 + cc ; m21 <- n21 + cc ; m22 <- n22 + cc
	
	#theta <- (n11*n22)/(n12*n21)
	theta <- (m11*m22)/(m12*m21)  # with continuity correction
	sighat <- sqrt(1/m11 + 1/m12 + 1/m21 + 1/m22)
	
	za <- qnorm(1-alpha/2)
	log.ci <- round(c(log(theta) - za * sighat, log(theta), log(theta) + za * sighat), 5)
	ori.ci <- round(exp(log.ci), 5)
	rev.or <- round((1/ori.ci)[c(3, 2, 1)], 5)
	
	Log.o.r. <- log(theta)
	ase <- sighat
	zz <- Log.o.r./ase
	pp <- 2*(1-pnorm(abs(zz)))
	
	summ <- data.frame(Log.Odds.Ratio=Log.o.r., ASE=ase, Z=zz, P=pp)
	
	o.r. <- data.frame(round(rbind(log.ci, ori.ci, rev.or), 4))
	names(o.r.) <- c('LWR', 'Point', 'UPR')
	
	
	r.s <- apply(mat, 2, sum)
	c.s <- apply(mat, 1, sum)	
	g.s <- sum(r.s)
	
	ans <- data.frame(cbind(rbind(mat, r.s), c(c.s, g.s)), row.names=c("  1  ", "  2  ", "Total"))
	names(ans) <- c(" A ", " B ", "Total")
	
	list(summ, o.r., ans)
	}
	
odds.ratio <- function(mat, alpha=.05, con.cor=T){
	temp <- odds.ratio2(mat, alpha, con.cor)
	list(temp[[1]], temp[[2]])}