R/smgvcr.R

smgvcr <-
function(m,nullv=rep(0,ncol(m)),SEED=TRUE,op=0,
nboot=500,plotit=TRUE){
#
# m is an n by p matrix
#
# Test hypothesis that estimand of the MGV estimator
# is equal to the null value, which defaults to zero vector.
# The level of the test is .05.
#
# Argument op: See function outmgv
#
if(SEED)set.seed(2)
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<-outmgv(mm,plotit=FALSE,op=op)$keep
val[j,]<-apply(mm[temp,],2,mean)
}
temp<-mgvar(rbind(val,nullv),op=op)
flag2<-is.na(temp)
if(sum(flag2)>0)temp[flag2]<-0
sig.level<-sum(temp[nboot+1]<temp[1:nboot])/nboot
if(ncol(m)==2 && plotit){
plot(val[,1],val[,2],xlab="VAR 1",ylab="VAR 2")
temp3<-mgvmean(m,op=op)
points(temp3[1],temp3[2],pch="+")
ic<-round((1-crit.level)*nboot)
temp<-mgvar(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.