ddep:

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

Usage

ddep(x, est = onestep, alpha = 0.05, grp = NA, nboot = 500, plotit = TRUE, SEED = TRUE, pr = TRUE, WT = TRUE, ...)

Arguments

x
est
alpha
grp
nboot
plotit
SEED
pr
WT
...

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 = onestep, alpha = 0.05, grp = NA, nboot = 500, 
    plotit = TRUE, SEED = TRUE, pr = TRUE, WT = TRUE, ...) 
{
    if (pr) 
        print("Warning: Might not be level robust if the number of groups is relatively large and n is small")
    if (pr) 
        print("Currently seems that rmmismcp is preferable")
    if (is.list(x)) {
        nv <- NA
        for (j in 1:length(x)) nv[j] <- length(x[[j]])
        if (var(nv) != 0) {
            stop("The groups are stored in list mode and appear to have different sample sizes")
        }
        temp <- matrix(NA, ncol = length(x), nrow = nv[1])
        for (j in 1:length(x)) temp[, j] <- x[[j]]
        x <- temp
    }
    J <- ncol(x)
    if (!is.na(grp[1])) {
        J <- length(grp)
        for (j in 1:J) temp[, j] <- x[, grp[j]]
        x <- temp
    }
    x <- elimna(x)
    bvec <- matrix(0, ncol = J, nrow = nboot)
    hval <- vector("numeric", J)
    if (SEED) 
        set.seed(2)
    print("Taking bootstrap samples. Please wait.")
    n <- nrow(x)
    totv <- apply(x, 2, est, na.rm = TRUE, ...)
    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, ...)
    if (!WT) {
        gv <- rep(mean(totv), J)
    }
    bplus <- nboot + 1
    center <- totv
    cmat <- var(bvec)
    if (WT) {
        wt = 1/diag(cmat)
        ut = sum(wt)
        gv <- rep(sum(wt * totv)/ut, J)
    }
    m1 <- rbind(bvec, gv)
    discen <- mahalanobis(m1, totv, cmat)
    if (plotit && ncol(x) == 2) {
        plot(bvec, xlab = "Group 1", ylab = "Group 2")
        temp.dis <- order(discen[1:nboot])
        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)
    }
    sig.level <- sum(discen[bplus] <= discen)/bplus
    list(p.value = sig.level, center = totv, grand.mean = gv)
  }

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