Here is my primitive and computationally intensive solution to this problem.
Thank you to Michal (see above) for making me aware of the nice polygon
function. I am sure that there are better solutions but this one works for
me at work and for publication.
##################################################################
#Passing Bablok Regression
#Bootstrapped and directly calculated CIs
#Written by Dan Holmes, Department of Pathology and Lab Medicine
#St. Paul's Hospital, Vancouver, BC
#Not gauranteed to produce correct results
#Use at your own risk ;)
##################################################################
library("boot")
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,alpha,indices){
#The data is selected by the set of indices that is passed
d<-data[indices,]
#The data for the PB.reg function is then taken from the resampled set, d
x<-d[,1]
y<-d[,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
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)
#Determine the slope b and the intercept a
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 by Passing and Bablok's method
C.gamma<-qnorm(1-alpha/2)*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 by Passing and Bablok's method
a.lower<-median(y-b.upper*x)
a.upper<-median(y-b.lower*x)
return(as.vector(c(a,b,a.lower,a.upper,b.lower,b.upper)))
}
#Determine bootstrapped CIs for the intercept and slope
PB.boot<-function(data,R,n.fit,alpha){
boot.results<-boot(data=data, alpha=alpha, statistic=PB.reg, R=R)
a.ci<-boot.ci(boot.results, type="bca",index=1) # intercept
b.ci<-boot.ci(boot.results, type="bca",index=2) # slope
a.vector<-boot.results$t[,1]
b.vector<-boot.results$t[,2]
x.points<-seq(min(data$x)-0.1*(max(data$x)-min(data$x)),max(data$x)+0.1*(max(data$x)-min(data$x)),length.out=n.fit)
y.points<-array(dim=R)
reg.CI.data<-matrix(data="NA",nrow=n.fit,ncol=3)
for (i in 1:n.fit){
for (j in 1:R){
#Determine the intersections of all bootstrapped regression lines with
vertical line x=a
y.points[j]<-a.vector[j]+b.vector[j]*x.points[i]
}
#Determine the central (1-alpha)/2 quantiles
lower.y<-quantile(y.points,probs=alpha/2)
upper.y<-quantile(y.points,probs=(1-alpha/2))
reg.CI.data[i,1:3]<-c(x.points[i],lower.y,upper.y)
}
return(list(reg.CI.data=reg.CI.data,a.ci=a.ci,b.ci=b.ci))
}
#Draw the regression and the Bland Altman plot
plotaroo<-function(data,xlab,ylab,R,n.fit,alpha){
reg<-PB.reg(data,alpha)
boot<-PB.boot(data,R,n.fit,alpha)
reg.CI.data<-boot$reg.CI.data
a.ci<-boot$a.ci
b.ci<-boot$b.ci
x<-data[,1]
y<-data[,2]
par(mfrow=c(1,2))
plot(x,y,main="Passing Bablok Regression",xlab=xlab,ylab=ylab)
#Colour in the confidence band with polygons.
for (i in 1:n.fit){
if (i<n.fit) {
vx<-c(reg.CI.data[i,1],reg.CI.data[i,1],reg.CI.data[i+1,1],reg.CI.data[i+1,1])
vy<-c(reg.CI.data[i,2],reg.CI.data[i,3],reg.CI.data[i+1,3],reg.CI.data[i+1,2])
polygon(vx,vy,fillOddEven=T,col="gray",border="gray")
}
}
#Replot the points 'cause the nasty polygons erased your points
points(x,y,pch=21,bg="gray")
#Plot the regression plot
abline(reg[1],reg[2],col="blue")
#Put regression info in the top left hand corner
correl<-paste("R = ",round(cor.test(x,y)$estimate,3),sep="")
slope<-paste("Slope = ",round(reg[2],2),"
[",round(reg[5],2),",",round(reg[6],2),"]",sep="")
intercept<-paste("Intercept = ",round(reg[1],2),"
[",round(reg[3],2),",",round(reg[4],2),"]",sep="")
text<-c(correl,slope,intercept)
legend("topleft",text,title="Regression Data",inset = .05)
lines(reg.CI.data[,1],reg.CI.data[,2],col="black",lwd=2,lty=3)
lines(reg.CI.data[,1],reg.CI.data[,3],col="black",lwd=2,lty=3)
#Plot the line of identity
abline(0,1,lty=3)
#Plot the Bland Altman Plot
diff<-(y-x)
avg<-(x+y)/2
plot(avg,diff,pch=21,bg="gray",ylim=c(-max(abs(diff)),max(abs(diff))),main="Bland
Altman Plot",xlab="Average",ylab="Difference")
abline(h=mean(diff),col="blue")
abline(h=c(-2,2)*sd(diff)+mean(diff),col="black",lty=2)
abline(h=0)
par(mfrow=c(1,1))
#Spit out the results in a readable format
reg.list<-list(intercept=reg[1],slope=reg[2],CI.intercept=c(reg[3:4]),CI.slope=c(reg[5:6]),CI.intercept.bootstrap=a.ci$bca[4:5],CI.slope.bootstrap=b.ci$bca[4:5])
return(reg.list)
}
#Example application to mock data
x<-seq(from=1,to=30,length.out=30)
y<-1.15*x-5+rnorm(30,0,0.25*x)
data<-as.data.frame(cbind(x,y))
x11(width=14)
#R is the number of bootstrap resampling events
#n.fit is the number of points with which to build the regressions
confidence band
#alpha is the level of confidence
plotaroo(data,xlab="Rhubarb by LCMSMS (nmol/L)",ylab="Rhubarb by ELISA
(nmol/L)",R=1000,n.fit=50,alpha=0.05)
--
View this message in context:
http://r.789695.n4.nabble.com/Plotting-confidence-bands-around-regression-line-tp2319909p3300704.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.
