R/sband.R

sband <-
function(x,y,
crit=1.36*sqrt((length(x)+length(y))/(length(x)*length(y))),flag=FALSE,plotit=TRUE,
sm=TRUE,op=1,xlab="First Group",ylab="Delta",pr=TRUE){
#
#  Compute a confidence band for the shift function.
#  Assuming two independent groups are being compared
#
#  The default critical value is the approximate .05 critical value.
#
#  If flag is F (false), the exact significance level is not computed.
#
#  If plotit=TRUE, a plot of the shift function is created, assuming that
#  the graphics window has already been activated.
#
#  This function removes all missing observations.
#
#  When plotting, the median of x is marked with a + and the two
#  quaratiles are marked with o.
#
#  sm=T, shift function is smoothed using:
#  op!=1, running interval smoother,
#  otherwise use lowess.
#
x<-x[!is.na(x)]  # Remove missing values from x.
y<-y[!is.na(y)]  # Remove missing values from y.
plotit<-as.logical(plotit)
flag<-as.logical(flag)
pc<-NA
if(flag){
if(pr)print("Computing the exact value of the probability coverage")
pc<-1-kssig(length(x),length(y),crit)
chk=sum(duplicated(x,y))
if(chk>0)pc=1-kstiesig(x,y,crit)
}
xsort<-sort(x)
ysort<-c(NA,sort(y))
l<-0
u<-0
ysort[length(y)+1+1]<-NA
for(ivec in 1:length(x))
{
isub<-max(0,ceiling(length(y)*(ivec/length(x)-crit)))
l[ivec]<-ysort[isub+1]-xsort[ivec]
isub<-min(length(y)+1,floor(length(y)*(ivec/length(x)+crit))+1)
u[ivec]<-ysort[isub+1]-xsort[ivec]
}
num<-length(l[l>0 & !is.na(l)])+length(u[u<0 & !is.na(u)])
qhat<-c(1:length(x))/length(x)
m<-matrix(c(qhat,l,u),length(x),3)
dimnames(m)<-list(NULL,c("qhat","lower","upper"))
if(plotit){
xsort<-sort(x)
ysort<-sort(y)
del<-0
for (i in 1:length(x)){
ival<-round(length(y)*i/length(x))
if(ival<=0)ival<-1
if(ival>length(y))ival<-length(y)
del[i]<-ysort[ival]-xsort[i]
}
xaxis<-c(xsort,xsort)
yaxis<-c(m[,1],m[,2])
allx<-c(xsort,xsort,xsort)
ally<-c(del,m[,2],m[,3])
temp2<-m[,2]
temp2<-temp2[!is.na(temp2)]
plot(allx,ally,type="n",ylab=ylab,xlab=xlab)
ik<-rep(F,length(xsort))
if(sm){
if(op==1){
ik<-duplicated(xsort)
del<-lowess(xsort,del)$y
}
if(op!=1)del<-runmean(xsort,del,pyhat=TRUE)
}
lines(xsort[!ik],del[!ik])
lines(xsort,m[,2],lty=2)
lines(xsort,m[,3],lty=2)
temp<-summary(x)
text(temp[3],min(temp2),"+")
text(temp[2],min(temp2),"o")
text(temp[5],min(temp2),"o")
}
list(m=m,crit=crit,numsig=num,pc=pc)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.