sppbi:

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

Usage

sppbi(J, K, x, est = tmean, JK = J * K, grp = c(1:JK), nboot = 500, SEED = TRUE, pr = TRUE, ...)

Arguments

J
K
x
est
JK
grp
nboot
SEED
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), nboot = 500, 
    SEED = TRUE, pr = TRUE, ...) 
{
    if (pr) 
        print("As of Oct. 2014, 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
    MJ <- (J^2 - J)/2
    MK <- (K^2 - K)/2
    JMK <- J * MK
    Jm <- J - 1
    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, ncol = J, nrow = nboot)
    print("Taking bootstrap samples. Please wait.")
    mvec <- NA
    it <- 0
    for (j in 1:J) {
        paste("Working on level ", j, " of Factor A")
        x <- matrix(NA, nrow = nvec[j], ncol = MK)
        im <- 0
        for (k in 1:K) {
            for (kk in 1:K) {
                if (k < kk) {
                  im <- im + 1
                  kp <- j * K + k - K
                  kpp <- j * K + kk - K
                  x[, im] <- xx[[kp]] - xx[[kpp]]
                  it <- it + 1
                  mvec[it] <- est(x[, im], ...)
                }
            }
        }
        data <- matrix(sample(nvec[j], size = nvec[j] * nboot, 
            replace = TRUE), nrow = nboot)
        bvec <- matrix(NA, ncol = MK, nrow = nboot)
        mat <- listm(x)
        for (k in 1:MK) {
            temp <- x[, k]
            bvec[, k] <- apply(data, 1, rmanogsub, temp, est, 
                ...)
        }
        if (j == 1) 
            bloc <- bvec
        if (j > 1) 
            bloc <- cbind(bloc, bvec)
    }
    MJMK <- MJ * MK
    con <- matrix(0, nrow = JMK, ncol = MJMK)
    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 = MK - 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 (MK > 1) {
        for (k in 2:MK) {
            con1 <- push(con1)
            con <- cbind(con, con1)
        }
    }
    bcon <- t(con) %*% t(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
    smat <- var(temp)
    if (sum(is.na(smat)) == 0) {
        chkrank <- qr(smat)$rank
        bcon <- rbind(bcon, vecz)
        if (chkrank == ncol(smat)) 
            dv <- mahalanobis(bcon, tvec, smat)
        if (chkrank < ncol(smat)) {
            smat <- ginv(smat)
            dv <- mahalanobis(bcon, tvec, smat, inverted = T)
        }
    }
    if (sum(is.na(smat)) > 0) 
        print("Computational Problem. Try est=tmean or use function spmcpi or tsplitbt")
    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.