# TODO: Add comment
# 
# from Bablok W, Passing H 1985, Application of statistical procedures in analytical instrument testing
# J. of Automatic Chemistry, Vol 7, No 2, pp 74-79, 1985
#
# Author: R.Voss -- E.Korsching Nov 21, 2011
###############################################################################

bablok.passing <- function(x, y, confidence=0.95, verbose=F)
{
	#ini
	xname=deparse(substitute(x))
	yname=deparse(substitute(y))
	xy <- data.frame(x=x,y=y)
	
	# remove NAs from x and y
	nax <- is.na(xy[,1])
	nay <- is.na(xy[,2])
	xy <- xy[!nax & !nay,,drop=F]
	
	# sort data 
	xy <- xy[order(xy[,1]),,drop=F]
	
	x <- xy[,1]
	y <- xy[,2]
	n <- length(x)
	df <- n-2
	# compute all slopes bi = (Yi-Yj)/(Xi-Xj)
	# sort slopes in ascending order
	# k = number of slopes < -1
	# b = slope which lies k positions higher than the median(bi)
	# a = median(Y - bX)
	xd <- outer(x, x, function(x,y){ x-y })
	yd <- outer(y, y, function(x,y){ x-y })
	bi <- yd/xd
	bi <- bi[upper.tri(bi)]
	bi <- bi[!is.na(bi)]
	bi.inf <- which(bi==is.infinite(bi))
	n.inf <- round(length(bi.inf)/2)      # suppose half of the infs are at the lower end
	bi <- bi[!is.infinite(bi)]
	bi <- sort(bi)
	k <- which(bi < -1)
	k <- length(k)
	med.bi.pos <- which(bi >= median(bi, na.rm=T))[1]
	
	nb <- med.bi.pos + k + n.inf
#	return( list(k=k, med.bi.pos=med.bi.pos, n=length(bi), nb=nb) )
	b <- bi[nb]
	bmed <- median(bi, na.rm=T)
	ai <- y-b*x
	a <- median(ai)
	se.b <- sqrt(var(bi, na.rm=F, use="everything") / length(bi))
	se.a <- sqrt(var(ai, na.rm=F, use="everything") / length(bi))
	ci.a <- c(a-qt(confidence, df)*se.a, a+qt(confidence, df)*se.a)
	ci.b <- c(b-qt(confidence, df)*se.b, b+qt(confidence, df)*se.b)
	
	res <- rbind( intercept=cbind(coefficient=a, "standard error"=se.a, "lower 95%-CL"=ci.a[1], "upper 95%-CL"=ci.a[2], df=df),
			slope=cbind(coefficient=b, "standard error"=se.b, "lower 95%-CL"=ci.b[1], "upper 95%-CL"=ci.b[2], df=df) )
	dimnames(res)[[1]] <- c("intercept", "slope")
	
	ei.1 <- y-a-b*x						# estimator for the theoretical residuals    ??check
	ei <- sum(ei.1^2)					# square sum
	
	#
	if(verbose){
		cat( "\n", "Bablok W, Passing H, J. of Automatic Chemistry, Vol 7, 1985" )
		cat( "\n", c( paste("Passing Bablok Regression : "), paste("X: ",xname,", Y: ",yname,", n=",length(x),sep="") ) )
		cat( "\n", paste("Coefficients of Y = ", round(a,4), " + ", round(b,4), " * X", sep=""), "\n" )
	}
	
	#
	return(list(res=res,ei=ei))
}

#aa <- bablok.passing(y=m.ercc.bc.o[,2], x=m.ercc.bc.o[,4])		#intercept: 1-1 slope: 2-1
# bablok.passing(y=c(3,5,5.5,6,8), x=c(2,4,4.2,4.5,5), verbose=T)