[R] Passing Bablok

[dtholmes]

Hi Benjamin,

I was frustrated by the same thing. I pulled the papers last week and wrote
this little function. Produces the same results as MethVal, a commercial
clinical chemistry package. I have not tested it thoroughly but I hope that
helps you.
Dan Holmes, MD, Vancouver


#########################################################
# 	R code for Passing Bablok Regression with			          #     
#   confidence intervals								  #
#	Ref: J Clin Chem Biochem Vol 26, 1988 pp 783-790	          #
# 	Written by Daniel Holmes,MD					          #
# 	Department of Pathology and Laboratory Medicine		  #
# 	University of British Columbia						  #
# 	St. Paul's Hospital 								  #
# 	1081 Burrard St. 							          #
#	Vancouver, BC									  #
#########################################################

isodd <- function(N){
if (N%%2==0) {
    ans<-FALSE
  	} else if (N%%2==1) {
    ans<-TRUE 
	} else {
	ans<-"Neither!"
  }
return(ans)
}

PB.reg<-function(data){
x<-data[,1]
y<-data[,2]
#then we proceed as usual
lx<-length(x)
l<-choose(lx,2)
S<-matrix(1:l,nrow=1,ncol=l)

for (i in 1:(lx-1)) {
	for (j in (i+1):lx) {
	# This expression generates the natural numbers from 
	# i and j	avoiding the use of another counter
	index<-(j-i)+(i-1)*(lx-i/2)
	S[index]<-(y[i]-y[j])/(x[i]-x[j])
	}
}

S.sort<-sort(S)
S.sort<-subset(S.sort,S.sort!=is.na(S.sort))
N<-length(S.sort)
neg<-length(subset(S.sort,S.sort<0))
K<-floor(neg/2)

if (isodd(N)) {
	b<-S.sort[(N+1)/2+neg/2]
} else {
	b<-sqrt(S.sort[N/2+K]*S.sort[N/2+K+1])
}
a<-median(y-b*x)
#CI of b
C.gamma<-qnorm(0.975)*sqrt(lx*(lx-1)*(2*lx+5)/18)
M1<-round((N-C.gamma)/2)
M2<-N-M1+1
b.lower<-S.sort[M1+K]
b.upper<-S.sort[M2+K]
#CI of a
a.lower<-median(y-b.upper*x)
a.upper<-median(y-b.lower*x)
result<-list(intercept=a,intercept.CI=c(a.lower,a.upper),slope=b,slope.CI=c(b.lower,b.upper))
return(result)
}

#Example application to mock data
x<-seq(0:50)
y<-5*x+rnorm(50,0,10)
data<-as.data.frame(cbind(x,y))
reg<-PB.reg(data)
reg

#Plot the regression
plot(x,y,pch=21,bg="gray",main="Passing Bablok Regression")
abline(reg$intercept,reg$slope,col="blue")
abline(reg$intercept.CI[1],reg$slope.CI[1],col="red",lty=2)
abline(reg$intercept.CI[2],reg$slope.CI[2],col="red",lty=2)
correl<-paste("R = ",round(cor.test(x,y)$estimate,3),sep="")
slope<-paste("Slope = ",round(reg$slope,2),"
[",round(reg$slope.CI[1],2),",",round(reg$slope.CI[2],2),"]",sep="")
intercept<-paste("Intercept = ",round(reg$intercept,2),"
[",round(reg$intercept.CI[1],2),",",round(reg$intercept.CI[2],2),"]",sep="")
text<-c(correl,slope,intercept)
legend("topleft",text,title="Regression Data",inset = .05)



-- 
View this message in context: http://r.789695.n4.nabble.com/Passing-Bablok-tp886121p3262090.html
Sent from the R help mailing list archive at Nabble.com.