sppba:

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

Usage

sppba(J, K, x, est = tmean, JK = J * K, grp = c(1:JK), avg = TRUE, nboot = 500, SEED = TRUE, MC = FALSE, MDIS = FALSE, pr = TRUE, ...)

Arguments

J
K
x
est
JK
grp
avg
nboot
SEED
MC
MDIS
pr
...

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (J, K, x, est = tmean, JK = J * K, grp = c(1:JK), avg = TRUE, 
    nboot = 500, SEED = TRUE, MC = FALSE, MDIS = FALSE, pr = TRUE, 
    ...) 
{
    if (pr) 
        print("As of Oct. 2014 the argument est defaults to tmean")
    library(MASS)
    if (is.matrix(x)) {
        y <- list()
        for (j in 1:ncol(x)) y[[j]] <- x[, j]
        x <- y
    }
    JK <- J * K
    data <- list()
    for (j in 1:length(x)) {
        data[[j]] <- x[[grp[j]]]
    }
    x <- data
    jp <- 1 - K
    kv <- 0
    kv2 <- 0
    for (j in 1:J) {
        jp <- jp + K
        xmat <- matrix(NA, ncol = K, nrow = length(x[[jp]]))
        for (k in 1:K) {
            kv <- kv + 1
            xmat[, k] <- x[[kv]]
        }
        xmat <- elimna(xmat)
        for (k in 1:K) {
            kv2 <- kv2 + 1
            x[[kv2]] <- xmat[, k]
        }
    }
    xx <- x
    if (SEED) 
        set.seed(2)
    nvec <- NA
    jp <- 1 - K
    for (j in 1:J) {
        jp <- jp + K
        nvec[j] <- length(x[[jp]])
    }
    bloc <- matrix(NA, nrow = J, ncol = nboot)
    print("Taking bootstrap samples. Please wait.")
    mvec <- NA
    ik <- 0
    for (j in 1:J) {
        paste("Working on level ", j, " of Factor A")
        x <- matrix(NA, nrow = nvec[j], ncol = K)
        for (k in 1:K) {
            ik <- ik + 1
            x[, k] <- xx[[ik]]
            if (!avg) 
                mvec[ik] <- est(xx[[ik]], ...)
        }
        tempv <- apply(x, 2, est, ...)
        data <- matrix(sample(nvec[j], size = nvec[j] * nboot, 
            replace = TRUE), nrow = nboot)
        bvec <- matrix(NA, ncol = K, nrow = nboot)
        for (k in 1:K) {
            temp <- x[, k]
            bvec[, k] <- apply(data, 1, rmanogsub, temp, est, 
                ...)
        }
        if (avg) {
            mvec[j] <- mean(tempv)
            bloc[j, ] <- apply(bvec, 1, mean)
        }
        if (!avg) {
            if (j == 1) 
                bloc <- bvec
            if (j > 1) 
                bloc <- cbind(bloc, bvec)
        }
    }
    if (avg) {
        d <- (J^2 - J)/2
        con <- matrix(0, J, d)
        id <- 0
        Jm <- J - 1
        for (j in 1:Jm) {
            jp <- j + 1
            for (k in jp:J) {
                id <- id + 1
                con[j, id] <- 1
                con[k, id] <- 0 - 1
            }
        }
    }
    if (!avg) {
        MJK <- K * (J^2 - J)/2
        JK <- J * K
        MJ <- (J^2 - J)/2
        cont <- matrix(0, nrow = J, ncol = MJ)
        ic <- 0
        for (j in 1:J) {
            for (jj in 1:J) {
                if (j < jj) {
                  ic <- ic + 1
                  cont[j, ic] <- 1
                  cont[jj, ic] <- 0 - 1
                }
            }
        }
        tempv <- matrix(0, nrow = K - 1, ncol = MJ)
        con1 <- rbind(cont[1, ], tempv)
        for (j in 2:J) {
            con2 <- rbind(cont[j, ], tempv)
            con1 <- rbind(con1, con2)
        }
        con <- con1
        if (K > 1) {
            for (k in 2:K) {
                con1 <- push(con1)
                con <- cbind(con, con1)
            }
        }
    }
    if (!avg) 
        bcon <- t(con) %*% t(bloc)
    if (avg) 
        bcon <- t(con) %*% (bloc)
    tvec <- t(con) %*% mvec
    tvec <- tvec[, 1]
    tempcen <- apply(bcon, 1, mean)
    vecz <- rep(0, ncol(con))
    bcon <- t(bcon)
    temp = bcon
    for (ib in 1:nrow(temp)) temp[ib, ] = temp[ib, ] - tempcen + 
        tvec
    bcon <- rbind(bcon, vecz)
    if (!MDIS) {
        if (!MC) 
            dv = pdis(bcon, center = tvec, na.rm = FALSE)
        if (MC) 
            dv = pdisMC(bcon, center = tvec, na.rm = FALSE)
        lbcon = length(elimna(bcon))
        bplus <- nboot + 1
        if (lbcon < bplus) {
            print(paste("Effective value for nboot is", lbcon - 
                1))
            nboot = lbcon - 1
        }
    }
    if (MDIS) {
        smat <- var(temp)
        bcon <- rbind(bcon, vecz)
        chkrank <- qr(smat)$rank
        if (chkrank == ncol(smat)) 
            dv <- mahalanobis(bcon, tvec, smat)
        if (chkrank < ncol(smat)) {
            smat <- ginv(smat)
            dv <- mahalanobis(bcon, tvec, smat, inverted = TRUE)
        }
    }
    bplus <- nboot + 1
    sig.level <- 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
    list(p.value = sig.level, psihat = tvec, con = con)
  }

musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.