R/lband.R

lband <-
function(x,y=NULL,alpha=.05,plotit=TRUE,sm=TRUE,op=1,ylab="delta",
xlab="x (first group)"){
#
#  Compute a confidence band for the shift function.
#  Assuming two dependent groups are being compared
#
#  See Lombard (2005, Technometrics, 47, 364-369)
#
#  if y=NA, assume x is a matrix with two columns or it has list mode
#
#  If plotit=TRUE, a plot of the shift function is created, assuming that
#  the graphics window has already been activated.
#
#  sm=T, plot of shift function is smoothed using:
#  expected frequency curve if op!=1
#  otherwise use S+ function lowess is used.
#
#  This function removes all missing observations.
#
#  When plotting, the median of x is marked with a + and the two
#  quartiles are marked with o.
#
if(!is.null(y[1]))x<-cbind(x,y)
if(is.list(x))x=matl(x)
if(ncol(x)!=2)stop("Should have two groups only")
m<-elimna(x)
y<-m[,2]
x<-m[,1]
n<-length(x)
crit<-nelderv2(m,1,lband.fun2,alpha=alpha)
plotit<-as.logical(plotit)
xsort<-sort(x)
ysort<-sort(y)
l<-0
u<-0
ysort[0]<-NA
ysort[n+1]<-NA
lsub<-c(1:n)-floor(sqrt(2*n)*crit)
usub<-c(1:n)+floor(sqrt(2*n)*crit)
for(ivec in 1:n){
isub<-max(0,lsub[ivec])
l[ivec]<-NA
if(isub>0)l[ivec]<-ysort[isub]-xsort[ivec]
isub<-min(n+1,usub[ivec])
u[ivec]<-NA
if(isub <= n)u[ivec]<-ysort[isub]-xsort[ivec]
}
num<-length(l[l>0 & !is.na(l)])+length(u[u<0 & !is.na(u)])
qhat<-c(1:n)/n
m<-cbind(qhat,l,u)
dimnames(m)<-list(NULL,c("qhat","lower","upper"))
if(plotit){
xsort<-sort(x)
ysort<-sort(y)
del<-0
for (i in 1:n)del[i]<-ysort[i]-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)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.