R/oancpb.R

oancpb <-
function(x1,y1,x2,y2,est=tmean,tr=.2,pts=NA,fr1=1,fr2=1,nboot=600,
alpha=.05,plotit=TRUE,SEED=TRUE,PRO=FALSE,...){
#
# Compare two independent  groups using an ancova method
# with a percentile bootstrap combined with a running interval
# smooth.
#
#  CURRENTLY SEEMS THAT THE R FUNCTION ancGLOB is better.
#
#  This function performs an omnibus test using data corresponding
#  to K design points  specified by the argument pts. If
#  pts=NA, K=5 points are chosen for you (see Introduction to Robust
#  Estimation and Hypothesis Testing.)
#  Null hypothesis is that conditional distribution of Y, given X for first
#  group,  minus the conditional distribution of Y, given X for second
#  group is equal to zero.
#  The strategy is to choose K specific X values
#  and then test the hypothesis that all K differences are zero.
#
#  If you want to choose specific X values, Use the argument
#  pts
#  Example: pts=c(1,3,5) will use X=1, 3 and 5.
#
#  For multiple comparisons using these J points, use ancpb
#
#  Assume data are in x1 y1 x2 and y2
#
# PRO=F, means Mahalanobis distance is used.
# PRO=T, projection distance is used.
#
#  fr1 and fr2 are the spans used to fit a smooth to the data.
#
stop('USE ancGLOB')
#
#
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)
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)]
gv1[[i]]<-temp1
gv1[[j]]<-temp2
}
#
loc<-NA
if(SEED)set.seed(2)
bvec<-matrix(NA,nrow=nboot,ncol=5)
for(j in 1:5){
k<-j+5
loc[j]<-est(gv1[[j]])-est(gv1[[k]])
xx<-matrix(sample(gv1[[j]],size=length(gv1[[j]])*nboot,replace=TRUE),
nrow=nboot)
yy<-matrix(sample(gv1[[k]],size=length(gv1[[k]])*nboot,replace=TRUE),
nrow=nboot)
bvec[,j]<-apply(xx,1,FUN=est,...)-apply(yy,1,FUN=est,...)
}
nullv<-rep(0,5)
if(!PRO){
mvec<-apply(bvec,2,FUN=mean)
m1<-var(t(t(bvec)-mvec+loc))
temp<-mahalanobis(rbind(bvec,nullv),loc,m1)
}
if(PRO){
temp<-pdis(rbind(bvec,nullv))
}
sig.level<-sum(temp[nboot+1]<temp[1:nboot])/nboot
}
if(!is.na(pts[1])){
npt<-length(pts)
n1<-1
n2<-1
vecn<-1
mat<-matrix(NA,nrow=2*length(pts),ncol=3)
for(i in 1:length(pts)){
n1[i]<-length(y1[near(x1,pts[i],fr1)])
n2[i]<-length(y2[near(x2,pts[i],fr2)])
}
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)<=10)print(paste("Warning, there are",length(temp1)," points corresponding to the design point X=",pts[i]))
if(length(temp2)<=10)print(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
}
loc<-NA
if(SEED)set.seed(2)
bvec<-matrix(NA,nrow=nboot,ncol=npt)
for(j in 1:npt){
k<-j+npt
loc[j]<-est(gv1[[j]])-est(gv1[[k]])
xx<-matrix(sample(gv1[[j]],size=length(gv1[[j]])*nboot,replace=TRUE),
nrow=nboot)
yy<-matrix(sample(gv1[[k]],size=length(gv1[[k]])*nboot,replace=TRUE),
nrow=nboot)
bvec[,j]<-apply(xx,1,FUN=est,...)-apply(yy,1,FUN=est,...)
}
nullv<-rep(0,npt)
if(!PRO){
mvec<-apply(bvec,2,FUN=mean)
m1<-var(t(t(bvec)-mvec+loc))
temp<-mahalanobis(rbind(bvec,nullv),loc,m1)
}
if(PRO)temp<-pdis(rbind(bvec,nullv))
sig.level<-sum(temp[nboot+1]<temp[1:nboot])/nboot
}
if(plotit)runmean2g(x1,y1,x2,y2,fr=fr1,est=est,...)
list(p.value=sig.level)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.