R/smeancr.R

smeancr <-
function(m,nullv=rep(0,ncol(m)),cop=3,MM=FALSE,SEED=NA,
nboot=500,plotit=TRUE,xlab="VAR 1",ylab="VAR 2",STAND=TRUE){
#
# m is an n by p matrix
#
# Test hypothesis that multivariate skipped estimators
# are all equal to the null value, which defaults to zero.
# The level of the test is .05.
#
# 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(is.na(SEED))set.seed(2)
if(!is.na(SEED))set.seed(SEED)
m<-elimna(m)
n<-nrow(m)
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
data<-matrix(sample(n,size=n*nboot,replace=TRUE),nrow=nboot)
val<-matrix(NA,ncol=ncol(m),nrow=nboot)
for(j in 1: nboot){
mm<-m[data[j,],]
temp<-outpro(mm,plotit=FALSE,cop=cop,STAND=STAND)$keep
val[j,]<-apply(mm[temp,],2,mean)
}
temp<-pdis(rbind(val,nullv))
sig.level<-sum(temp[nboot+1]<temp[1:nboot])/nboot
if(ncol(m)==2 && plotit){
plot(val[,1],val[,2],xlab=xlab,ylab=ylab)
temp3<-smean(m,cop=cop,STAND=STAND)
points(temp3[1],temp3[2],pch="+")
ic<-round((1-crit.level)*nboot)
temp<-pdis(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.