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.
