R/rqtest.R

rqtest <-
function(x,y,qval=.5,nboot=200,alpha=.05,SEED=TRUE,xout=FALSE,outfun=outpro,...){
#
#   Omnibus test when using  a quantile regression estimator
#
x<-as.matrix(x)
if(xout){
x<-as.matrix(x)
flag<-outfun(x,...)$keep
x<-x[flag,]
y<-y[flag]
}
x<-as.matrix(x)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples. Please wait.")
data<-matrix(sample(length(y),size=length(y)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,rqtest.sub,x,y,qval=qval)
# 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.
p<-ncol(x)
if(p==1)stop("Use qregci when p=1")
n<-length(y)
np<-p+1
bvec<-t(bvec)
semat<-var(bvec[,2:np])
temp<-rqfit(x,y,qval=qval)$coef[2:np]
temp<-as.matrix(temp)
test<-t(temp)%*%solve(semat)%*%temp
test<-test*(n-p)/((n-1)*p)
p.value<-1-pf(test,p,n-p)
# Determine adjusted critical level, if possible.
adjusted.alpha=NULL
b1=NULL
if(n<=60){
if(alpha==.1){
if(p==2){
b1<-0-0.001965
b0<-.2179
}
if(p==3){
b1<-0-.003
b0<-.2814
}
if(p==4){
b1<-0-.0058
b0<-.4478
}
if(p==5){
b1<-0-.00896
b0<-.6373
}
if(p>=6){
b1<-0-.0112
b0<-.7699
}}
if(alpha==.05){
if(p==2){
b1<-0-0.001173
b0<-.1203
}
if(p==3){
b1<-0-.00223
b0<-.184
}
if(p==4){
b1<-0-.00476
b0<-.3356
}
if(p==5){
b1<-0-.0063
b0<-.425
}
if(p==6){
b1<-0-.00858
b0<-.5648
}}
if(alpha==.025){
if(p==2){
b1<-0-0.00056
b0<-.05875
}
if(p==3){
b1<-0-.00149
b0<-.1143
}
if(p==4){
b1<-0-.00396
b0<-.2624
}
if(p==5){
b1<-0-.00474
b0<-.3097
}
if(p==6){
b1<-0-.0064
b0<-.4111
}}
if(alpha==.01){
if(p==2){
b1<-0-0.00055
b0<-.043
}
if(p==3){
b1<-0-.00044
b0<-.0364
}
if(p==4){
b1<-0-.0024
b0<-.1546
}
if(p==5){
b1<-0-.00248
b0<-.159
}
if(p==6){
b1<-0-.00439
b0<-.2734
}}
if(!is.null(b1))adjusted.alpha<-b1*n+b0
adjusted.alpha<-max(alpha,adjusted.alpha)
}
list(test.stat=test,p.value=p.value,adjusted.alpha=adjusted.alpha)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.