R/man2pb.R

man2pb <-
function(x,y,alpha=.05,nboot=NA,crit=NA){
#
#   Two-sample Behrens-Fisher problem.
#
#   For each of two independent groups,
#   have P measures for each subject. The goal is to compare the 20%
#   trimmed means of the first group to the trimmed means for the second;
#   this is done for each of the  P measures.
#
#   The percentile bootstrap method is used to
#   compute a .95, or .975, or .99 confidence interval.
#
#   Only 20% trimming is allowed.
#
#   x contains the data for the first group; it
#    can be an n by J matrix or it can have list mode.
#   y contains the data for the second group.
#
#   Vectors with missing values are eliminated from the analysis.
#
if(!is.list(x) && !is.matrix(x))stop("Data must be stored in a matrix or in list mode.")
if(!is.list(y) && !is.matrix(y))stop("Data must be stored in a matrix or in list mode.")
if(is.list(x)){
# put the data in an n by p matrix
matx<-matrix(0,length(x[[1]]),length(x))
for (j in 1:length(x))matx[,j]<-x[[j]]
}
if(is.list(y)){
# put the data in an n by p matrix
maty<-matrix(0,length(y[[1]]),length(y))
for (j in 1:length(y))maty[,j]<-y[[j]]
}
if(is.matrix(x)){
matx<-x
}
if(is.matrix(y)){
maty<-y
}
if(ncol(matx)!=ncol(maty))stop("The number of variables for group 1 is not equal to the number for group 2")
if(sum(is.na(matx)>=1))matx<-elimna(matx)
if(sum(is.na(maty)>=1))maty<-elimna(maty)
J<-ncol(matx)
connum<-ncol(matx)
if(is.na(nboot)){
if(ncol(matx)<=4)nboot<-2000
if(ncol(matx)>4)nboot<-5000
}
#
#  Determine critical value
#
if(ncol(matx)==2){
if(alpha==.05)crit<-.0125
if(alpha==.025)crit<-.0060
if(alpha==.01)crit<-.0015
}
if(ncol(matx)==3){
if(alpha==.05)crit<-.007
if(alpha==.025)crit<-.003
if(alpha==.01)crit<-.001
}
if(ncol(matx)==4){
if(alpha==.05)crit<-.0055
if(alpha==.025)crit<-.0020
if(alpha==.01)crit<-.0005
}
if(ncol(matx)==5){
if(alpha==.05)crit<-.0044
if(alpha==.025)crit<-.0016
if(alpha==.01)crit<-.0005
}
if(ncol(matx)==6){
if(alpha==.05)crit<-.0038
if(alpha==.025)crit<-.0018
if(alpha==.01)crit<-.0004
}
if(ncol(matx)==7){
if(alpha==.05)crit<-.0028
if(alpha==.025)crit<-.0010
if(alpha==.01)crit<-.0002
}
if(ncol(matx)==8){
if(alpha==.05)crit<-.0026
if(alpha==.025)crit<-.001
if(alpha==.01)crit<-.0002
}
if(ncol(matx)>8){
# Use an approximation of the critical value
if(alpha==.025)warning("Can't determine a critical value when alpha=.025 and the number of groups exceeds 8.")
nmin<-min(nrow(matx),nrow(maty))
if(alpha==.05){
if(nmin<100)wval<-smmcrit(60,ncol(matx))
if(nmin>=100)wval<-smmcrit(300,ncol(matx))
wval<-0-wval
crit<-pnorm(wval)
}
if(alpha==.01){
if(nmin<100)wval<-smmcrit01(60,ncol(matx))
if(nmin>=100)wval<-smmcrit01(300,ncol(matx))
wval<-0-wval
crit<-pnorm(wval)
}
}
if(is.na(crit))warning("Critical values can be determined for alpha=.05, .025 and .01 only")
icl<-ceiling(crit*nboot)
icu<-ceiling((1-crit)*nboot)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples. Please wait.")
bootx<-bootdep(matx,tr=.2,nboot)
booty<-bootdep(maty,tr=.2,nboot)
        #
        # Now have an nboot by J matrix of bootstrap values.
        #
test<-1
for (j in 1:connum){
test[j]<-sum(bootx[,j]<booty[,j])/nboot
if(test[j]>.5)test[j]<-1-test[j]
}
output <- matrix(0, connum, 5)
        dimnames(output) <- list(NULL, c("variable #", "psihat", "test",
                "ci.lower", "ci.upper"))
        tmeanx <- apply(matx, 2, mean, trim = .2)
        tmeany <- apply(maty, 2, mean, trim = .2)
        psi <- 1
        for(ic in 1:connum) {
                output[ic, 2] <- tmeanx[ic]-tmeany[ic]
                output[ic, 1] <- ic
                output[ic, 3] <- test[ic]
                temp <- sort(bootx[,ic]-booty[,ic])
print(length(temp))
                output[ic, 4] <- temp[icl]
                output[ic, 5] <- temp[icu]
        }
        list(output = output, crit.value = crit)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.