med2way:

Usage Arguments Examples

Usage

1
med2way(J, K, x, grp = c(1:p), alpha = 0.05, p = J * K)

Arguments

J
K
x
grp
alpha
p

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
67
68
69
70
71
72
73
74
75
##---- 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 (J, K, x, grp = c(1:p), alpha = 0.05, p = J * K) 
{
    print("Suggestion: Use the function med2way or m2way instead, especially with tied values")
    if (is.matrix(x)) 
        x <- listm(x)
    if (!is.list(x)) 
        stop("Data are not stored in a matrix or in list mode")
    if (p != length(x)) {
        print("Warning: The number of groups in your data is not equal to JK")
    }
    xbar <- 0
    h <- 0
    d <- 0
    R <- 0
    W <- 0
    d <- 0
    r <- 0
    w <- 0
    nuhat <- 0
    omegahat <- 0
    DROW <- 0
    DCOL <- 0
    xtil <- matrix(0, J, K)
    aval <- matrix(0, J, K)
    for (j in 1:p) {
        xbar[j] <- median(x[[grp[j]]])
        h[j] <- length(x[[grp[j]]])
        d[j] <- msmedse(x[[grp[j]]])^2
    }
    d <- matrix(d, J, K, byrow = T)
    xbar <- matrix(xbar, J, K, byrow = T)
    h <- matrix(h, J, K, byrow = T)
    for (j in 1:J) {
        R[j] <- sum(xbar[j, ])
        nuhat[j] <- (sum(d[j, ]))^2/sum(d[j, ]^2/(h[j, ] - 1))
        r[j] <- 1/sum(d[j, ])
        DROW[j] <- sum(1/d[j, ])
    }
    for (k in 1:K) {
        W[k] <- sum(xbar[, k])
        omegahat[k] <- (sum(d[, k]))^2/sum(d[, k]^2/(h[, k] - 
            1))
        w[k] <- 1/sum(d[, k])
        DCOL[k] <- sum(1/d[, k])
    }
    D <- 1/d
    for (j in 1:J) {
        for (k in 1:K) {
            xtil[j, k] <- sum(D[, k] * xbar[, k]/DCOL[k]) + sum(D[j, 
                ] * xbar[j, ]/DROW[j]) - sum(D * xbar/sum(D))
            aval[j, k] <- (1 - D[j, k] * (1/sum(D[j, ]) + 1/sum(D[, 
                k]) - 1/sum(D)))^2/(h[j, k] - 3)
        }
    }
    Rhat <- sum(r * R)/sum(r)
    What <- sum(w * W)/sum(w)
    Ba <- sum((1 - r/sum(r))^2/nuhat)
    Bb <- sum((1 - w/sum(w))^2/omegahat)
    Va <- sum(r * (R - Rhat)^2)/((J - 1) * (1 + 2 * (J - 2) * 
        Ba/(J^2 - 1)))
    Vb <- sum(w * (W - What)^2)/((K - 1) * (1 + 2 * (K - 2) * 
        Bb/(K^2 - 1)))
    sig.A <- 1 - pf(Va, J - 1, 9999999)
    sig.B <- 1 - pf(Vb, K - 1, 9999999)
    Vab <- sum(D * (xbar - xtil)^2)
    dfinter <- (J - 1) * (K - 1)
    sig.AB <- 1 - pchisq(Vab, dfinter)
    list(test.A = Va, p.val.A = sig.A, test.B = Vb, p.val.B = sig.B, 
        test.AB = Vab, p.val.AB = sig.AB)
  }

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