R/mee.R

mee <-
function(x,y,alpha=.05){
#
#  For two independent groups, compute a 1-\alpha confidence interval
#  for p=P(X<Y) using Mee's method, which assumes there are no ties.
#  If ties are detected among the pooled observations, a warning message is
#  printed. The type I error probability might exceed the nominal level by
#  an unacceptable amount. Also, mee(x,y) might give  different  results than
#  mee(y,x). (One might reject and the other might not.)
#
x<-x[!is.na(x)]  # Remove missing values from x
y<-y[!is.na(y)]  # Remove missing values from y
xy<-c(x,y)
tiexy<-sum(abs(c(1:length(xy))-sort(rank(xy))))
if(tiexy > 0){print("Warning: Tied values detected")
print("so even if distributions are identical,")
print("P(X<Y) is not necessarily equal to .5")
}
u<-outer(x,y,FUN="<")
p1<-0
p2<-0
for (j in 1:length(y)){
temp<-outer(u[,j],u[,j])
p1<-p1+sum(temp)-sum(u[,j]*u[,j])
}
for (i in 1: length(x)){
temp<-outer(u[i,],u[i,])
p2<-p2+sum(temp)-sum(u[i,]*u[i,])
}
p<-sum(u)/(length(x)*length(y))
p1<-p1/(length(x)*length(y)*(length(x)-1))
p2<-p2/(length(x)*length(y)*(length(y)-1))
b1<-(p1-p^2)/(p-p^2)
b2<-(p2-p^2)/(p-p^2)
A<-((length(x)-1)*b1+1)/(1-1/length(y))+((length(y)-1)*b2+1)/(1-1/length(x))
nhat<-length(x)*length(y)/A
crit<-(qnorm(1-alpha/2))^2/nhat
D<-sqrt(crit*(p*(1-p)+.25*crit))
low<-(p+.5*crit-D)/(1+crit)
hi<-(p+.5*crit+D)/(1+crit)
list(phat=p,ci=c(low,hi))
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.