wilcox_R: What the package does (short line)

bbwmcp <-
function(J, K, L, x, tr = 0.2, JKL = J * K*L, con = 0,
 alpha = 0.05, grp =c(1:JKL), nboot = 599, SEED = TRUE, ...)
{
        #
        # A bootstrap-t for multiple comparisons among
        # all main effects and interactions
#         for a between-by-between-within design.
        # 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.
        #
        #  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
}

        conM = con3way(J,K,L)
 p <- J*K*L
if(p>length(x))stop("JKL is less than the Number of groups")
JK=J*K
        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.
                data[[j]] = data[[j]] - mean(data[[j]], tr = tr)
        }
ilow=0-L
iup=0
for(j in 1:JK){
ilow <- ilow + L
 iup = iup + L
sel <- c(ilow:iup)
xx[sel]=listm(elimna(matl(xx[sel])))
 v[sel, sel] <- covmtrim(xx[sel], tr)
                }
A=lindep(xx,conM$conA,cmat=v,tr=tr)$test.stat
B=lindep(xx,conM$conB,cmat=v,tr=tr)$test.stat
C=lindep(xx,conM$conC,cmat=v,tr=tr)$test.stat
AB=lindep(xx,conM$conAB,cmat=v,tr=tr)$test.stat
AC=lindep(xx,conM$conAC,cmat=v,tr=tr)$test.stat
BC=lindep(xx,conM$conBC,cmat=v,tr=tr)$test.stat
ABC=lindep(xx,conM$conABC,cmat=v,tr=tr)$test.stat
        x <- data
        jp <- 1 - K
        kv <- 0
        if(SEED)
                set.seed(2)
        # set seed of random number generator so that
        #             results can be duplicated.
        testA = NA
        testB = NA
testC=NA
        testAB = NA
        testAC = NA
        testBC = NA
        testABC = NA
        bsam = list()
        bdat = list()
aboot=matrix(NA,nrow=nboot,ncol=ncol(conM$conA))
bboot=matrix(NA,nrow=nboot,ncol=ncol(conM$conB))
cboot=matrix(NA,nrow=nboot,ncol=ncol(conM$conC))
abboot=matrix(NA,nrow=nboot,ncol=ncol(conM$conAB))
acboot=matrix(NA,nrow=nboot,ncol=ncol(conM$conAC))
bcboot=matrix(NA,nrow=nboot,ncol=ncol(conM$conBC))
abcboot=matrix(NA,nrow=nboot,ncol=ncol(conM$conABC))
#        for(j in 1:JK)
#                nvec[j] = length(x[[j]])
        for(ib in 1:nboot) {
                ilow <- 1 - L
                iup = 0
 for(j in 1:JK) {
 ilow <- ilow + L
 iup = iup + L
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-L
iup=0
for(j in 1:JK){
ilow <- ilow + L
 iup = iup + L
sel <- c(ilow:iup)
 v[sel, sel] <- covmtrim(bsam[sel], tr)
                }
temp=abs(lindep(bsam,conM$conA, cmat=v,tr=tr)$test.stat[,4])
aboot[ib,]=temp
testA[ib] = max(temp)
temp=abs(lindep(bsam,conM$conB,cmat=v,tr=tr)$test.stat[,4])
bboot[ib,]=temp
testB[ib]= max(temp)

temp=abs(lindep(bsam,conM$conC,cmat=v,tr=tr)$test.stat[,4])
cboot[ib,]=temp
testC[ib]= max(temp)

temp=abs(lindep(bsam,conM$conAC,cmat=v,tr=tr)$test.stat[,4])
acboot[ib,]=temp
testAC[ib]= max(temp)

temp=abs(lindep(bsam,conM$conBC,cmat=v,tr=tr)$test.stat[,4])
bcboot[ib,]=temp
testBC[ib]= max(temp)

temp=abs(lindep(bsam,conM$conAB,cmat=v,tr=tr)$test.stat[,4])
testAB[ib] = max(temp)
abboot[ib,]=temp

temp=abs(lindep(bsam,conM$conABC,cmat=v,tr=tr)$test.stat[,4])
abcboot[ib,]=temp
testABC[ib]= max(temp)

        }
pbA=NA
pbB=NA
pbC=NA
pbAB=NA
pbAC=NA
pbBC=NA
pbABC=NA
for(j in 1:ncol(aboot))pbA[j]=mean((abs(A[j,4])<aboot[,j]))
for(j in 1:ncol(bboot))pbB[j]=mean((abs(B[j,4])<bboot[,j]))
for(j in 1:ncol(cboot))pbC[j]=mean((abs(C[j,4])<cboot[,j]))
for(j in 1:ncol(abboot))pbAB[j]=mean((abs(AB[j,4])<abboot[,j]))
for(j in 1:ncol(acboot))pbAC[j]=mean((abs(AC[j,4])<acboot[,j]))
for(j in 1:ncol(bcboot))pbBC[j]=mean((abs(BC[j,4])<bcboot[,j]))
for(j in 1:ncol(abcboot))pbABC[j]=mean((abs(ABC[j,4])<abcboot[,j]))
        critA = sort(testA)
        critB = sort(testB)
        critC = sort(testC)
        critAB = sort(testAB)
        critAC = sort(testAC)
        critBC = sort(testBC)
        critABC = sort(testABC)
        ic <- floor((1 - alpha) * nboot)
        critA = critA[ic]
        critB = critB[ic]
        critC = critC[ic]
        critAB = critAB[ic]
        critAC = critAC[ic]
        critBC = critBC[ic]
        critABC = critABC[ic]
critA=matrix(critA,ncol=1,nrow=nrow(A))
dimnames(critA)=list(NULL,c("crit.val"))
p.value=pbA
A=cbind(A,critA,p.value)

critB=matrix(critB,ncol=1,nrow=nrow(B))
dimnames(critB)=list(NULL,c("crit.val"))
p.value=pbB
B=cbind(B,critB,p.value)

critC=matrix(critC,ncol=1,nrow=nrow(C))
dimnames(critC)=list(NULL,c("crit.val"))
p.value=pbC
C=cbind(C,critC,p.value)

critAB=matrix(critAB,ncol=1,nrow=nrow(AB))
dimnames(critAB)=list(NULL,c("crit.val"))
p.value=pbAB
AB=cbind(AB,critAB,p.value)

critAC=matrix(critAC,ncol=1,nrow=nrow(AC))
dimnames(critAC)=list(NULL,c("crit.val"))
p.value=pbAC
AC=cbind(AC,critAC,p.value)


critBC=matrix(critBC,ncol=1,nrow=nrow(BC))
dimnames(critBC)=list(NULL,c("crit.val"))
p.value=pbBC
BC=cbind(BC,critBC,p.value)

critABC=matrix(critABC,ncol=1,nrow=nrow(ABC))
dimnames(critABC)=list(NULL,c("crit.val"))
p.value=pbABC
ABC=cbind(ABC,critABC,p.value)

list(Fac.A=A,Fac.B=B,Fac.C=C,Fac.AB=AB,Fac.AC=AC,Fac.BC=BC,Fac.ABC=ABC)
}