R/qhomtv2.R

qhomtv2 <-
function(x,y,nboot=100,alpha=.05,qval=c(.2,.8),plotit=TRUE,SEED=TRUE){
#
#   Test hypothesis of homoscedasticiy  by
#   computing a confidence interval for beta_1-beta_2, the
#   difference between the slopes of the qval[2] and qval[1]
#   regression slopes, where qval[1] and qval[2] are
#   the quantile regression slopes
#   estimated via the Koenker-Bassett method.
#   So by default, use the .8 quantile slope minus the
#   the .2 quantile slope.
#
print("Note: adjusted confidence intervals are used;")
print("they can differ from p-values")
print("FWE is not controlled")
if(length(qval)!=2)stop("Argument qval should have 2 values")
x<-as.matrix(x)
p<-ncol(x)
xy<-elimna(cbind(x,y))
x<-xy[,1:p]
x<-as.matrix(x)
pp<-p+1
y<-xy[,pp]
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)
bvec1<-apply(data,1,qhomtsub2,x,y,qval[1]) # A p by nboot matrix
bvec2<-apply(data,1,qhomtsub2,x,y,qval[2])
bvec1<-as.matrix(bvec1)
bvec2<-as.matrix(bvec2)
if(p==1){
bvec1<-t(bvec1)
bvec2<-t(bvec2)
}
se<-NA
for(j in 1:p)se[j]<-sqrt(var(bvec1[j,]-bvec2[j,]))
temp1<-rqfit(x,y,qval[1])$coef[2:pp]
temp2<-rqfit(x,y,qval[2])$coef[2:pp]
crit<-qnorm(1-alpha/2)
crit.ad<-NA
dif<-temp2-temp1
regci<-NA
regci1<-dif-crit*se
regci2<-dif+crit*se
sig.level<-2*(1-pnorm(abs(dif)/se))
regci.ad<-NA
if(alpha==.05 && qval[1]==.2 && qval[2]==.8)
crit.ad<-qnorm(0-.09/sqrt(length(y))+.975)
output<-matrix(NA,ncol=7,nrow=p)
dimnames(output)<-list(NULL,c("Est 1","Est 2","Dif","SE",
"ci.lower","ci.upper","p.value"))
output[,1]<-temp1
output[,2]<-temp2
output[,3]<-dif
output[,4]<-se
output[,5]<-regci1
output[,6]<-regci2
output[,7]<-sig.level
ci.ad<-c(dif-crit.ad*se,dif+crit.ad*se)
output
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.