R/qsm.R

qsm <-
function(x,y,qval=c(.2,.5,.8),fr=.8,plotit=TRUE,scat=TRUE,pyhat=FALSE,eout=FALSE,xout=FALSE,outfun=out,op=TRUE,LP=TRUE,
xlab='X',ylab='Y'){
#
# running  interval smoother for the quantiles stored in
# qval
#
# fr controls amount of smoothing
# op=T, use Harrell-Davis estimator
# op=F, use single order statistic
#
#  LP=TRUE: The initial smooth is smoothed again using LOESS
#
plotit<-as.logical(plotit)
scat<-as.logical(scat)
m<-cbind(x,y)
if(ncol(m)!=2)stop("Must have exactly one predictor. For more than one, use qhdsm.")
m<-elimna(m)
x<-m[,1]
y<-m[,2]
if(eout && xout)stop("Not allowed to have eout=xout=T")
if(eout){
flag<-outfun(m,plotit=FALSE)$keep
m<-m[flag,]
}
if(xout){
flag<-outfun(x)$keep
m<-m[flag,]
}
x<-m[,1]
y<-m[,2]
rmd<-c(1:length(x))
if(pyhat)outval<-matrix(NA,ncol=length(qval),nrow=length(x))
if(scat)plot(x,y,xlab=xlab,ylab=ylab)
if(!scat)plot(x,y,type="n",xlab=xlab,ylab=ylab)
for(it in 1:length(qval)){
if(!op)for(i in 1:length(x))rmd[i]<-qest(y[near(x,x[i],fr)],q=qval[it])
if(op)for(i in 1:length(x))rmd[i]<-hd(y[near(x,x[i],fr)],q=qval[it])
if(pyhat)outval[,it]<-rmd
points(x,rmd,type="n")
sx<-sort(x)
xorder<-order(x)
sysm<-rmd[xorder]
if(LP)sysm=lplot(sx,sysm,pyhat=TRUE,plotit=FALSE)$yhat.values
lines(sx,sysm)
}
if(pyhat)output<-outval
if(!pyhat)output<-"Done"
list(output=output)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.