1 |
J |
|
K |
|
x |
|
grp |
|
alpha |
|
p |
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)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.