R/qhomt.R

qhomt <-
function(x,y,nboot=100,alpha=.05,qval=c(.2,.8),plotit=TRUE,SEED=TRUE,
xlab="X",ylab="Y",xout=FALSE,outfun=outpro,pr=TRUE,...){
#
#   Test hypothesis that the error term is  homogeneous 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.
#
if(length(qval)!=2)stop("Argument qval should have 2 values exactly")
x<-as.matrix(x)
if(ncol(x)!=1)stop("Only one predictor is allowed; use qhomtv2")
xy<-elimna(cbind(x,y))
x<-xy[,1]
x<-as.matrix(x)
y<-xy[,2]
if(xout){
flag<-outfun(x,...)$keep
x<-as.matrix(x)
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.
if(pr)print("Taking bootstrap samples. Please wait.")
data<-matrix(sample(length(y),size=length(y)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,qhomtsub,x,y,qval) # An nboot vector.
se<-sqrt(var(bvec))
temp<-qplotreg(x,y,qval=qval,plotit=plotit,xlab=xlab,ylab=ylab)
crit<-qnorm(1-alpha/2)
crit.ad<-NA
dimnames(temp)=c(NULL,NULL)
dif<-temp[2,2]-temp[1,2]
regci<-NA
regci[1]<-dif-crit*se
regci[2]<-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)
ci.ad<-c(dif-crit.ad*se,dif+crit.ad*se)
list(slope.bottom=temp[1,2],slope.top=temp[2,2],
dif.est=dif,dif.ci=regci,sig.level=sig.level,se=se,adjusted.ci=ci.ad)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.