mcppb20:

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

Usage

mcppb20(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) {
        print("A critical value must be specified when")
        stop("the amount of trimming differs from .2")
    }
    if (WIN) {
        if (tr < 0.2) {
            print("Warning: When Winsorizing, the amount")
            print("of trimming should be at least .2")
        }
        if (win > tr) 
            stop("Amount of Winsorizing must <= amount of trimming")
        if (min(tempn) < 15) {
            print("Warning: Winsorizing with sample sizes")
            print("less than 15 can result in poor control")
            print("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) {
        print("Warning: Suggest using nboot=5000")
        print("when the number of 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)) {
        print("The number of groups does not match")
        stop("the number of contrast coefficients.")
    }
    bvec <- matrix(NA, nrow = J, ncol = nboot)
    set.seed(2)
    print("Taking bootstrap samples. Please wait.")
    for (j in 1:J) {
        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) + 0.5 * 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] <- 2 * test[d]
        psihat[d, 2] <- testit$psihat[1, 2]
        psihat[d, 3] <- testit$test[1, 4]
    }
    list(psihat = psihat, crit.p.value = 2 * crit, con = con)
  }

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