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. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.