aov2depth:

Usage Arguments Examples

Usage

1
aov2depth(x1, x2, est = tmean, nboot = 500, SEED = TRUE, nmin = 12, CR = FALSE, xlab = " DIF 1", ylab = "DIF 2", zlab = "DIF 3", alpha = 0.05, ...)

Arguments

x1
x2
est
nboot
SEED
nmin
CR
xlab
ylab
zlab
alpha
...

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
##---- 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 (x1, x2, est = tmean, nboot = 500, SEED = TRUE, nmin = 12, 
    CR = FALSE, xlab = " DIF 1", ylab = "DIF 2", zlab = "DIF 3", 
    alpha = 0.05, ...) 
{
    if (is.matrix(x1) || is.data.frame(x1)) 
        x1 = listm(x1)
    if (is.matrix(x2) || is.data.frame(x2)) 
        x2 = listm(x2)
    J = length(x1)
    if (J != length(x2)) 
        stop("x1 and x2 should have same number of groups")
    if (SEED) 
        set.seed(2)
    for (j in 1:J) {
        x1[[j]] = na.omit(x1[[j]])
        x2[[j]] = na.omit(x2[[j]])
    }
    n1 = mapply(x1, FUN = length)
    n2 = mapply(x2, FUN = length)
    bplus = nboot + 1
    bvec1 = matrix(NA, nrow = nboot, ncol = J)
    bvec2 = matrix(NA, nrow = nboot, ncol = J)
    for (j in 1:J) {
        data1 = matrix(sample(x1[[j]], size = n1[j] * nboot, 
            replace = TRUE), nrow = nboot)
        data2 = matrix(sample(x2[[j]], size = n2[j] * nboot, 
            replace = TRUE), nrow = nboot)
        bvec1[, j] = apply(data1, 1, est, ...)
        bvec2[, j] = apply(data2, 1, est, ...)
    }
    difb = bvec1 - bvec2
    est1 = mapply(x1, FUN = est, ...)
    est2 = mapply(x2, FUN = est, ...)
    dif = est1 - est2
    m1 = var(difb)
    nullvec = rep(0, J)
    difz = rbind(difb, nullvec)
    dis = mahalanobis(difz, dif, m1)
    sig = sum(dis[bplus] <= dis)/bplus
    if (CR) {
        dis2 <- order(dis[1:nboot])
        dis <- sort(dis)
        critn <- floor((1 - alpha) * nboot)
        if (J == 2) {
            plot(difb[, 1], difb[, 2], xlab = xlab, ylab = ylab)
            points(0, 0, pch = 0)
            xx <- difb[dis2[1:critn], ]
            xord <- order(xx[, 1])
            xx <- xx[xord, ]
            temp <- chull(xx)
            lines(xx[temp, ])
            lines(xx[c(temp[1], temp[length(temp)]), ])
        }
        if (J == 3) {
            scatterplot3d(difb[dis2[1:critn], ], xlab = xlab, 
                ylab = ylab, zlab = zlab, tick.marks = TRUE)
        }
    }
    list(p.value = sig, est1 = est1, est2 = est2, dif = dif, 
        n1 = n1, n2 = n2)
  }

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