R/stsregp1.R

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

stsregp1 <-
function(x,y,sc=pbvar,xout=FALSE,outfun=out,...){
#
# Compute the S-type modification of
# the Theil-Sen regression estimator.
# Only a single predictor is allowed in this version
#
xy=elimna(cbind(x,y))
p=ncol(as.matrix(x))
if(p!=1)stop("Current version is limited to one predictor")
p1=p+1
x=xy[,1:p]
y=xy[,p1]
x<-as.matrix(x)
if(xout){
x<-as.matrix(x)
flag<-outfun(x,plotit=plotit,...)$keep
x<-x[flag,]
y<-y[flag]
x<-as.matrix(x)
}
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<-v1/v2
allvar<-NA
for(i in 1:length(slope))allvar[i]<-sc(y-slope[i]*x,...)
temp<-order(allvar)
coef<-0
coef[2]<-slope[temp[1]]
coef[1]<-median(y)-coef[2]*median(x)
res<-y-coef[2]*x-coef[1]
list(coef=coef,residuals=res)
}