mcppb:

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

Usage

mcppb(x, crit = NA, con = 0, tr = 0.2, alpha = 0.05, nboot = 2000, grp = NA, WIN = FALSE, win = 0.1)

Arguments

x
crit
con
tr
alpha
nboot
grp
WIN
win

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, crit = NA, con = 0, tr = 0.2, alpha = 0.05, nboot = 2000, 
    grp = NA, WIN = FALSE, win = 0.1) 
{
    con <- as.matrix(con)
    if (is.matrix(x)) {
        xx <- list()
        for (i in 1:ncol(x)) {
            xx[[i]] <- x[, i]
        }
        x <- xx
    }
    if (!is.list(x)) 
        stop("Data must be stored in list mode or in matrix mode.")
    if (!is.na(sum(grp))) {
        xx <- list()
        for (i in 1:length(grp)) xx[[i]] <- x[[grp[1]]]
        x <- xx
    }
    J <- length(x)
    tempn <- 0
    for (j in 1:J) {
        temp <- x[[j]]
        temp <- temp[!is.na(temp)]
        tempn[j] <- length(temp)
        x[[j]] <- temp
    }
    Jm <- J - 1
    d <- ifelse(sum(con^2) == 0, (J^2 - J)/2, ncol(con))
    if (is.na(crit) && tr != 0.2) 
        stop("A critical value must be specified when\nthe amount of trimming differs from .2")
    if (WIN) {
        if (tr < 0.2) 
            warning("When Winsorizing, the amount of trimming should be at least\n.2")
        if (win > tr) 
            stop("Amount of Winsorizing must <= amount of trimming")
        if (min(tempn) < 15) {
            warning("Winsorizing with sample sizes less than 15 can")
            warning(" result in poor control over the probability of a Type I error")
        }
        for (j in 1:J) {
            x[[j]] <- winval(x[[j]], win)
        }
    }
    if (is.na(crit)) {
        if (d == 1) 
            crit <- alpha/2
        if (d == 2 && alpha == 0.05 && nboot == 1000) 
            crit <- 0.014
        if (d == 2 && alpha == 0.05 && nboot == 2000) 
            crit <- 0.014
        if (d == 3 && alpha == 0.05 && nboot == 1000) 
            crit <- 0.009
        if (d == 3 && alpha == 0.05 && nboot == 2000) 
            crit <- 0.0085
        if (d == 3 && alpha == 0.025 && nboot == 1000) 
            crit <- 0.004
        if (d == 3 && alpha == 0.025 && nboot == 2000) 
            crit <- 0.004
        if (d == 3 && alpha == 0.01 && nboot == 1000) 
            crit <- 0.001
        if (d == 3 && alpha == 0.01 && nboot == 2000) 
            crit <- 0.001
        if (d == 4 && alpha == 0.05 && nboot == 2000) 
            crit <- 0.007
        if (d == 5 && alpha == 0.05 && nboot == 2000) 
            crit <- 0.006
        if (d == 6 && alpha == 0.05 && nboot == 1000) 
            crit <- 0.004
        if (d == 6 && alpha == 0.05 && nboot == 2000) 
            crit <- 0.0045
        if (d == 6 && alpha == 0.025 && nboot == 1000) 
            crit <- 0.002
        if (d == 6 && alpha == 0.025 && nboot == 2000) 
            crit <- 0.0015
        if (d == 6 && alpha == 0.01 && nboot == 2000) 
            crit <- 5e-04
        if (d == 10 && alpha == 0.05 && nboot <= 2000) 
            crit <- 0.002
        if (d == 10 && alpha == 0.05 && nboot == 3000) 
            crit <- 0.0023
        if (d == 10 && alpha == 0.025 && nboot <= 2000) 
            crit <- 5e-04
        if (d == 10 && alpha == 0.025 && nboot == 3000) 
            crit <- 0.001
        if (d == 15 && alpha == 0.05 && nboot == 2000) 
            crit <- 0.0016
        if (d == 15 && alpha == 0.025 && nboot == 2000) 
            crit <- 5e-04
        if (d == 15 && alpha == 0.05 && nboot == 5000) 
            crit <- 0.0026
        if (d == 15 && alpha == 0.025 && nboot == 5000) 
            crit <- 6e-04
    }
    if (is.na(crit) && alpha == 0.05) 
        crit <- 0.0268660714 * (1/d) - 0.0003321429
    if (is.na(crit)) 
        crit <- alpha/(2 * d)
    if (d > 10 && nboot < 5000) 
        warning("Suggest using nboot=5000 when the number\nof contrasts exceeds 10.")
    icl <- round(crit * nboot) + 1
    icu <- round((1 - crit) * nboot)
    if (sum(con^2) == 0) {
        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
            }
        }
    }
    psihat <- matrix(0, ncol(con), 6)
    dimnames(psihat) <- list(NULL, c("con.num", "psihat", "se", 
        "ci.lower", "ci.upper", "p.value"))
    if (nrow(con) != length(x)) 
        stop("The number of groups does not match the number\n of  contrast coefficients.")
    bvec <- matrix(NA, nrow = J, ncol = nboot)
    set.seed(2)
    print("Taking bootstrap samples. Please wait.")
    for (j in 1:J) {
        print(paste("Working on group ", j))
        data <- matrix(sample(x[[j]], size = length(x[[j]]) * 
            nboot, replace = TRUE), nrow = nboot)
        bvec[j, ] <- apply(data, 1, mean, tr)
    }
    test <- NA
    for (d in 1:ncol(con)) {
        top <- 0
        for (i in 1:J) {
            top <- top + con[i, d] * bvec[i, ]
        }
        test[d] <- sum((top > 0))/nboot
        test[d] <- min(test[d], 1 - test[d])
        top <- sort(top)
        psihat[d, 4] <- top[icl]
        psihat[d, 5] <- top[icu]
    }
    for (d in 1:ncol(con)) {
        psihat[d, 1] <- d
        testit <- lincon(x, con[, d], tr, pr = FALSE)
        psihat[d, 6] <- test[d]
        psihat[d, 2] <- testit$psihat[1, 2]
        psihat[d, 3] <- testit$test[1, 4]
    }
    print("Reminder: To control FWE, reject if the p-value is less than")
    print("the crit.p.value listed in the output.")
    list(psihat = psihat, crit.p.value = crit, con = con)
  }

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