R/twomanbt.R

twomanbt <-
function(x,y,tr=.2,alpha=.05,nboot=599){
#
#   Two-sample Behrens-Fisher problem.
#
#   For each of two independent groups,
#   have p measures for each subject. The goal is to compare the
#   trimmed means of the first measure, the trimmed means for the second
#   and so on.   So there are a total of p comparisons between the two
#   groups, one for each measure.
#
#   The percentile t bootstrap method is used to
#   compute a .95 confidence interval.
#
#   By default, 20% trimming is used with B=599 bootstrap samples.
#
#   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.
#
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 one is not equal to the number for group 2")
if(sum(is.na(mat)>=1))stop("Missing values are not allowed.")
J<-ncol(mat)
connum<-ncol(matx)
bvec<-matrix(0,connum,nboot)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
xcen<-matrix(0,nrow(matx),ncol(matx))
ycen<-matrix(0,nrow(maty),ncol(maty))
for (j in 1:connum)xcen[,j]<-matx[,j]-mean(matx[,j],tr) #Center data
for (j in 1:connum)ycen[,j]<-maty[,j]-mean(maty[,j],tr) #Center data
print("Taking bootstrap samples. Please wait.")
bootx<-sample(nrow(matx),size=nrow(matx)*nboot,replace=TRUE)
booty<-sample(nrow(maty),size=nrow(maty)*nboot,replace=TRUE)
matval<-matrix(0,nrow=nboot,ncol=connum)
for (j in 1:connum){
datax<-matrix(xcen[bootx,j],ncol=nrow(matx))
datay<-matrix(ycen[booty,j],ncol=nrow(maty))
paste("Working on variable", j)
top<- apply(datax, 1., mean, tr) - apply(datay, 1., mean, tr)
botx <- apply(datax, 1., trimse, tr)
boty <- apply(datay, 1., trimse, tr)
matval[,j]<-abs(top)/sqrt(botx^2. + boty^2.)
}
bvec<-apply(matval,1,max)
icrit<-round((1-alpha)*nboot)
bvec<-sort(bvec)
crit<-bvec[icrit]
psihat<-matrix(0,ncol=4,nrow=connum)
dimnames(psihat)<-list(NULL,c("con.num","psihat","ci.lower","ci.upper"))
test<-matrix(0,ncol=3,nrow=connum)
dimnames(test)<-list(NULL,c("con.num","test","se"))
for(j in 1:ncol(matx)){
temp<-yuen(matx[,j],maty[,j],tr=tr)
test[j,1]<-j
test[j,2]<-abs(temp$test)
test[j,3]<-temp$se
psihat[j,1]<-j
psihat[j,2]<-mean(matx[,j],tr)-mean(maty[,j])
psihat[j,3]<-mean(matx[,j],tr)-mean(maty[,j])-crit*temp$se
psihat[j,4]<-mean(matx[,j],tr)-mean(maty[,j])+crit*temp$se
}
list(psihat=psihat,teststat=test,critical.value=crit)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.