# TODO: Add comment
# 
# Author: E.Korsching Aug 8, 2013
###############################################################################


## Bland-Altmann (1986) - qualitative Betrachtung für die Übereinstimmun  (Pearson Korrelation r alleine z.T. 'mangelhaft' für Übereinstimmung)

bland.altmann.ek <- function(x1,x2=NULL,full){
	# x1, x2 : two measurement vectors  or  one measurement matrix and all versus all columns
	plotBA <- function(x1,x2,a.cor,a.mittel,a.diff,a.upp95u,a.upplim,a.upp95l,a.mdiff,a.low95u,a.lowlim,a.low95l,full){
		if(full){
			tmp <- max(x1,x2)
			plot(x=x1, y=x2, xlab="measurement 1", ylab="measurement 2", xlim=c(0,tmp), ylim=c(0,tmp))
			abline(a=0, b=1, lty=1)
			mtext(text=paste("r =",round(a.cor,3),sep=" "), side=3, line=1)
			plot(x=a.mittel, y=a.diff, xlab="mean (of pair)", ylab="difference", ylim=c(a.low95l,a.upp95u))
			abline(h=a.upp95u, lty=3)
			abline(h=a.upplim, lty=2)
			abline(h=a.upp95l, lty=3)
			abline(h=a.mdiff, lty=2)
			abline(h=a.low95u, lty=3)
			abline(h=a.lowlim, lty=2)
			abline(h=a.low95l, lty=3)
		}else{
			plot(x=a.mittel, y=a.diff, xlab="", ylab="", ylim=c(a.low95l,a.upp95u), xaxt="n", yaxt="n")
			abline(h=a.upp95u, lty=3)
			abline(h=a.upplim, lty=2)
			abline(h=a.upp95l, lty=3)
			abline(h=a.mdiff, lty=2)
			abline(h=a.low95u, lty=3)
			abline(h=a.lowlim, lty=2)
			abline(h=a.low95l, lty=3)
		}
	}
	calcBA <- function(x1,x2){
		a.cor <- cor(x1,x2)
		a.diff <- x1-x2							# vec
		a.mittel <- (x1+x2)/2					# vec  mittlere Groesse des Messwertpaares
		a.mdiff <- mean(a.diff)
		a.sdiff <- sd(a.diff)
		a.upplim <- a.mdiff + 2*a.sdiff
		a.lowlim <- a.mdiff - 2*a.sdiff
		a.n <- length(a.diff)
		a.tval <- qt(0.025, a.n-1, lower.tail=F)
		a.upp95u <- a.upplim + a.tval * sqrt(a.sdiff^2/a.n)
		a.upp95l <- a.upplim - a.tval * sqrt(a.sdiff^2/a.n)
		a.low95u <- a.lowlim + a.tval * sqrt(a.sdiff^2/a.n)
		a.low95l <- a.lowlim - a.tval * sqrt(a.sdiff^2/a.n)
		return(list(a.cor,a.mittel,a.diff,a.upp95u,a.upplim,a.upp95l,a.mdiff,a.low95u,a.lowlim,a.low95l))
	}
	
	if(is.null(x2)){
		nc <- ncol(x1)
		if(!full){
			par(mfrow=c(nc-1,nc-1), mar=c(0.1,0.1,0.1,0.1))
			for(i in 1:(nc-1)){
				k <- i+1
				for(j in k:nc){
					tmp <- calcBA(x1[ ,i], x1[ ,j])
					plotBA(x1[ ,i], x1[ ,j],tmp[[1]],tmp[[2]],tmp[[3]],tmp[[4]],tmp[[5]],tmp[[6]],tmp[[7]],tmp[[8]],tmp[[9]],tmp[[10]],F)
				}
				if(i!=(nc-1)){ for(l in 1:(k-1)){ plot.new() } }	#advance to next plot slot - diagonal arrangement
			}
		}else{
			par(mfrow=c(1,2), mar=c(6,5,5,3))
			for(i in 1:(nc-1)){
				k <- i+1
				for(j in k:nc){
					tmp <- calcBA(x1[ ,i], x1[ ,j])
					plotBA(x1[ ,i], x1[ ,j],tmp[[1]],tmp[[2]],tmp[[3]],tmp[[4]],tmp[[5]],tmp[[6]],tmp[[7]],tmp[[8]],tmp[[9]],tmp[[10]],T)
				}
			}
		}
	}else{
		tmp <- calcBA(x1,x2)
		par(mfrow=c(1,2))
		plotBA(x1,x2,tmp[[1]],tmp[[2]],tmp[[3]],tmp[[4]],tmp[[5]],tmp[[6]],tmp[[7]],tmp[[8]],tmp[[9]],tmp[[10]],T)
	}
}


#aa<-rnorm(20,mean=10,sd=5)
#bland.altmann.ek(x1=aa,x2=0.95*aa + rnorm(20,mean=0,sd=2),full=F)