I'm using Wilcox's R functions (specifically, regplot) to plot a Theil-Sen estimator with a single predictor.

However, regplot doesn't plot a confidence band.

How do I calculate upper and lower percentile bootstrap confidence limits for each y value so that I can plot these limits?


1 Answer 1


The following R code is copied verbatim from Appendix C, page C-20 of U.S. EPA's Statistical Analysis of Groundwater Monitoring Data at RCRA Facilities; Unified Guidance; March 2009 (EPA 530/R-09-007; https://nepis.epa.gov/Exe/ZyPURL.cgi?Dockey=P10055GQ.TXT). The code was written by Kirk M. Cameron, Ph.D. of MacStat Consulting, Ltd.

# R script for Theil-Sen Confidence band
# Compute bootstrapped confidence band around Theil-Sen trend line
# user inputs: list of x-values, list of y-values, desired confidence level
# Note: replace numbers in parentheses below with specific x and y values
# corresponding to data-specific ordered pairs
# x-values should be numeric values representing sampling dates or events
# y-values should be concentration values corresponding to these dates or events
# Script produces a plot of the Theil-Sen trend line, the confidence band around the trend,
# and an overlay of the actual data values
x= c(89.6,90.1,90.8,91.1,92.1,93.1,94.1,95.6,96.1,96.3)
y= c(56,53,51,55,52,60,62,59,61,63)
conf = .90
elimna= function(m){
# remove any rows of data having missing values
m= as.matrix(m)
ikeep= c(1:nrow(m))
for(i in 1:nrow(m)) if (sum(is.na(m[i,])>=1)) ikeep[i]= 0
elimna= m[ikeep[ikeep>=1],]
theilsen2= function(x,y){
# Compute the Theil-Sen regression estimator
# Do not compute residuals in this version
# Assumes missing pairs already removed
ord= order(x)
xs= x[ord]
ys= y[ord]
vec1= outer(ys,ys,"-")
vec2= outer(xs,xs,"-")
v1= vec1[vec2>0]
v2= vec2[vec2>0]
slope= median(v1/v2)
coef= 0
coef[1]= median(y)-slope*median(x)
coef[2]= slope
nb= 1000
temp= matrix(c(x,y),ncol=2)
temp= elimna(temp) #remove any pairs with missing values
x= temp[,1]
y= temp[,2]
n= length(x)
ord= order(x)
cut= min(x) + (0:100)*(max(x)-min(x))/100 #compute 101 cut pts
t0= theilsen2(x,y) #compute trend line on original data
tmp= matrix(nrow=nb,ncol=101)
for (i in 1:nb) {
idx= sample(ord,n,rep=T)
xboot= x[idx]
yboot= y[idx]
tboot= theilsen2(xboot,yboot)
tmp[i,]= tboot$coef[1] + cut*tboot$coef[2]
lb= 0; ub= 0
for (i in 1:101){
lb[i]= quantile(tmp[,i],c((1-conf)/2))
ub[i]= quantile(tmp[,i],c((1+conf)/2))
tband= list(xcut=cut,lo=lb,hi=ub,ths0=t0)
yt= tband$ths0$coef[1] + tband$ths0$coef[2]*tband$xcut

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.