medr:

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

Usage

medr(x, est = median, alpha = 0.05, nboot = 500, grp = NA, op = 1, MM = FALSE, cop = 3, pr = TRUE, SEED = TRUE, ...)

Arguments

x
est
alpha
nboot
grp
op
MM
cop
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, est = median, alpha = 0.05, nboot = 500, grp = NA, 
    op = 1, MM = FALSE, cop = 3, pr = TRUE, SEED = TRUE, ...) 
{
    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(grp)) {
        xx <- list()
        for (i in 1:length(grp)) xx[[i]] <- x[[grp[1]]]
        x <- xx
    }
    J <- length(x)
    mvec <- NA
    for (j in 1:J) {
        temp <- x[[j]]
        temp <- temp[!is.na(temp)]
        x[[j]] <- temp
        mvec[j] <- est(temp, ...)
    }
    Jm <- J - 1
    d <- (J^2 - J)/2
    data <- list()
    bvec <- matrix(NA, ncol = d, nrow = nboot)
    if (SEED) 
        set.seed(2)
    if (pr) 
        print("Taking bootstrap samples. Please wait.")
    for (it in 1:nboot) {
        for (j in 1:J) data[[j]] <- sample(x[[j]], size = length(x[[j]]), 
            replace = TRUE)
        dval <- 0
        for (j in 1:J) {
            for (k in 1:J) {
                if (j < k) {
                  dval <- dval + 1
                  bvec[it, dval] <- loc2dif(data[[j]], data[[k]], 
                    est = est, ...)
                }
            }
        }
    }
    output <- matrix(NA, nrow = d, ncol = 3)
    dimnames(output) <- list(NULL, c("Group", "Group", "psihat"))
    tvec <- NA
    dval <- 0
    for (j in 1:J) {
        for (k in 1:J) {
            if (j < k) {
                dval <- dval + 1
                output[dval, 1] <- j
                output[dval, 2] <- k
                tvec[dval] <- loc2dif(x[[j]], x[[k]], est = est, 
                  ...)
                output[dval, 3] <- tvec[dval]
            }
        }
    }
    tempcen <- apply(bvec, 1, mean)
    vecz <- rep(0, d)
    smat <- var(bvec - tempcen + tvec)
    temp <- bvec - tempcen + tvec
    bcon <- rbind(bvec, vecz)
    if (op == 1) 
        dv <- mahalanobis(bcon, tvec, smat)
    if (op == 2) {
        smat <- cov.mcd(temp)$cov
        dv <- mahalanobis(bcon, tvec, smat)
    }
    if (op == 3) {
        print("Computing p-value. Might take a while with op=3")
        dv <- pdis(bcon, MM = MM, cop = cop, center = tvec)
    }
    if (op == 4) 
        dv <- pdis.for(bcon, MM = MM, cop = cop, pr = FALSE, 
            center = tvec)
    bplus <- nboot + 1
    sig.level <- 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
    if (op == 4) 
        print(sig.level)
    list(sig.level = sig.level, output = output)
  }

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