R/pbcan.R

pbcan <-
function(x,nboot=1000,grp=NA,est=onestep,...){
#
#   Test the hypothesis that J independent groups have
#   equal measures of location using the percentile bootstrap method.
#   in conjunction with a partially centering technique.
#
#   The data are assumed to be stored in x
#   which either has list mode or is a matrix.  In the first case
#   x[[1]] contains the data for the first group, x[[2]] the data
#   for the second group, etc. Length(x)=the number of groups = J.
#   If stored in a matrix, the columns of the matrix correspond
#   to groups.
#
#   est is the measure of location and defaults to an M-estimator
#   ... can be used to set optional arguments associated with est
#
#   The argument grp can be used to analyze a subset of the groups
#   Example: grp=c(1,3,5) would compare groups 1, 3 and 5.
#
#   Missing values are allowed.
#
if(is.matrix(x))x<-listm(x)
if(!is.list(x))stop("Data must be stored in list mode or in matrix mode.")
if(!is.na(sum(grp))){  # Only analyze specified groups.
xx<-list()
for(i in 1:length(grp))xx[[i]]<-x[[grp[1]]]
x<-xx
}
J<-length(x)
tempn<-0
vecm<-0
for(j in 1:J){
temp<-x[[j]]
temp<-temp[!is.na(temp)] # Remove missing values.
tempn[j]<-length(temp)
x[[j]]<-temp
vecm[j]<-est(x[[j]],...)
}
xcen<-list()
flag<-rep(T,J)
for(j in 1:J){
flag[j]<-F
temp<-mean(vecm[flag])
xcen[[j]]<-x[[j]]-temp
flag[j]<-T
}
icrit<-round((1-alpha)*nboot)
bvec<-matrix(NA,nrow=J,ncol=nboot)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples. Please wait.")
for(j in 1:J){
paste("Working on group ",j)
data<-matrix(sample(xcen[[j]],size=length(x[[j]])*nboot,replace=TRUE),nrow=nboot)
bvec[j,]<-apply(data,1,est,...) # Bootstrapped values for jth group
}
vvec<-NA
for(j in 1:J){
vvec[j]<-sum((bvec[j,]-vecm[j])^2)/(nboot-1)
}
dis<-NA
for(i in 1:nboot){
dis[i]<-sum((bvec[,i]-vecm)^2/vvec)
}
tvec<-sum((0-vecm)^2/vvec)
dis<-sort(dis)
print(tvec)
print(dis[icrit])
print(vecm)
sig<-1-sum((tvec>=dis))/nboot
list(p.value=sig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.