R/qregci_1.R

Qregci <-
function(x,y,nboot=100,alpha=.05,
qval=.5,q=NULL,SEED=TRUE,pr=TRUE,xout=FALSE,outfun=outpro,...){
#
#  Test the hypothesis that the quantile regression slopes are zero.
#  Can use the .5 quantile regression line only,
#  the .2 and .8 quantile regression lines, or
#  the .2, .5 and .8 quantile regression lines.
#  In the latter two cases, FWE is controlled for alpha=.1, .05, .025 and .01.
#
if(!is.null(q))qval=q
xx<-elimna(cbind(x,y))
np<-ncol(xx)
p<-np-1
y<-xx[,np]
x<-xx[,1:p]
x<-as.matrix(x)
if(xout){
x<-as.matrix(x)
flag<-outfun(x,plotit=FALSE,...)$keep
x<-x[flag,]
y<-y[flag]
}
x<-as.matrix(x)
n<-length(y)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
if(pr)print("Taking bootstrap samples. Please wait.")
data<-matrix(sample(length(y),size=length(y)*nboot,replace=T),nrow=nboot)
# determine critical value.
crit<-NA
if(alpha==.1)crit<-1.645-1.19/sqrt(n)
if(alpha==.05)crit<-1.96-1.37/sqrt(n)
if(alpha==.025)crit<-2.24-1.18/sqrt(n)
if(alpha==.01)crit<-2.58-1.69/sqrt(n)
crit.fwe<-crit
if(length(qval)==2 || p==2){
if(alpha==.1)crit.fwe<-1.98-1.13/sqrt(n)
if(alpha==.05)crit.fwe<-2.37-1.56/sqrt(n)
if(alpha==.025)crit.fwe<-2.60-1.04/sqrt(n)
if(alpha==.01)crit.fwe<-3.02-1.35/sqrt(n)
}
if(length(qval)==3 || p==3){
if(alpha==.1)crit.fwe<-2.145-1.31/sqrt(n)
if(alpha==.05)crit.fwe<-2.49-1.49/sqrt(n)
if(alpha==.025)crit.fwe<-2.86-1.52/sqrt(n)
if(alpha==.01)crit.fwe<-3.42-1.85/sqrt(n)
}
if(is.na(crit.fwe)){
print("Could not determine a critical value")
print("Only alpha=.1, .05, .025 and .01 are allowed")
}
if(p==1){
bvec<-apply(data,1,Qindbt.sub,x,y,q=qval)
estsub<-NA
for(i in 1:length(qval)){
estsub[i]<-Qreg(x,y,q=qval[i])$coef[2]
}
if(is.matrix(bvec))se.val<-sqrt(apply(bvec,1,FUN=var))
if(!is.matrix(bvec))se.val<-sqrt(var(bvec))
test<-abs(estsub)/se.val
ci.mat<-matrix(nrow=length(qval),ncol=3)
dimnames(ci.mat)<-list(NULL,c("Quantile","ci.lower","ci.upper"))
ci.mat[,1]<-qval
ci.mat[,2]<-estsub-crit*se.val
ci.mat[,3]<-estsub+crit*se.val
}
if(p>1){
if(length(qval)>1){
print("With p>1 predictors,only the first qval value is used")
}
bvec<-apply(data,1,regboot,x,y,regfun=Qreg,qval=qval[1])
se.val<-sqrt(apply(bvec,1,FUN=var))
estsub<-Qreg(x,y,q=qval[1])$coef
test<-abs(estsub)/se.val
ci.mat<-matrix(nrow=np,ncol=3)
dimnames(ci.mat)<-list(NULL,c("Predictor","ci.lower","ci.upper"))
ci.mat[,1]<-c(0:p)
ci.mat[,2]<-estsub-crit*se.val
ci.mat[,3]<-estsub+crit*se.val
}
list(n=length(y),test=test,se.val=se.val,crit.val=crit,crit.fwe=crit.fwe,est.values=estsub,ci=ci.mat)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.