R/bwwmcppb.sub.R

bwwmcppb.sub <-
function(J, K,L, x, est=tmean, JKL = J * K*L, con = 0,
 alpha = 0.05, grp =c(1:JKL), nboot = 500, bhop=FALSE,SEED = TRUE, ...){
        #
        # A percentile bootstrap for multiple comparisons 
        # for all main effects and interactions.
        # The analysis is done by generating bootstrap samples and
        # using an appropriate linear contrast.
        #
        #  The R variable x is assumed to contain the raw
        #  data stored in list mode or in a matrix.
        #  If in list mode, x[[1]] contains the data
        #  for the first level of both factors: level 1,1.
        #  x[[2]] is assumed to contain the data for level 1 of the
        #  first factor and level 2 of the second: level 1,2
        #  x[[K]] is the data for level 1,K
        #  x[[K+1]] is the data for level 2,1, x[[2K]] is level 2,K, etc.
        #
#
#  J independent groups, KL dependent groups
#

       #  If the data are in a matrix, column 1 is assumed to
        #  correspond to x[[1]], column 2 to x[[2]], etc.
        #
        #  When in list mode x is assumed to have length JK, the total number
        #  groups being tested, but a subset of the data can be analyzed
        #  using grp
        #
        if(is.matrix(x)) {
                y <- list()
                for(j in 1:ncol(x))
                        y[[j]] <- x[, j]
x=y
}
#        nvec <- NA
#for(j in 1:length(x))nvec[j]=length(x[[j]])
ncon=ncol(con)
 p <- J*K*L
if(p>length(x))stop("JKL is less than the Number of groups")
JK=J*K
KL=K*L
#        v <- matrix(0, p, p)
        data <- list()
xx=list()
        for(j in 1:length(x)) {
#                data[[j]] <- x[[grp[j]]]
xx[[j]]=x[[grp[j]]] # save input data
#                # Now have the groups in proper order.
        }
ilow=1-KL
iup=0
for(j in 1:J){
ilow <- ilow + KL
 iup = iup + KL
sel <- c(ilow:iup)
xx[sel]=listm(elimna(matl(xx[sel])))
}

        jp <- 1 - KL
        kv <- 0
        if(SEED)
                set.seed(2)
        # set seed of random number generator so that
        #             results can be duplicated.
        # Next determine the n_j values
        testA = NA
        bsam = list()
        bdat = list()
aboot=matrix(NA,nrow=nboot,ncol=ncol(con))
tvec=NA
tvec=linhat(x,con,est=est,...)
        for(ib in 1:nboot) {
                ilow <- 1 - KL
                iup = 0
 for(j in 1:J) {
 ilow <- ilow + KL
 iup = iup + KL
nv=length(x[[ilow]])
bdat[[j]] = sample(nv, size = nv, replace =T)
for(k in ilow:iup){
 bsam[[k]] = x[[k]][bdat[[j]]]
}
}
ilow=0-KL
iup=0
aboot[ib,]=linhat(bsam,con=con,est=est,...)
}
pbA=NA
for(j in 1:ncol(aboot)){
pbA[j]=mean(aboot[,j]>0)
pbA[j]=2*min(c(pbA[j],1-pbA[j]))
}
# Determine critical values
if(!bhop){
if(alpha==.05){
dvec<-c(.05,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
if(ncol(con) > 10){
avec<-.05/c(11:(ncol(con)))
dvec<-c(dvec,avec)
}}
if(alpha==.01){
dvec<-c(.01,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
if(con > 10){
avec<-.01/c(11:ncol(con))
dvec<-c(dvec,avec)
}}
if(alpha != .05 && alpha != .01){
dvec<-alpha/c(1:ncol(con))
}
}
if(bhop)dvec<-(ncol(con)-c(1:ncol(con))+1)*alpha/ncol(con)
outputA<-matrix(0,ncol(con),6)
dimnames(outputA)<-list(NULL,c("con.num","psihat","p.value","p.crit",
"ci.lower","ci.upper"))
test=pbA
temp2<-order(0-test)
zvec<-dvec[1:ncon]
sigvec<-(test[temp2]>=zvec)
outputA[temp2,4]<-zvec
icl<-round(dvec[ncon]*nboot/2)+1
icu<-nboot-icl-1
outputA[,2]<-tvec
for (ic in 1:ncol(con)){
outputA[ic,1]<-ic
outputA[ic,3]<-test[ic]
temp<-sort(aboot[,ic])
outputA[ic,5]<-temp[icl]
outputA[ic,6]<-temp[icu]
}
outputA
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.