R/ancpb.R

ancpb <-
function(x1,y1,x2,y2,est=hd,pts=NA,fr1=1,fr2=1,nboot=NA,alpha=.05,xout=FALSE,outfun=outpro,plotit=TRUE,LP=TRUE,xlab='X',ylab='Y',pch1='*',pch2='+',...){
#
# Compare two independent  groups using an ancova method
# with a percentile bootstrap combined with a running interval
# smooth.
#
#  Assume data are in x1 y1 x2 and y2
#  Comparisons are made at the design points contained in the vector
#  pts
#
flag.est=FALSE
if(identical(est,onestep))flag.est=TRUE
if(flag.est)LP=FALSE   # Get an error when using onestep in conjunction with LP=T
if(identical(est,mom))flag.est=TRUE
xy1=elimna(cbind(x1,y1))
x1=xy1[,1]
y1=xy1[,2]
xy2=elimna(cbind(x2,y2))
x2=xy2[,1]
y2=xy2[,2]
if(xout){
flag<-outfun(x1,...)$keep
x1<-x1[flag]
y1<-y1[flag]
flag<-outfun(x2,...)$keep
x2<-x2[flag]
y2<-y2[flag]
}
npt<-5
gv1<-vector("list")
if(is.na(pts[1])){
isub<-c(1:5)  # Initialize isub
test<-c(1:5)
xorder<-order(x1)
y1<-y1[xorder]
x1<-x1[xorder]
xorder<-order(x2)
y2<-y2[xorder]
x2<-x2[xorder]
n1<-1
n2<-1
vecn<-1
for(i in 1:length(x1))n1[i]<-length(y1[near(x1,x1[i],fr1)])
for(i in 1:length(x1))n2[i]<-length(y2[near(x2,x1[i],fr2)])
for(i in 1:length(x1))vecn[i]<-min(n1[i],n2[i])
sub<-c(1:length(x1))
isub[1]<-min(sub[vecn>=12])
isub[5]<-max(sub[vecn>=12])
isub[3]<-floor((isub[1]+isub[5])/2)
isub[2]<-floor((isub[1]+isub[3])/2)
isub[4]<-floor((isub[3]+isub[5])/2)
mat<-matrix(NA,5,3)
dimnames(mat)<-list(NULL,c("X","n1","n2"))
for (i in 1:5){
j<-i+5
temp1<-y1[near(x1,x1[isub[i]],fr1)]
temp2<-y2[near(x2,x1[isub[i]],fr2)]
temp1<-temp1[!is.na(temp1)]
temp2<-temp2[!is.na(temp2)]
mat[i,1]<-x1[isub[i]]
mat[i,2]<-length(temp1)
mat[i,3]<-length(temp2)
gv1[[i]]<-temp1
gv1[[j]]<-temp2
}
I1<-diag(npt)
I2<-0-I1
con<-rbind(I1,I2)
if(flag.est)test<-pbmcp(gv1,alpha=alpha,nboot=nboot,est=est,con=con,...)
if(!flag.est)test<-linconpb(gv1,alpha=alpha,nboot=nboot,est=est,con=con,...)
}
#
if(!is.na(pts[1])){
npt<-length(pts)
n1<-1
n2<-1
vecn<-1
for(i in 1:length(pts)){
n1[i]<-length(y1[near(x1,pts[i],fr1)])
n2[i]<-length(y2[near(x2,pts[i],fr2)])
}
mat<-matrix(NA,length(pts),3)
dimnames(mat)<-list(NULL,c("X","n1","n2"))
gv<-vector("list",2*length(pts))
for (i in 1:length(pts)){
j<-i+npt
temp1<-y1[near(x1,pts[i],fr1)]
temp2<-y2[near(x2,pts[i],fr2)]
temp1<-temp1[!is.na(temp1)]
temp2<-temp2[!is.na(temp2)]
mat[i,1]<-pts[i]
if(length(temp1)<=5)paste("Warning, there are",length(temp1)," points corresponding to the design point X=",pts[i])
if(length(temp2)<=5)paste("Warning, there are",length(temp2)," points corresponding to the design point X=",pts[i])
mat[i,2]<-length(temp1)
mat[i,3]<-length(temp2)
gv1[[i]]<-temp1
gv1[[j]]<-temp2
}
I1<-diag(npt)
I2<-0-I1
con<-rbind(I1,I2)
if(flag.est)test<-pbmcp(gv1,alpha=alpha,nboot=nboot,est=est,con=con,...)
if(!flag.est)test<-linconpb(gv1,alpha=alpha,nboot=nboot,est=est,con=con,...)
}
if(plotit){
runmean2g(x1,y1,x2,y2,fr=fr1,est=est,LP=LP,xlab=xlab,ylab=ylab,pch1=pch1,pch2=pch2,...)
}
list(mat=mat,output=test$output,con=test$con,num.sig=test$num.sig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.