R/bbwtrimbt.R

In musto101/wilcox_R: What the package does (short line)

bbwtrimbt <-
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, ...)
{
        #
        # Bootstrap-t for omniubs tests associated with
        # all main effects and interactions
        # for a between-by-between-within design.
        #
        #  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 for all three factors: levels 1,1,1.
        #  x[[2]] is assumed to contain the data for level 1 of the
        #  first factor and second factor and level 2 of the third: level 1,1,2
        #  x[[K]] is the data for level 1,1,L
        #  x[[L+1]] is the data for level 1,2,1, x[[2K]] is level 1,2,2, 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.data.frame(x))data=as.matrix(x)
        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)
#                }
test.stat=bbwtrim(J,K,L,xx,tr=tr)
        x <- data  # Centered 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=NA
bboot=NA
cboot=NA
abboot=NA
acboot=NA
bcboot=NA
abcboot=NA
nvec=NA
        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]]]
}
}
temp=bbwtrim(J,K,L,bsam,tr=tr)
aboot[ib]=temp$Qa
bboot[ib]=temp$Qb
cboot[ib]=temp$Qc
acboot[ib]=temp$Qac
bcboot[ib]=temp$Qbc
abboot[ib]=temp$Qab
abcboot[ib]=temp$Qabc
}}
pbA=NA
pbB=NA
pbC=NA
pbAB=NA
pbAC=NA
pbBC=NA
pbABC=NA
pbA=mean(test.stat$Qa[1,1]<aboot)
pbB=mean(test.stat$Qb[1,1]<bboot)
pbC=mean(test.stat$Qc[1,1]<cboot)
pbAB=mean(test.stat$Qab[1,1]<abboot)
pbAC=mean(test.stat$Qac[1,1]<acboot)
pbBC=mean(test.stat$Qbc[1,1]<bcboot)
pbABC=mean(test.stat$Qabc[1,1]<abcboot)
list(p.value.A=pbA,p.value.B=pbB,p.value.C=pbC,p.value.AB=pbAB,
p.value.AC=pbAC,p.value.BC=pbBC,p.value.ABC=pbABC)
}