R/ancdet2C.R

ancdet2C <-
function(x1,y1,x2,y2,fr1=1,fr2=1,tr=.2,
alpha=.05,plotit=TRUE,op=FALSE,pts=NA,sm=FALSE,FRAC=.5,
pr=TRUE,xout=FALSE,outfun=outpro,MC=FALSE,
p.crit=NULL,nreps=2000,SEED=TRUE,FAST=TRUE,ticktype='detail',
xlab='X1',ylab='X2',zlab='Y',pch1='*',pch2='+',...){
#
#  Like ancdet, only two covariate values can be used in conjunction
#  with method MC2 for pick the covariate values. That is,
#  using the deepest half of the covariate values. 
#
#  politit=TRUE. Plot covariate points. Significant points are indicated by
#  pch='+'
#
if(ncol(as.matrix(x1))!=2)stop('Two covariates only can be used')
if(is.null(p.crit))set.seed(2)
xy=elimna(cbind(x1,y1)) 
x1=xy[,1:2]
y1=xy[,3]
xy=elimna(cbind(x2,y2))
x2=xy[,1:2]
y2=xy[,3]
if(min(length(y1),length(y2))<50)stop('The minimum sample size must be greater than or equal to  50')
if(xout){
flag<-outfun(x1,plotit=FALSE,...)$keep
x1<-x1[flag,]
y1<-y1[flag]
flag<-outfun(x2,plotit=FALSE,...)$keep
x2<-x2[flag,]
y2<-y2[flag]
}
if(FAST){
if(FRAC==.5){
if(is.null(p.crit)){
if(alpha==.05){
nv=c(50,  55,  60,  70,  80, 100, 200, 300, 400, 500, 600,800)
pv=c(0.004585405, 0.003199894, 0.002820089, 0.002594342, 0.002481210, 0.001861313,
  0.001419821, 0.001423000, 0.001313700, 0.001351900, 0.001075, 0.00095859)
 n1= length(y1)
 n2=length(y2)
p.crit=(lplot.pred(1/nv,pv,1/n1)$yhat+lplot.pred(1/nv,pv,1/n2)$yhat)/2
if(max(n1,n2)>max(nv)){
p.crit=min(pv)
print('Warning: p.crit has not been computed exactly for sample sizes greater than 800')
if(n1>800)p.crit1=regYhat(1/pv[8:12,1],pv[8:12,2],1/n1)
if(n1<=800)p.crit1=lplot.pred(1/nv,pv,1/n1)$yhat
if(n2>800)p.crit1=regYhat(1/pv[8:12,1],pv[8:12,2],1/n2)
if(n2<=800)p.crit1=lplot.pred(1/nv,pv,1/n2)$yhat
p.crit=(p.crit1+p.crit2)/2
}
}}}}
res1=ancov2COV(x1,y1,x2,y2,DETAILS=TRUE,pr=FALSE,FRAC=FRAC)
if(is.null(p.crit))p.crit=ancdet2C.pv(length(y1),length(y2),MC=MC,nreps=nreps,
SEED=SEED)
LL=length(ncol(res1$all.results))
if(LL==1)num.sig=sum(res1$all.results[,3]<=p.crit)
if(LL==0)num.sig=NA
sig.points=NA
if(LL==1){
flag=res1$all.results[,3]<=p.crit
sig.points=res1$all.points.used[flag,1:2]
}
if(plotit){
if(!op){
if(pr)print('To plot the estimated differences for the covariate points used, set op=TRUE')
if(LL==0)plot(res1$all.points.used[,1],res1$all.points.used[,2],xlab=xlab,ylab=ylab)
if(LL==1){
plot(res1$all.points.used[,1],res1$all.points.used[,2],type='n',xlab=xlab,ylab=ylab)
points(res1$all.points.used[!flag,1],res1$all.points.used[!flag,2],pch=pch1)
points(res1$all.points.used[flag,1],res1$all.points.used[flag,2],pch=pch2)
}}
if(op)
lplot(res1$all.points.used[,1:2],res1$all.points.used[,3],xlab=xlab,ylab=ylab,zlab=zlab,ticktype=ticktype)
}
list(num.sig=num.sig,p.crit=p.crit,points.used=cbind(res1$all.points.used[,1:3],res1$all.results),sig.points=sig.points)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.