rmmcppb:

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

Usage

rmmcppb(x, y = NULL, alpha = 0.05, con = 0, est = onestep, plotit = TRUE, dif = TRUE, grp = NA, nboot = NA, BA = FALSE, hoch = FALSE, xlab = "Group 1", ylab = "Group 2", pr = TRUE, SEED = TRUE, SR = FALSE, ...)

Arguments

x
y
alpha
con
est
plotit
dif
grp
nboot
BA
hoch
xlab
ylab
pr
SEED
SR
...

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, dif = TRUE, grp = NA, nboot = NA, BA = FALSE, 
    hoch = FALSE, xlab = "Group 1", ylab = "Group 2", pr = TRUE, 
    SEED = TRUE, SR = FALSE, ...) 
{
    if (hoch) 
        SR = FALSE
    if (SR) {
        okay = FALSE
        if (identical(est, onestep)) 
            okay = TRUE
        if (identical(est, mom)) 
            okay = TRUE
        SR = okay
    }
    if (dif) {
        if (pr) 
            print("dif=T, so analysis is done on difference scores")
        temp <- rmmcppbd(x, y = y, alpha = 0.05, con = con, est, 
            plotit = plotit, grp = grp, nboot = nboot, hoch = TRUE, 
            ...)
        output <- temp$output
        con <- temp$con
    }
    if (!dif) {
        if (pr) {
            print("dif=F, so analysis is done on marginal distributions")
            if (!BA) 
                print("With M-estimator or MOM, suggest using BA=T and hoch=T")
        }
        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]
        mat <- elimna(mat)
        x <- mat
        J <- ncol(mat)
        xcen <- x
        for (j in 1:J) xcen[, j] <- x[, j] - est(x[, j], ...)
        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)
        psihatcen <- matrix(0, connum, nboot)
        bvec <- matrix(NA, ncol = J, nrow = nboot)
        bveccen <- 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, ...)
            bveccen[ib, ] <- apply(xcen[data[ib, ], ], 2, est, 
                ...)
        }
        test <- 1
        bias <- NA
        for (ic in 1:connum) {
            psihat[ic, ] <- apply(bvec, 1, bptdpsi, con[, ic])
            psihatcen[ic, ] <- apply(bveccen, 1, bptdpsi, con[, 
                ic])
            bias[ic] <- sum((psihatcen[ic, ] > 0))/nboot - 0.5
            ptemp <- (sum(psihat[ic, ] > 0) + 0.5 * sum(psihat[ic, 
                ] == 0))/nboot
            if (BA) 
                test[ic] <- ptemp - 0.1 * bias[ic]
            if (!BA) 
                test[ic] <- ptemp
            test[ic] <- min(test[ic], 1 - test[ic])
            test[ic] <- max(test[ic], 0)
        }
        test <- 2 * test
        ncon <- ncol(con)
        dvec <- alpha/c(1:ncon)
        if (SR) {
            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[2] <- alpha
            }
        }
        if (hoch) 
            dvec <- alpha/c(1:ncon)
        dvecba <- dvec
        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]
        if (BA) 
            zvec <- dvecba[1:ncon]
        sigvec <- (test[temp2] >= zvec)
        output <- matrix(0, connum, 6)
        dimnames(output) <- list(NULL, c("con.num", "psihat", 
            "p.value", "p.sig", "ci.lower", "ci.upper"))
        tmeans <- apply(mat, 2, est, ...)
        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(alpha * 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.