R/tshd.R

In musto101/wilcox_R: What the package does (short line)

tshd <-
function(x,y,HD=TRUE,plotit=FALSE,xlab='X',ylab='Y',OPT=FALSE){
#
# Compute the Theil-Sen regression estimator.
# Only a single predictor is allowed in this version
#
#  HD=TRUE, use Harrell-Davis for slopes
#  HD=FALSE, use usual median
#
#  OPT=TRUE, compute the intercept using median(y)-b_1median(X) 
#  OPT=FALSE compute the intercept using median of y-b_1X
#
#
temp<-matrix(c(x,y),ncol=2)
temp<-elimna(temp)     # Remove any pairs with missing values
x<-temp[,1]
y<-temp[,2]
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]
if(!HD)slope<-median(v1/v2,na.rm=TRUE)
if(HD)slope<-hd(v1/v2,na.rm=TRUE)
res=y-slope*x
if(!OPT)int=hd(res)
if(OPT)int=hd(y,na.rm=TRUE)-slope*hd(x,na.rm=TRUE)
coef=c(int,slope)
if(plotit){
plot(x,y,xlab=xlab,ylab=ylab)
abline(coef)
}
list(coef=coef)
}