R/regse.R

regse <-
function(x,y,xout=FALSE,regfun=tsreg,outfun=outpro,nboot=200,SEED=TRUE,...){
#
#  Estimate the standard errors and 
#  covariance matrix associated with the estimates of
#  the regression parameters based on the estimator indicated by the
#  argument
#  regfun:  default is Theil--Sen. 
#  So the diagonal elements of the matrix returned by this function
#  are the squared standard errors of the intercept estimator, etc.
#
#  Function returns
# param.estimates: the estimate	of the intercept and slopes
# covar: the covariance matrix	associated with the estimator used
# s.e.:	 the standard errors.
#

if(SEED)set.seed(2)
x<-as.matrix(x)
p1<-ncol(x)+1
p<-ncol(x)
xy<-cbind(x,y)
xy<-elimna(xy)
x<-xy[,1:p]
y<-xy[,p1]
nrem=length(y)
estit=regfun(x,y,xout=xout,...)$coef
if(xout){
m<-cbind(x,y)
flag<-outfun(x,plotit=FALSE,...)$keep
m<-m[flag,]
x<-m[,1:p]
y<-m[,p1]
}
x<-as.matrix(x)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
data<-matrix(sample(length(y),size=length(y)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,regboot,x,y,regfun,xout=FALSE,...)
#Leverage points already removed.
# bvec is a p+1 by nboot matrix. The first row
#                     contains the bootstrap intercepts, the second row
#                     contains the bootstrap values for first predictor, etc.
sqe=var(t(bvec))
list(param.estimates=estit,covar=sqe,s.e.=sqrt(diag(sqe)))
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.