R/tsp1reg.R

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

tsp1reg <-
function(x,y,plotit=FALSE,HD=FALSE,OPT=TRUE){
#
# Compute the Theil-Sen regression estimator.
# Only a single predictor is allowed in this version
#
#  OPT=TRUE, compute the intercept using median(y)-beta_1median(X)
#  OPT=FALSE compute the intercept using median of y-beta_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)
if(OPT){
if(!HD)coef<-median(y,na.rm=TRUE)-slope*median(x,na.rm=TRUE)
if(HD)coef<-hd(y,na.rm=TRUE)-slope*hd(x,na.rm=TRUE)
}
if(!OPT){
if(!HD)coef<-median(y-slope*x,na.rm=TRUE)
if(HD)coef<-hd(y-slope*x,na.rm=TRUE)
}
names(coef)<-"Intercept"
coef<-c(coef,slope)
if(plotit){
plot(x,y,xlab="X",ylab="Y")
abline(coef)
}
res<-y-slope*x-coef[1]
list(coef=coef,residuals=res)
}