rmmismcp:

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

Usage

rmmismcp(x, y = NA, alpha = 0.05, con = 0, est = tmean, plotit = TRUE, grp = NA, nboot = 500, SEED = TRUE, xlab = "Group 1", ylab = "Group 2", pr = FALSE, ...)

Arguments

x
y
alpha
con
est
plotit
grp
nboot
SEED
xlab
ylab
pr
...

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 = NA, alpha = 0.05, con = 0, est = tmean, plotit = TRUE, 
    grp = NA, nboot = 500, SEED = TRUE, xlab = "Group 1", ylab = "Group 2", 
    pr = FALSE, ...) 
{
    if (!is.na(y[1])) 
        x <- cbind(x, y)
    if (is.list(x)) 
        x = matl(x)
    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
    }
    J <- ncol(x)
    Jm <- J - 1
    flag.con = F
    if (sum(con^2) == 0) {
        flag.con = T
        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)
    n <- nrow(x)
    crit.vec <- alpha/c(1:d)
    connum <- ncol(con)
    if (SEED) 
        set.seed(2)
    xbars <- apply(x, 2, est, na.rm = TRUE)
    psidat <- NA
    bveccen <- matrix(NA, ncol = J, nrow = nboot)
    for (ic in 1:connum) psidat[ic] <- sum(con[, ic] * xbars)
    psihat <- matrix(0, connum, nboot)
    psihatcen <- matrix(0, connum, nboot)
    bvec <- matrix(NA, ncol = J, nrow = nboot)
    if (pr) 
        print("Taking bootstrap samples. Please wait.")
    data <- matrix(sample(n, size = n * nboot, replace = TRUE), 
        nrow = nboot)
    for (ib in 1:nboot) {
        bvec[ib, ] <- apply(x[data[ib, ], ], 2, est, na.rm = TRUE, 
            ...)
    }
    test <- 1
    for (ic in 1:connum) {
        for (ib in 1:nboot) {
            psihat[ic, ib] = sum(con[, ic] * bvec[ib, ])
        }
        matcon = c(0, psihat[ic, ])
        dis = mean((psihat[ic, ] < 0)) + 0.5 * mean((psihat[ic, 
            ] == 0))
        test[ic] <- 2 * min(c(dis, 1 - dis))
    }
    ncon <- ncol(con)
    dvec <- alpha/c(1:ncon)
    if (plotit && ncol(bvec) == 2) {
        z <- c(0, 0)
        one <- c(1, 1)
        plot(rbind(bvec, z, one), xlab = xlab, ylab = ylab, type = "n")
        points(bvec)
        totv <- apply(x, 2, est, na.rm = TRUE, ...)
        cmat <- var(bvec)
        dis <- mahalanobis(bvec, totv, cmat)
        temp.dis <- order(dis)
        ic <- round((1 - alpha) * nboot)
        xx <- bvec[temp.dis[1:ic], ]
        xord <- order(xx[, 1])
        xx <- xx[xord, ]
        temp <- chull(xx)
        lines(xx[temp, ])
        lines(xx[c(temp[1], temp[length(temp)]), ])
        abline(0, 1)
    }
    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", 
        "crit.sig", "ci.lower", "ci.upper"))
    tmeans <- apply(x, 2, est, na.rm = TRUE, ...)
    psi <- 1
    output[temp2, 4] <- zvec
    for (ic in 1:ncol(con)) {
        output[ic, 2] <- sum(con[, ic] * tmeans)
        output[ic, 1] <- ic
        output[ic, 3] <- test[ic]
        temp <- sort(psihat[ic, ])
        icl <- round(output[ic, 4] * nboot/2) + 1
        icu <- nboot - (icl - 1)
        output[ic, 5] <- temp[icl]
        output[ic, 6] <- temp[icu]
    }
    if (!flag.con) {
    }
    if (flag.con) {
        CC = (J^2 - J)/2
        test <- matrix(NA, CC, 7)
        dimnames(test) <- list(NULL, c("Group", "Group", "psi.hat", 
            "p.value", "p.crit", "ci.low", "ci.upper"))
        jcom <- 0
        for (j in 1:J) {
            for (k in 1:J) {
                if (j < k) {
                  jcom <- jcom + 1
                  test[jcom, 1] = j
                  test[jcom, 2] = k
                  test[jcom, 3:5] = output[jcom, 2:4]
                  test[jcom, 6:7] = output[jcom, 5:6]
                  con = NULL
                }
            }
        }
    }
    if (!flag.con) 
        test = output
    if (flag.con) 
        num.sig <- sum(test[, 4] <= test[, 5])
    if (!flag.con) 
        num.sig <- sum(test[, 3] <= test[, 4])
    list(output = test, con = con, num.sig = num.sig)
  }

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