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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.