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

Usage Arguments Examples

1	spmcpi(J, K, x, est = tmean, JK = J * K, grp = c(1:JK), alpha = 0.05, nboot = NA, SEED = TRUE, pr = TRUE, SR = FALSE, ...)

`J`
`K`
`x`
`est`
`JK`
`grp`
`alpha`
`nboot`
`SEED`
`pr`
`SR`
`...`

##---- 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), alpha = 0.05, 
    nboot = NA, SEED = TRUE, pr = TRUE, SR = FALSE, ...) 
{
    if (is.matrix(x)) {
        y <- list()
        for (j in 1:ncol(x)) y[[j]] <- x[, j]
        x <- y
    }
    if (pr) 
        print("As of Sept. 2005, est defaults to tmean")
    JK <- J * K
    if (JK != length(x)) {
        print("Something is wrong.")
        paste(" Expected ", JK, " groups but x contains ", length(x), 
            "groups instead.")
        stop()
    }
    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]])
    }
    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)
        }
    }
    d <- ncol(con)
    if (is.na(nboot)) {
        if (d <= 4) 
            nboot <- 1000
        if (d > 4) 
            nboot <- 5000
    }
    connum <- d
    psihat <- matrix(0, connum, nboot)
    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)
    }
    test <- 1
    for (ic in 1:connum) {
        psihat[ic, ] <- apply(bloc, 1, bptdpsi, con[, ic])
        test[ic] <- (sum(psihat[ic, ] > 0) + 0.5 * sum(psihat[ic, 
            ] == 0))/nboot
        test[ic] <- min(test[ic], 1 - test[ic])
    }
    ncon <- ncol(con)
    dvec <- alpha/c(1:ncon)
    if (SR) {
        okay = FALSE
        if (identical(est, onestep)) 
            okay = TRUE
        if (identical(est, mom)) 
            okay = TRUE
        if (!okay) 
            stop("For estimators other than onestep and mom, use SR=FALSE")
        if (alpha == 0.05) {
            dvec <- c(0.025, 0.025, 0.0169, 0.0127, 0.0102, 0.00851, 
                0.0073, 0.00639, 0.00568, 0.00511)
            if (ncon > 10) {
                avec <- 0.05/c(11:ncon)
                dvec <- c(dvec, avec)
            }
        }
        if (alpha == 0.01) {
            dvec <- c(0.005, 0.005, 0.00334, 0.00251, 0.00201, 
                0.00167, 0.00143, 0.00126, 0.00112, 0.00101)
            if (ncon > 10) {
                avec <- 0.01/c(11:ncon)
                dvec <- c(dvec, avec)
            }
        }
        if (alpha != 0.05 && alpha != 0.01) {
            dvec <- alpha/c(1:ncon)
            dvec[1] <- alpha/2
        }
    }
    temp2 <- order(0 - test)
    ncon <- ncol(con)
    zvec <- dvec[1:ncon]
    sigvec <- (test[temp2] >= zvec)
    output <- matrix(0, connum, 6)
    dimnames(output) <- list(NULL, c("con.num", "psihat", "p.value", 
        "p.crit", "ci.lower", "ci.upper"))
    tmeans <- mvec
    psi <- 1
    for (ic in 1:ncol(con)) {
        output[ic, 2] <- sum(con[, ic] * tmeans)
        output[ic, 1] <- ic
        output[ic, 3] <- test[ic]
        output[temp2, 4] <- zvec
        temp <- sort(psihat[ic, ])
        icl <- round(dvec[ncon] * nboot) + 1
        icu <- nboot - (icl - 1)
        output[ic, 5] <- temp[icl]
        output[ic, 6] <- temp[icu]
    }
    output[, 3] <- 2 * output[, 3]
    if (SR) 
        output[, 4] <- 2 * output[, 4]
    num.sig <- sum(output[, 3] <= output[, 4])
    list(output = output, con = con, num.sig = num.sig)
  }