rmrvar:

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

Usage

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

Arguments

x
y
alpha
con
est
plotit
grp
hoch
nboot
xlab
ylab
pr
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 = NA, alpha = 0.05, con = 0, est = pbvar, plotit = FALSE, 
    grp = NA, hoch = TRUE, nboot = NA, xlab = "Group 1", ylab = "Group 2", 
    pr = TRUE, SEED = TRUE, ...) 
{
    if (!is.na(y[1])) 
        x = cbind(x, y)
    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\nnumber of groups.")
        mat <- x
    }
    if (is.matrix(x)) 
        mat <- x
    if (!is.na(sum(grp))) 
        mat <- mat[, grp]
    mat <- elimna(mat)
    x <- mat
    J <- ncol(mat)
    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)) {
        if (d <= 4) 
            nboot <- 1000
        if (d > 4) 
            nboot <- 5000
    }
    n <- nrow(mat)
    crit.vec <- alpha/c(1:d)
    connum <- ncol(con)
    if (SEED) 
        set.seed(2)
    xbars <- apply(mat, 2, est)
    psidat <- NA
    for (ic in 1:connum) psidat[ic] <- sum(con[, ic] * xbars)
    psihat <- matrix(0, connum, nboot)
    bvec <- matrix(NA, ncol = J, nrow = nboot)
    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, ...)
    }
    test <- 1
    bias <- NA
    for (ic in 1:connum) {
        psihat[ic, ] <- apply(bvec, 1, bptdpsi, con[, ic])
        test[ic] <- 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)
        dvecba <- c(0.05, 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)
        dvecba <- c(0.01, 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)
        dvecba <- dvec
        dvec[1] <- alpha/2
    }
    if (hoch) 
        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, ...)
        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", "est.var", "p.value", 
        "crit.p.value", "ci.lower", "ci.upper"))
    tmeans <- apply(mat, 2, est, ...)
    psi <- 1
    for (ic in 1:ncol(con)) {
        output[ic, 2] <- sum(con[, ic] * tmeans)
        output[ic, 1] <- ic
        output[ic, 3] <- test[ic]
        output[temp2, 4] <- zvec
        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]
    }
    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.