discstep:

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

Usage

discstep(x, nboot = 500, alpha = 0.05, SEED = TRUE)

Arguments

x
nboot
alpha
SEED

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 (x, nboot = 500, alpha = 0.05, SEED = TRUE) 
{
    if (is.matrix(x) || is.data.frame(x)) 
        x = listm(x)
    library(mc2d)
    if (SEED) 
        set.seed(2)
    vals = lapply(x, unique)
    vals = sort(elimna(list2vec(vals)))
    K = length(unique(vals))
    n = lapply(x, length)
    n = list2vec(n)
    J = length(x)
    if (J == 2) 
        stop("For 2 groups use disc2com")
    if (J > 5) 
        stop("Designed for 5 groups or less")
    com = modgen(J)
    ntest = length(com)
    jp1 = J + 1
    com = com[jp1:length(com)]
    ntest = length(com)
    mout = matrix(NA, nrow = ntest, ncol = 3)
    dimnames(mout) = list(NULL, c("Groups", "p-value", "p.crit"))
    test = NULL
    for (i in 1:ntest) {
        test[i] = discANOVA.sub(x[com[[i]]])$test
        nmod = length(com[[i]]) - 1
        temp = c(nmod:0)
        mout[i, 1] = sum(com[[i]] * 10^temp)
    }
    mout[, 3] = alpha
    xx = list()
    pv = NA
    jm2 = J - 2
    mout[, 3] = alpha
    TB = matrix(NA, nrow = nboot, ncol = ntest)
    step1 = discANOVA.sub(x)
    C1 = step1$C1
    HT = NULL
    for (i in 1:K) HT[i] = mean(C1[i, ])
    for (ib in 1:nboot) {
        xx = list()
        for (j in 1:J) {
            temp = rmultinomial(n[j], 1, HT)
            xx[[j]] = which(temp[1, ] == 1)
            for (i in 2:n[j]) xx[[j]][i] = which(temp[i, ] == 
                1)
        }
        for (k in 1:ntest) TB[ib, k] = discANOVA.sub(xx[com[[k]]])$test
    }
    for (k in 1:ntest) {
        mout[k, 2] = 1 - mean(test[k] > TB[, k]) - 0.5 * mean(test[k] == 
            TB[, k])
        pnum = length(com[[k]])
        pe = 1 - (1 - alpha)^(pnum/J)
        if (length(com[[k]]) <= jm2) 
            mout[k, 3] = pe
    }
    list(results = mout[nrow(mout):1, ])
  }

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