Statistical Applications for Environmental Analysis and Risk Assessment

Appendix ASoftware Scripts

A.1 R Script for Confidence Band Around Theil–Sen Line

# Unified Guidance R software package, v1.0
# Trends/nonparametric
# Author -- Kirk Cameron, MacStat Consulting, Ltd
# Last Update -- Jun 26, 2013, rv01
# function to compute bootstrapped ptwise confidence band around Theil-Sen line and
# then to estimate the simultaneous coverage;
# assumes any missing pairs have already been removed from (x,y)
# calls function theilsen.line()
# nb = number of bootstrap replications; conf = confidence level
# side= type of confidence band desired (lo= lower one-sided, hi= upper one-sided, two= two-sided)
SodiumData = read.table("SodiumConc.dat", header = TRUE)
x = SodiumData$Date
y = SodiumData$Sodium
ug.theilsen.band= function(x,y,conf=0.95,nb=5000,side='hi') {
 side= match.arg(side)
 n= length(x); ncut= 101
# compute 101 cut pts along range of x; compute trend estimate on original data
 cut= min(x) + (0:(ncut-1))*(max(x)-min(x))/(ncut-1)
 t0= theilsen.line(x,y)
 tmp= matrix(nrow=nb,ncol=ncut)
# compute nb bootstrap replicate Theilsen lines
 for (i in 1:nb) {
 idx= sample(1:n,n,rep=T)
 xboot= x[idx]; yboot= y[idx]
 tboot= theilsen.line(xboot,yboot)
 tmp[i,]= tboot[1] + cut*tboot[2]
 }
# compute quantiles of boot matrix at each cut; choose quantiles to account for type of
# band desired (hi,lo,two-sided)
 bb= matrix(nrow=ncut,ncol=2)
 if (side=='two') {
 for (i in 1:ncut) bb[i,]= quantile(tmp[,i],p=c((1-conf)/2,(1+conf)/2))
 } else {
 for (i in 1:ncut) bb[i,]= ...