R/aov2depth.R

aov2depth <-
function(x1,x2,est=tmean,nboot=500,SEED=TRUE,nmin=12,CR=FALSE,
xlab=' DIF 1',ylab='DIF 2',zlab='DIF 3',alpha=.05,...){
#
# 2 by K ANOVA independent group (K levels not necessarily independent and 
#                                 not completely dependent
#
#   Main effect Factor A only
#
# Strategy: Use depth of zero based on estimated
# differences for each column  of the K levels of Factor B
# That is, testing no main effects for Factor A in 
# a manner that takes into account the pattern of the
# measures of location rather then simply averaging 
# across columns.
#
#  x1 can be a matrix with K columns corrspoding to groups, ditto for x2
#  Or x1 and x2 can have list mode.
#   Assuming x1 and x2 contain data for indepedendent groups.
#
if(is.matrix(x1)||is.data.frame(x1))x1=listm(x1)
if(is.matrix(x2)||is.data.frame(x2))x2=listm(x2)
J=length(x1)
if(J!=length(x2))stop('x1 and x2 should have same number of groups')
if(SEED)set.seed(2)
for(j in 1:J){
x1[[j]]=na.omit(x1[[j]])
x2[[j]]=na.omit(x2[[j]])
}
n1=mapply(x1,FUN=length)
n2=mapply(x2,FUN=length)
bplus=nboot+1
bvec1=matrix(NA,nrow=nboot,ncol=J)
bvec2=matrix(NA,nrow=nboot,ncol=J)
for(j in 1:J){
data1=matrix(sample(x1[[j]],size=n1[j]*nboot,replace=TRUE),nrow=nboot)
data2=matrix(sample(x2[[j]],size=n2[j]*nboot,replace=TRUE),nrow=nboot)
bvec1[,j]=apply(data1,1,est,...)
bvec2[,j]=apply(data2,1,est,...)
}
difb=bvec1-bvec2
est1=mapply(x1,FUN=est,...)
est2=mapply(x2,FUN=est,...)
dif=est1-est2
m1=var(difb)
nullvec=rep(0,J)
difz=rbind(difb,nullvec)
dis=mahalanobis(difz,dif,m1)
sig=sum(dis[bplus]<=dis)/bplus
if(CR){
dis2<-order(dis[1:nboot])
dis<-sort(dis)
critn<-floor((1-alpha)*nboot)
if(J==2){
plot(difb[,1],difb[,2],xlab=xlab,ylab=ylab)
points(0,0,pch=0)
xx<-difb[dis2[1:critn],]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,])
lines(xx[c(temp[1],temp[length(temp)]),])
}
if(J==3){
scatterplot3d(difb[dis2[1:critn],],xlab=xlab,ylab=ylab,zlab=zlab,tick.marks=TRUE)
}

}
list(p.value=sig,est1=est1,est2=est2,dif=dif,n1=n1,n2=n2)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.