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.