lsqs3:

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

Usage

lsqs3(x, y, plotit = TRUE, cop = 2, ap.dep = FALSE, v2 = FALSE, pv = FALSE, SEED = TRUE, nboot = 1000, ypch = "o", xpch = "+")

Arguments

x
y
plotit
cop
ap.dep
v2
pv
SEED
nboot
ypch
xpch

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, plotit = TRUE, cop = 2, ap.dep = FALSE, v2 = FALSE, 
    pv = FALSE, SEED = TRUE, nboot = 1000, ypch = "o", xpch = "+") 
{
    if (is.list(x)) 
        x <- matl(x)
    if (is.list(y)) 
        y <- matl(y)
    x <- elimna(x)
    y <- elimna(y)
    x <- as.matrix(x)
    y <- as.matrix(y)
    nx = nrow(x)
    ny = nrow(y)
    if (ncol(x) != ncol(y)) 
        stop("Number of variables not equal")
    disyx <- NA
    disxy <- NA
    if (ncol(x) == 2) {
        if (plotit) {
            plot(rbind(x, y), type = "n", xlab = "VAR 1", ylab = "VAR 2")
            points(x, pch = xpch)
            points(y, pch = ypch)
            if (nrow(x) > 50) {
                if (!ap.dep) {
                  print("If execution time is high, might use ap.dep=F")
                }
                if (!ap.dep) 
                  temp <- depth2(x, plotit = FALSE)
                if (ap.dep) 
                  temp <- fdepth(x, plotit = FALSE, cop = cop)
            }
            if (!ap.dep) 
                temp <- depth2(x, plotit = FALSE)
            if (ap.dep) 
                temp <- fdepth(x, plotit = FALSE, cop = cop)
            flag <- (temp >= median(temp))
            xx <- x[flag, ]
            xord <- order(xx[, 1])
            xx <- xx[xord, ]
            temp <- chull(xx)
            xord <- order(xx[, 1])
            xx <- xx[xord, ]
            temp <- chull(xx)
            lines(xx[temp, ])
            lines(xx[c(temp[1], temp[length(temp)]), ])
            if (ap.dep) 
                temp <- fdepth(y, plotit = FALSE, cop = cop)
            if (!ap.dep) 
                temp <- depth2(y, plotit = FALSE)
            if (!ap.dep) 
                temp <- depth2(y, plotit = FALSE)
            if (!ap.dep) 
                temp <- fdepth(y, plotit = FALSE)
            flag <- (temp >= median(temp))
            xx <- y[flag, ]
            xord <- order(xx[, 1])
            xx <- xx[xord, ]
            temp <- chull(xx)
            flag <- (temp >= median(temp))
            xord <- order(xx[, 1])
            xx <- xx[xord, ]
            temp <- chull(xx)
            lines(xx[temp, ], lty = 2)
            lines(xx[c(temp[1], temp[length(temp)]), ], lty = 2)
        }
        tempyx <- NA
        tempxy <- NA
        if (ap.dep) 
            tempyx <- fdepth(x, y, plotit = FALSE, cop = cop)
        if (!ap.dep) 
            tempyx <- depth2(x, y, plotit = FALSE)
        if (ap.dep) 
            tempxy <- fdepth(y, x, plotit = FALSE, cop = cop)
        tempxy <- depth2(y, x, plotit = FALSE)
    }
    if (ncol(x) == 1) {
        tempyx <- unidepth(as.vector(x), as.vector(y))
        tempxy <- unidepth(as.vector(y), as.vector(x))
    }
    if (ncol(x) > 2) {
        if (!v2) {
            tempxy <- fdepth(y, x, plotit = FALSE, cop = cop)
            tempyx <- fdepth(x, y, plotit = FALSE, cop = cop)
        }
        if (v2) {
            tempxy <- fdepthv2(y, x, plotit = FALSE, cop = cop)
            tempyx <- fdepthv2(x, y, plotit = FALSE, cop = cop)
        }
    }
    qhatxy <- mean(tempxy)
    qhatyx <- mean(tempyx)
    qhat <- max(c(qhatxy, qhatyx))
    n1 <- nrow(x)
    n2 <- nrow(y)
    nv <- (3 * min(c(n1, n2)) + max(c(n1, n2)))/4
    if (ncol(x) == 1) 
        crit <- 0.2536 - 0.4578/sqrt(nv)
    if (ncol(x) == 2) 
        crit <- 0.1569 - 0.3/sqrt(nv)
    if (ncol(x) == 3) 
        crit <- 0.0861 - 0.269/sqrt(nv)
    if (ncol(x) == 4) 
        crit <- 0.054 - 0.1568/sqrt(nv)
    if (ncol(x) == 5) 
        crit <- 0.0367 - 0.0968/sqrt(nv)
    if (ncol(x) == 6) 
        crit <- 0.0262 - 0.0565/sqrt(nv)
    if (ncol(x) == 7) 
        crit <- 0.0174 - 0.0916/sqrt(nv)
    if (ncol(x) > 7) 
        crit <- 0.013
    rej <- "Fail to reject"
    if (qhat <= crit) 
        rej <- "Reject"
    testv = NULL
    pval = NULL
    if (pv) {
        if (SEED) 
            set.seed(2)
        rej = "NULL"
        for (i in 1:nboot) testv[i] = lsqs3.sub(rmul(n1, ncol(x)), 
            rmul(n2, ncol(x)), cop = cop, ap.dep = ap.dep, v2 = v2, 
            )$test
        pval = mean(qhat >= testv)
    }
    list(n1 = nx, n2 = ny, avg.depth.of.x.in.y = qhatxy, avg.depth.of.y.in.x = qhatyx, 
        test = qhat, crit = crit, Decision = rej, p.value = pval)
  }

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