1 |
x |
|
alpha |
|
nboot |
|
grp |
|
est |
|
con |
|
bhop |
|
SEED |
|
PDIS |
|
J |
|
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 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 | ##---- 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, nboot = 1000, grp = NA, est = tmean,
con = 0, bhop = FALSE, SEED = TRUE, PDIS = FALSE, J = NULL,
p = NULL, ...)
{
con <- as.matrix(con)
if (is.matrix(x) || is.data.frame(x)) {
if (is.null(J) && is.null(p))
stop("Specify J or P")
x = MAT2list(x, p = p, J = J)
}
if (!is.list(x))
stop("Data must be stored in list mode.")
if (!is.na(sum(grp))) {
xx <- list()
for (i in 1:length(grp)) xx[[i]] <- x[[grp[i]]]
x <- xx
}
J <- length(x)
nullvec = rep(0, ncol(x[[1]]))
bplus = nboot + 1
tempn <- 0
mvec <- list
for (j in 1:J) {
x[[j]] <- elimna(x[[j]])
}
Jm <- J - 1
if (sum(con^2) == 0) {
ncon <- (J^2 - J)/2
con <- matrix(0, J, ncon)
id <- 0
for (j in 1:Jm) {
jp <- j + 1
for (k in jp:J) {
id <- id + 1
con[j, id] <- 1
con[k, id] <- 0 - 1
}
}
}
ncon <- ncol(con)
if (nrow(con) != J)
stop("Something is wrong with con; the number of rows does not match the number of groups.")
if (!bhop) {
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)
}
}
if (bhop)
dvec <- (ncon - c(1:ncon) + 1) * alpha/ncon
bvec <- array(NA, c(J, nboot, ncol(x[[1]])))
if (SEED)
set.seed(2)
nvec = lapply(x, nrow)
for (j in 1:J) {
data <- matrix(sample(nvec[[j]], size = nvec[[j]] * nboot,
replace = TRUE), nrow = nboot)
bvec[j, , ] <- apply(data, 1, linconMpb.sub, x[[j]],
est, ...)
}
test <- NA
for (d in 1:ncon) {
tv = matrix(0, nboot, ncol(x[[1]]))
estit = rep(0, ncol(x[[1]]))
for (j in 1:J) {
tv = tv + con[j, d] * bvec[j, , ]
estit = estit + con[j, d] * apply(x[[j]], 2, est,
...)
}
if (!PDIS)
m1 = cov(tv)
tv = rbind(tv, nullvec)
if (!PDIS)
dv = mahalanobis(tv, center = estit, m1)
if (PDIS)
dv = pdis(tv, center = estit)
test[d] = 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
}
output <- matrix(0, ncon, 3)
dimnames(output) <- list(NULL, c("con.num", "p.value", "p.crit"))
temp2 <- order(0 - test)
zvec <- dvec[1:ncon]
sigvec <- (test[temp2] >= zvec)
output[temp2, 3] <- zvec
for (ic in 1:ncol(con)) {
output[ic, 1] <- ic
output[ic, 2] <- test[ic]
}
num.sig <- sum(output[, 2] <= output[, 3])
list(output = output, con = con, num.sig = num.sig)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.