R/smean2v2.R

smean2v2 <-
function(m1,m2,nullv=rep(0,ncol(m1)),cop=3,MM=FALSE,SEED=NA,
nboot=500,plotit=TRUE,MC=FALSE,STAND=TRUE){
#
# m is an n by p matrix
#
# For two independent groups,
# test hypothesis that multivariate skipped estimators
# are all equal.
#
# The level of the test is .05.
#
# Skipped estimator is used, i.e.,
# eliminate outliers using a projection method
# That is, determine center of data using:
#
# cop=1 Donoho-Gasko median,
# cop=2 MCD,
# cop=3 marginal medians.
# cop=4 MVE
#
# For each point
# consider the line between it and the center
# project all points onto this line, and
# check for outliers using
#
# MM=F, a boxplot rule.
# MM=T, rule based on MAD and median
#
# Repeat this for all points. A point is declared
# an outlier if for any projection it is an outlier
# using a modification of the usual boxplot rule.
#
# Eliminate any outliers and compute means
#  using remaining data.
#
if(ncol(m1) != ncol(m2)){
stop("Number of variables in group 1 does not equal the number in group 2.")
}
if(is.na(SEED))set.seed(2)
if(!is.na(SEED))set.seed(SEED)
m1<-elimna(m1)
m2<-elimna(m2)
n1<-nrow(m1)
n2<-nrow(m2)
n<-min(c(n1,n2))
crit.level<-.05
if(n<=120)crit.level<-.045
if(n<=80)crit.level<-.04
if(n<=60)crit.level<-.035
if(n<=40)crit.level<-.03
if(n<=30)crit.level<-.025
if(n<=20)crit.level<-.02
val<-matrix(NA,ncol=ncol(m1),nrow=nboot)
est=smean(m1)-smean(m2)
for(j in 1: nboot){
data1<-sample(n1,size=n1,replace=TRUE)
data2<-sample(n2,size=n2,replace=TRUE)
mm1<-m1[data1,]
temp<-outpro(mm1,plotit=FALSE,cop=cop,STAND=STAND)$keep
v1<-apply(mm1[temp,],2,mean)
mm2<-m2[data2,]
temp<-outpro(mm2,plotit=FALSE,cop=cop,STAND=STAND)$keep
v2<-apply(mm2[temp,],2,mean)
val[j,]<-v1-v2
}
if(!MC)temp<-pdis(rbind(val,nullv))
if(MC)temp<-pdisMC(rbind(val,nullv))
sig.level<-sum(temp[nboot+1]<temp[1:nboot])/nboot
if(ncol(m1)==2 && plotit){
plot(val[,1],val[,2],xlab="VAR 1",ylab="VAR 2")
if(!MC)temp3<-smean(m1,cop=cop)-smean(m2,cop=cop)
if(MC)temp3<-smeanMC(m1,cop=cop)-smeanMC(m2,cop=cop)
points(temp3[1],temp3[2],pch="+")
ic<-round((1-crit.level)*nboot)
if(!MC)temp<-pdis(val)
if(MC)temp<-pdisMC(val)
temp.dis<-order(temp)
xx<-val[temp.dis[1:ic],]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,])
lines(xx[c(temp[1],temp[length(temp)]),])
}
list(p.value=sig.level,crit.level=crit.level)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.