R/tsregNW.R

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

tsregNW <-
function(x,y,xout=FALSE,outfun=out,iter=10,varfun=pbvar,
corfun=pbcor,plotit=FALSE,tol=.0001,...){
#
#  Compute Theil-Sen regression estimator
#
#  Use Gauss-Seidel algorithm
#  when there is more than one predictor
#
#
x<-as.matrix(x)
xx<-cbind(x,y)
xx<-elimna(xx)
x<-xx[,1:ncol(x)]
x<-as.matrix(x)
y<-xx[,ncol(x)+1]
temp<-NA
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)
}
if(ncol(x)==1){
temp1<-tsp1reg(x,y)
coef<-temp1$coef
res<-temp1$res
}
if(ncol(x)>1){
for(p in 1:ncol(x)){
temp[p]<-tsp1reg(x[,p],y)$coef[2]
}
res<-y-x%*%temp
alpha<-median(res)
r<-matrix(NA,ncol=ncol(x),nrow=nrow(x))
tempold<-temp
for(it in 1:iter){
for(p in 1:ncol(x)){
r[,p]<-y-x%*%temp-alpha+temp[p]*x[,p]
temp[p]<-tsp1reg(x[,p],r[,p],plotit=FALSE)$coef[2]
}
if(max(abs(temp-tempold))<tol)break
alpha<-median(y-x%*%temp)
tempold<-temp
}
coef<-c(alpha,temp)
res<-y-x%*%temp-alpha
}
yhat<-y-res
stre=NULL
temp=varfun(y)
#if(temp==0)print('Warning: When computing strength of association, measure of variation=0')
e.pow=NULL
if(temp>0){
e.pow<-varfun(yhat)/varfun(y)
if(!is.na(e.pow)){
if(e.pow>=1)e.pow<-corfun(yhat,y)$cor^2
e.pow=as.numeric(e.pow)
stre=sqrt(e.pow)
}}
res=NULL
list(coef=coef,residuals=res,Strength.Assoc=stre,Explanatory.Power=e.pow)
}