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

Usage Arguments Examples

1	mulrank(J, K, x, grp = c(1:p), p = J * K)

`J`
`K`
`x`
`grp`
`p`

##---- 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, grp = c(1:p), p = J * K) 
{
    newx = list()
    GV = matrix(c(1:p), ncol = K, byrow = TRUE)
    if (is.list(x)) {
        temp = NA
        jk = 0
        for (j in 1:J) {
            temp = elimna(matl(x[GV[j, ]]))
            for (k in 1:K) {
                jk = jk + 1
                newx[[jk]] = temp[, k]
            }
        }
        x = NA
        x = newx
    }
    if (is.matrix(x)) {
        x = elimna(x)
        x <- listm(x)
    }
    xx <- list()
    nvec <- NA
    for (j in 1:p) {
        xx[[j]] <- x[[grp[j]]]
        nvec[j] <- length(xx[[j]])
    }
    Nrow = nvec[GV[, 1]]
    v <- matrix(0, p, p)
    Ja <- matrix(1, J, J)
    Ia <- diag(1, J)
    Pa <- Ia - Ja/J
    Jb <- matrix(1, K, K)
    Ib <- diag(1, K)
    Pb <- Ib - Jb/K
    cona <- kron(Pa, Ib)
    xr <- list()
    N <- 0
    jj = 0
    for (k in 1:K) {
        temp <- x[[k]]
        jk <- k
        for (j in 2:J) {
            jj = jj + 1
            jk <- jk + K
            temp <- c(temp, x[[jk]])
        }
        N <- length(temp)
        pr <- rank(temp)
        xr[[k]] <- pr[1:nvec[k]]
        top <- nvec[k]
        jk <- k
        bot <- 1
        for (j in 2:J) {
            jk <- jk + K
            bot <- bot + nvec[jk]
            top <- top + nvec[jk]
            xr[[jk]] <- pr[bot:top]
        }
    }
    phat <- NA
    botk <- 0
    for (j in 1:J) {
        for (k in 1:K) {
            botk <- botk + 1
            phat[botk] <- (mean(xr[[botk]]) - 0.5)/N
        }
    }
    klow <- 1 - K
    kup <- 0
    for (j in 1:J) {
        klow <- klow + K
        kup <- kup + K
        sel <- c(klow:kup)
        v[sel, sel] <- covmtrim(xr[klow:kup], tr = 0)/N
    }
    qhat <- matrix(phat, J, K, byrow = TRUE)
    test <- N * t(phat) %*% cona %*% phat/sum(diag(cona %*% v))
    nu1 <- sum(diag(cona %*% v))^2/sum(diag(cona %*% v %*% cona %*% 
        v))
    sig.level <- 1 - pf(test, nu1, 1e+06)
    list(test.stat = test[1, 1], nu1 = nu1, p.value = sig.level, 
        N = N, q.hat = qhat)
  }