R/lsqs3.R

lsqs3 <-
function(x,y,plotit=TRUE,cop=2,ap.dep=FALSE,v2=FALSE,pv=FALSE,SEED=TRUE,nboot=1000,ypch="o",xpch="+"){
#
#  Compute the typical depth of x in y,
#  Compute the typical depth of y in x,
#  use the maximum of the two typical depths
#  as a test statistic.
#  This method is designed to be sensitive to
#  shifts in location.
#
# Use Tukey's depth; bivariate case only.
#
# cop=2 use MCD location estimator when
# computing depth with function fdepth
# cop=3 uses medians
# cop=3 uses MVE
#
#  xpch="+" means when plotting the data, data from the first
#  group are indicated by a +
#  ypch="o" are data from the second group
#
if(is.list(x))x<-matl(x)
if(is.list(y))y<-matl(y)
x<-elimna(x)
y<-elimna(y)
x<-as.matrix(x)
y<-as.matrix(y)
nx=nrow(x)
ny=nrow(y)
if(ncol(x) != ncol(y))stop("Number of variables not equal")
disyx<-NA # depth of y in x
disxy<-NA # depth of x in y
#
if(ncol(x)==2){
if(plotit){
plot(rbind(x,y),type="n",xlab="VAR 1",ylab="VAR 2")
points(x,pch=xpch)
points(y,pch=ypch)
if(nrow(x)>50){
if(!ap.dep){
print("If execution time is high, might use ap.dep=F")
}
if(!ap.dep)temp<-depth2(x,plotit=FALSE)
if(ap.dep)temp<-fdepth(x,plotit=FALSE,cop=cop)
}
if(!ap.dep)temp<-depth2(x,plotit=FALSE)
if(ap.dep)temp<-fdepth(x,plotit=FALSE,cop=cop)
flag<-(temp>=median(temp))
xx<-x[flag,]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,])
lines(xx[c(temp[1],temp[length(temp)]),])
if(ap.dep)temp<-fdepth(y,plotit=FALSE,cop=cop)
if(!ap.dep)temp<-depth2(y,plotit=FALSE)
if(!ap.dep)temp<-depth2(y,plotit=FALSE)
if(!ap.dep)temp<-fdepth(y,plotit=FALSE)
flag<-(temp>=median(temp))
xx<-y[flag,]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
flag<-(temp>=median(temp))
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,],lty=2)
lines(xx[c(temp[1],temp[length(temp)]),],lty=2)
}
tempyx<-NA
tempxy<-NA
if(ap.dep)tempyx<-fdepth(x,y,plotit=FALSE,cop=cop)
if(!ap.dep)tempyx<-depth2(x,y,plotit=FALSE)
if(ap.dep)tempxy<-fdepth(y,x,plotit=FALSE,cop=cop)
tempxy<-depth2(y,x,plotit=FALSE)
}
if(ncol(x)==1){
tempyx<-unidepth(as.vector(x),as.vector(y))
tempxy<-unidepth(as.vector(y),as.vector(x))
}
if(ncol(x)>2){
if(!v2){
tempxy<-fdepth(y,x,plotit=FALSE,cop=cop)
tempyx<-fdepth(x,y,plotit=FALSE,cop=cop)
}
if(v2){
tempxy<-fdepthv2(y,x,plotit=FALSE,cop=cop)
tempyx<-fdepthv2(x,y,plotit=FALSE,cop=cop)
}}
qhatxy<-mean(tempxy)
qhatyx<-mean(tempyx)
qhat<-max(c(qhatxy,qhatyx))
n1<-nrow(x)
n2<-nrow(y)
nv<-(3*min(c(n1,n2))+max(c(n1,n2)))/4
if(ncol(x)==1)crit<-.2536-.4578/sqrt(nv)
if(ncol(x)==2)crit<-.1569-.3/sqrt(nv)
if(ncol(x)==3)crit<-.0861-.269/sqrt(nv)
if(ncol(x)==4)crit<-.054-.1568/sqrt(nv)
if(ncol(x)==5)crit<-.0367-.0968/sqrt(nv)
if(ncol(x)==6)crit<-.0262-.0565/sqrt(nv)
if(ncol(x)==7)crit<-.0174-.0916/sqrt(nv)
if(ncol(x)>7)crit<-.013
rej<-"Fail to reject"
if(qhat<=crit)rej<-"Reject"
testv=NULL
pval=NULL
if(pv){
if(SEED)set.seed(2)
rej="NULL"
for(i in 1:nboot)testv[i]=lsqs3.sub(rmul(n1,ncol(x)),rmul(n2,ncol(x)),cop=cop,ap.dep=ap.dep,v2=v2,)$test
pval=mean(qhat>=testv)
}
list(n1=nx,n2=ny,avg.depth.of.x.in.y=qhatxy,avg.depth.of.y.in.x=qhatyx,test=qhat,crit=crit,Decision=rej,p.value=pval)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.