R/bmp.R

In musto101/wilcox_R: What the package does (short line)

bmp <-
function(x,y,alpha=.05,crit=NA,plotit=FALSE,pop=0,fr=.8,xlab="",ylab=""){
#
# Brunner and Munzel (2000) heteroscedastic analog of WMW test.
#
# plotit=TRUE causes a plot of the difference scores to be created
x<-x[!is.na(x)]  # Remove any missing values
y<-y[!is.na(y)]
n1<-length(x)
n2<-length(y)
N<-n1+n2
n1p1<-n1+1
flag1<-c(1:n1)
flag2<-c(n1p1:N)
R<-rank(c(x,y))
R1<-mean(R[flag1])
R2<-mean(R[flag2])
Rg1<-rank(x)
Rg2<-rank(y)
S1sq<-sum((R[flag1]-Rg1-R1+(n1+1)/2)^2)/(n1-1)
S2sq<-sum((R[flag2]-Rg2-R2+(n2+1)/2)^2)/(n2-1)
sig1<-S1sq/n2^2
sig2<-S2sq/n1^2
se<-sqrt(N)*sqrt(N*(sig1/n1+sig2/n2))
bmtest<-(R2-R1)/se
phat<-(R2-(n2+1)/2)/n1
dhat<-1-2*phat
df<-(S1sq/n2 + S2sq/n1)^2/((S1sq/n2)^2/(n1-1)+(S2sq/n1)^2/(n2-1))
sig<-2 * (1 - pt(abs(bmtest),df))
if(is.na(crit))vv<-qt(alpha/2,df)
if(!is.na(crit))vv<-crit
ci.p<-c(phat+vv*se/N,phat-vv*se/N)
if(plotit){
msave<-outer(x,y,FUN="-")
if(pop==0){
if(length(x)*length(y)>2500){
print("Product of sample sizes exceeds 2500.")
print("Execution time might be high when plotting and when using pop=1")
print("If this is case, might consider changing the argument pop or using plotit=F")
}
akerd(as.vector(msave),fr=fr)
}
if(pop==1)rdplot(as.vector(msave),fr=fr,xlab=xlab,ylab=ylab)
if(pop==2)kdplot(as.vector(msave),rval=rval,xlab=xlab,ylab=ylab)
if(pop==3)boxplot(as.vector(msave))
if(pop==4)stem(as.vector(msave))
if(pop==5)hist(as.vector(msave))
if(pop==6)skerd(as.vector(msave),xlab=xlab,ylab=ylab)
}
list(test.stat=bmtest,phat=phat,dhat=dhat,p.value=sig,ci.p=ci.p,df=df)
}