1 |
x |
|
alpha |
|
bop |
|
nboot |
|
pr |
|
q |
|
SEED |
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 76 77 78 79 80 81 82 83 84 85 86 87 88 89 | ##---- 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, alpha = 0.05, bop = FALSE, nboot = 100, pr = TRUE,
q = 0.5, SEED = TRUE)
{
if (is.data.frame(x))
x = as.matrix(x)
if (!is.matrix(x))
x <- matl(x)
if (!is.matrix(x))
stop("Data must be stored in a matrix or in list mode.")
J <- ncol(x)
xbar <- vector("numeric", J)
x <- elimna(x)
df <- nrow(x) - 1
nval <- nrow(x)
for (j in 1:J) {
if (!bop)
xbar[j] <- qest(x[, j], q = q)
if (bop)
xbar[j] <- median(x[, j])
}
CC <- (J^2 - J)/2
ncon <- CC
if (alpha == 0.05) {
dvec <- c(0.05, 0.025, 0.0169, 0.0127, 0.0102, 0.00851,
0.0073, 0.00639, 0.00568, 0.00511)
if (ncon > 10) {
avec <- 0.05/c(11:ncon)
dvec <- c(dvec, avec)
}
}
if (alpha == 0.01) {
dvec <- c(0.01, 0.005, 0.00334, 0.00251, 0.00201, 0.00167,
0.00143, 0.00126, 0.00112, 0.00101)
if (ncon > 10) {
avec <- 0.01/c(11:ncon)
dvec <- c(dvec, avec)
}
}
if (alpha != 0.05 && alpha != 0.01)
dvec <- alpha/c(1:ncon)
psihat <- matrix(0, CC, 5)
dimnames(psihat) <- list(NULL, c("Group", "Group", "psihat",
"ci.lower", "ci.upper"))
test <- matrix(NA, CC, 6)
dimnames(test) <- list(NULL, c("Group", "Group", "test",
"p-value", "p.crit", "se"))
if (bop)
se.val <- bootdse(x, nboot = nboot, pr = pr)
temp1 <- 0
jcom <- 0
for (j in 1:J) {
for (k in 1:J) {
if (j < k) {
jcom <- jcom + 1
if (!bop)
temp <- qdtest(x[, j], x[, k], q = q, bop = bop)
if (bop)
temp <- qdtest(x[, j], x[, k], se.val = se.val[jcom])
sejk <- temp$se
test[jcom, 6] <- sejk
test[jcom, 3] <- temp$test.stat
test[jcom, 4] <- temp$p.value
if (length(x[, j]) < 20)
test[jcom, 4] <- mrm1way(x[, c(j, k)], q = q,
SEED = SEED)$p.value
psihat[jcom, 1] <- j
psihat[jcom, 2] <- k
test[jcom, 1] <- j
test[jcom, 2] <- k
psihat[jcom, 3] <- (xbar[j] - xbar[k])
}
}
}
temp1 <- test[, 4]
temp2 <- order(0 - temp1)
zvec <- dvec[1:ncon]
test[temp2, 5] <- zvec
psihat[, 4] <- psihat[, 3] - qt(1 - test[, 5]/2, df) * test[,
6]
psihat[, 5] <- psihat[, 3] + qt(1 - test[, 5]/2, df) * test[,
6]
num.sig <- sum(test[, 4] <= test[, 5])
list(test = test, psihat = psihat, num.sig = num.sig)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.