## km.plot

km.plot <- function(km, xla='X Label', yla='Y Label', ma='Title', n.risk.at=F, xaxis.at=0, xaxis.v=0,
	col.p=1, col.d=NA, lty1=1, lty2=2, lwd1=2, lwd2=1, gr.names=LETTERS[1:5], axis.font.size=1, reorder=FALSE, ...){

		if(is.na(col.d[1])){ col.d <- col.p }

	n.group <- ng <- length(km$strata)

	if(n.group>9){ stop(' Too many groups. ')}
	simple <- !n.risk.at
	
	if(!simple){nra <- xaxis.at}

	sk <- summary(km)
	gr <- c(sk$strata) ; col.p <- rep(col.p,ng) ; col.d <- rep(col.d,ng) ; lty1 <- rep(lty1,ng) ; lty2 <- rep(lty2,ng)
	lwd1 <- rep(lwd1,ng) ; lwd2 <- rep(lwd2,ng) ; grn <- rep(gr.names,ng)
	
	out <- as.list(1:ng)
	for(i in 1:ng){
		temp <- matrix(0,ncol=5, nrow=sum(gr==i))
		temp[,1] <- sk$time[gr==i]
		temp[,2] <- sk$surv[gr==i]
		temp[,3] <- sk$lower[gr==i] ; temp[,4] <- sk$upper[gr==i]
		temp[,3][is.na(temp[,3])] <- min(temp[,3],na.rm=T)
		temp[,4][is.na(temp[,4])] <- min(temp[,4],na.rm=T)
		temp[,5] <- sk$n.risk[gr==i]
		out[[i]] <- temp
		}
	
	#layout(matrix(1)) ; 
	par(mar=c(5.1,4.1,4.1,2.1))
	
	if(!simple){
		layout(matrix(c(4,3,1,2),2), c(1,3), c(3,1))
		par(mar=c(0,0,4,2))}

	plot(km, conf.int=F, xaxt='n', yaxt='n', col=col.p, xlab=xla, ylab=yla, lty=lty1, lwd=lwd1, 
		main=ma, xlim=c(0,max(km$time)*1.02), cex.lab=axis.font.size)
	axis(side=2, at=seq(0,1,by=.2), labels=seq(0,1,by=.2), las=1)
	if(simple){axis(side=1, at=xaxis.at, labels=xaxis.v)}

	max.time <- rep(0,ng)
	for(i in 1:ng){
		max.time[i] <- km$time[cumsum(km$strata)[i]]
		lines(c(out[[i]][,1],max.time[i]), c(out[[i]][,3],min(out[[i]][,3])), col=col.d[i], lty=lty2[i], lwd=lwd2[i], type='s')
		lines(c(out[[i]][,1],max.time[i]), c(out[[i]][,4],min(out[[i]][,4])), col=col.d[i], lty=lty2[i], lwd=lwd2[i], type='s')
	}
	
	if(!simple){
	par(mar=c(6,0,0,2))
	plot(0, bty=']',type='n',xlab=xla,ylab='',xaxt='n',yaxt='n',ylim=c(1,ng+1.5),xlim=c(0,max(km$time)*1.02), cex.lab=axis.font.size)
	axis(side=1,at=xaxis.at, labels=xaxis.v)
	
	midpoint <- max(km$time)/2 ; midp <- rep(0,ng)
	for(i in 1:ng){midp[i] <- min(out[[i]][,2][out[[i]][,1]<midpoint])}
	
	ypos <- 1:ng ; if(reorder){ ypos <- rank(midp) }

	gr.id <- rep(1:ng, km$strata)
		nrisk <- matrix(0,ncol=ng,nrow=length(xaxis.at))
	for(i in 1:ng){
		nr <- km$n.risk[gr.id==i] ; ti <- km$time[gr.id==i]
		temp.nr <- c(max(nr),nr,0)
		temp.ti <- c(-1,ti,max(ti))
		for(j in 1:length(xaxis.at)){
			nrisk[j,i] <- temp.nr[sum(temp.ti<=xaxis.at[j])+1]}

	text(xaxis.at, .5+ypos[i], nrisk[,i])
	}

	par(mar=c(6,4,0,0))
	plot(0,0, type='n', bty='n', xlab='',ylab='',xaxt='n', yaxt='n',ylim=c(1,ng+1.5), xlim=c(0,1))
	text(0,ng+1.2, 'Number at risk', pos=4)
	
	for(i in 1:ng){
		text(0,.5+ypos[i], labels=grn[i], pos=4)
		lines(c(.4,.8), rep(0+ypos[i],2), lty=lty1[i],lwd=lwd1[i],col=col.p[i])
	}
	
	par(mar=c(0,4,4,2))
	plot(c(0,1),c(0,1),type='n', xaxt='n', yaxt='n', xlab='', ylab='', bty='n')
	axis(4, at=.5, labels=yla, tick=F, font=1, pos=0, cex.axis=axis.font.size)
}

	out
}

## 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]')
}

## 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
	}

tigol <- function(x){ exp(x) / (1+exp(x)) }



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]])}