rmmcppbd:

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

Usage

rmmcppbd(x, y = NULL, alpha = 0.05, con = 0, est = onestep, plotit = TRUE, grp = NA, nboot = NA, hoch = TRUE, SEED = TRUE, ...)

Arguments

x
y
alpha
con
est
plotit
grp
nboot
hoch
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, y = NULL, alpha = 0.05, con = 0, est = onestep, 
    plotit = TRUE, grp = NA, nboot = NA, hoch = TRUE, SEED = TRUE, 
    ...) 
{
    if (!is.null(y[1])) 
        x <- cbind(x, y)
    if (!is.list(x) && !is.matrix(x)) 
        stop("Data must be stored in a matrix or in list mode.")
    if (is.list(x)) {
        if (is.matrix(con)) {
            if (length(x) != nrow(con)) 
                stop("The number of rows in con is not equal to the number of groups.")
        }
    }
    if (is.list(x)) {
        mat <- matl(x)
    }
    if (is.matrix(x) && is.matrix(con)) {
        if (ncol(x) != nrow(con)) 
            stop("The number of rows in con is not equal to the number of groups.")
        mat <- x
    }
    if (is.matrix(x)) 
        mat <- x
    if (!is.na(sum(grp))) 
        mat <- mat[, grp]
    x <- mat
    mat <- elimna(mat)
    x <- mat
    J <- ncol(mat)
    n = nrow(mat)
    if (n >= 80) 
        hoch = T
    Jm <- J - 1
    if (sum(con^2) == 0) {
        d <- (J^2 - J)/2
        con <- matrix(0, J, d)
        id <- 0
        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
            }
        }
    }
    d <- ncol(con)
    if (is.na(nboot)) {
        nboot <- 5000
        if (d <= 10) 
            nboot <- 3000
        if (d <= 6) 
            nboot <- 2000
        if (d <= 4) 
            nboot <- 1000
    }
    n <- nrow(mat)
    crit.vec <- alpha/c(1:d)
    connum <- ncol(con)
    xx <- x %*% con
    xx <- as.matrix(xx)
    if (SEED) 
        set.seed(2)
    psihat <- matrix(0, connum, nboot)
    bvec <- matrix(NA, ncol = connum, nrow = nboot)
    data <- matrix(sample(n, size = n * nboot, replace = TRUE), 
        nrow = nboot)
    if (ncol(xx) == 1) {
        for (ib in 1:nboot) psihat[1, ib] <- est(xx[data[ib, 
            ]], ...)
    }
    if (ncol(xx) > 1) {
        for (ib in 1:nboot) psihat[, ib] <- apply(xx[data[ib, 
            ], ], 2, est, ...)
    }
    test <- 1
    for (ic in 1:connum) {
        test[ic] <- (sum(psihat[ic, ] > 0) + 0.5 * sum(psihat[ic, 
            ] == 0))/nboot
        test[ic] <- min(test[ic], 1 - test[ic])
    }
    test <- 2 * test
    ncon <- ncol(con)
    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[2] <- alpha/2
    }
    if (hoch) 
        dvec <- alpha/(2 * c(1:ncon))
    dvec <- 2 * dvec
    if (plotit && connum == 1) {
        plot(c(psihat[1, ], 0), xlab = "", ylab = "Est. Difference")
        points(psihat[1, ])
        abline(0, 0)
    }
    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 <- apply(xx, 2, est, ...)
    psi <- 1
    icl <- round(dvec[ncon] * nboot/2) + 1
    icu <- nboot - icl - 1
    for (ic in 1:ncol(con)) {
        output[ic, 2] <- tmeans[ic]
        output[ic, 1] <- ic
        output[ic, 3] <- test[ic]
        output[temp2, 4] <- zvec
        temp <- sort(psihat[ic, ])
        output[ic, 5] <- temp[icl]
        output[ic, 6] <- temp[icu]
    }
    num.sig <- sum(output[, 3] <= output[, 4])
    list(output = output, con = con, num.sig = num.sig)
  }

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