1 |
x |
|
x2 |
|
cil |
|
crit |
|
con |
|
tr |
|
alpha |
|
iter |
|
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 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 | ##---- 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, x2 = NA, cil = NA, crit = NA, con = 0, tr = 0.2,
alpha = 0.05, iter = 10000, SEED = TRUE)
{
x3 <- x2
if (is.matrix(x))
x <- listm(x)
if (!is.list(x))
stop("Data must be stored in list mode or in matrix mode.")
J <- length(x)
tempn <- 0
svec <- NA
for (j in 1:J) {
temp <- x[[j]]
temp <- temp[!is.na(temp)]
tempn[j] <- length(temp)
x[[j]] <- temp
svec[j] <- winvar(temp, tr = tr)/(1 - 2 * tr)^2
}
tempt <- floor((1 - 2 * tr) * tempn)
A <- sum(1/(tempt - 1))
df <- J/A
print(paste("If using the tables of Studentized range distribution,"))
print(paste("the degrees of freedom are:", df))
if (!is.list(x2) && !is.matrix(x2)) {
x2 <- list()
for (j in 1:J) x2[[j]] <- NA
}
if (is.na(cil))
stop("To proceed, you must specify the maximum length of the confidence intervals.")
if (is.na(crit)) {
print("Approximating critical value")
crit <- trange(tempn - 1, alpha = alpha, iter = iter,
SEED = SEED)
print(paste("The critical value is ", crit))
}
if (con[1] == 0) {
Jm <- J - 1
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)
avec <- NA
for (i in 1:ncon) {
temp <- con[, i]
avec[i] <- sum(temp[temp > 0])
}
dvec <- (cil/(2 * crit * avec))^2
d <- max(dvec)
n.vec <- NA
for (j in 1:J) {
n.vec[j] <- max(tempn[j], floor(svec[j]/d) + 1)
print(paste("Need an additional ", n.vec[j] - tempn[j],
" observations for group", j))
}
if (is.matrix(x2))
x2 <- listm(x2)
temp2 <- n.vec - tempn
if (!is.list(x3) && !is.matrix(x3) && sum(temp2) > 0)
stop("No second stage data supplied; this function is terminating")
if (length(x) != length(x2))
warning("Number of groups in first stage data does not match the number in the second stage.")
ci.mat <- NA
if (!is.na(x2[1]) || sum(temp2) == 0) {
xtil <- NA
nvec2 <- NA
for (j in 1:J) {
nvec2[j] <- 0
temp <- x2[[j]]
if (!is.na(temp[1]))
nvec2[j] <- length(x2[[j]])
if (nvec2[j] < n.vec[j] - tempn[j])
warning(paste("The required number of observations for group",
j, " in the second stage is ", n.vec[j] - tempn[j],
" but only ", nvec2[j], " are available"))
xtil[j] <- mean(c(x[[j]], x2[[j]]), tr = tr, na.rm = TRUE)
}
ci.mat <- matrix(0, ncol = 3, nrow = ncon)
dimnames(ci.mat) <- list(NULL, c("con.num", "ci.low",
"ci.high"))
for (ic in 1:ncon) {
ci.mat[ic, 1] <- ic
bvec <- con[, ic] * sqrt(svec/(nvec2 + tempn))
A <- sum(bvec[bvec > 0])
C <- 0 - sum(bvec[bvec < 0])
D <- max(A, C)
ci.mat[ic, 2] <- sum(con[, ic] * xtil) - crit * D
ci.mat[ic, 3] <- sum(con[, ic] * xtil) + crit * D
}
}
list(ci.mat = ci.mat, con = con)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.