R/mdiflcr.R

mdiflcr <-
function(m1,m2,tr=.5,nullv=rep(0,ncol(m1)),plotit=TRUE,
SEED=TRUE,pop=1,fr=.8,nboot=600){
#
# For two independent groups, let D=X-Y.
# Let theta_D be median of marginal distributions
# Goal: Test theta_D=0
#
# This is a multivariate analog of Wilcoxon-Mann-Whitney method
# Only alpha=.05 can be used.
#
# When plotting:
# pop=1 Use scatterplot
# pop=2 Use expected frequency curve.
# pop=3 Use adaptive kernel density
#
if(!is.matrix(m1))stop("m1 is not a matrix")
if(!is.matrix(m2))stop("m2 is not a matrix")
if(ncol(m1)!=ncol(m2))stop("number of columns for m1 and m2 are not equal")
n1<-nrow(m1)
n2<-nrow(m2)
if(SEED)set.seed(2)
data1 <- matrix(sample(n1, size = n1 * nboot, replace = T), nrow = nboot)
data2 <- matrix(sample(n2, size = n2 * nboot, replace = T), nrow = nboot)
bcon <- matrix(NA, ncol = ncol(m1), nrow = nboot)
for(j in 1:nboot)bcon[j,]<-mdifloc(m1[data1[j,],],m2[data2[j,],],est=lloc,tr=tr)
tvec<-mdifloc(m1,m2,est=lloc,tr=tr)
tempcen <- apply(bcon, 1, mean)
smat <- var(bcon - tempcen + tvec)
temp <- bcon - tempcen + tvec
bcon <- rbind(bcon, nullv)
dv <- mahalanobis(bcon, tvec, smat)
bplus <- nboot + 1
sig.level <- 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
if(plotit && ncol(m1)==2){
if(pop==2)rdplot(mdif,fr=fr)
if(pop==1){
plot(mdif[,1],mdif[,2],xlab="VAR 1",ylab="VAR 2",type="n")
points(mdif[,1],mdif[,2],pch=".")
points(center[1],center[2],pch="o")
points(0,0,pch="+")
}
if(pop==3)akerdmul(mdif,fr=fr)
}
list(p.value=sig.level,center=tvec)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.