1 |
x1 |
|
x2 |
|
est |
|
nboot |
|
SEED |
|
nmin |
|
CR |
|
xlab |
|
ylab |
|
zlab |
|
alpha |
|
... |
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 | ##---- 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 (x1, x2, est = tmean, nboot = 500, SEED = TRUE, nmin = 12,
CR = FALSE, xlab = " DIF 1", ylab = "DIF 2", zlab = "DIF 3",
alpha = 0.05, ...)
{
if (is.matrix(x1) || is.data.frame(x1))
x1 = listm(x1)
if (is.matrix(x2) || is.data.frame(x2))
x2 = listm(x2)
J = length(x1)
if (J != length(x2))
stop("x1 and x2 should have same number of groups")
if (SEED)
set.seed(2)
for (j in 1:J) {
x1[[j]] = na.omit(x1[[j]])
x2[[j]] = na.omit(x2[[j]])
}
n1 = mapply(x1, FUN = length)
n2 = mapply(x2, FUN = length)
bplus = nboot + 1
bvec1 = matrix(NA, nrow = nboot, ncol = J)
bvec2 = matrix(NA, nrow = nboot, ncol = J)
for (j in 1:J) {
data1 = matrix(sample(x1[[j]], size = n1[j] * nboot,
replace = TRUE), nrow = nboot)
data2 = matrix(sample(x2[[j]], size = n2[j] * nboot,
replace = TRUE), nrow = nboot)
bvec1[, j] = apply(data1, 1, est, ...)
bvec2[, j] = apply(data2, 1, est, ...)
}
difb = bvec1 - bvec2
est1 = mapply(x1, FUN = est, ...)
est2 = mapply(x2, FUN = est, ...)
dif = est1 - est2
m1 = var(difb)
nullvec = rep(0, J)
difz = rbind(difb, nullvec)
dis = mahalanobis(difz, dif, m1)
sig = sum(dis[bplus] <= dis)/bplus
if (CR) {
dis2 <- order(dis[1:nboot])
dis <- sort(dis)
critn <- floor((1 - alpha) * nboot)
if (J == 2) {
plot(difb[, 1], difb[, 2], xlab = xlab, ylab = ylab)
points(0, 0, pch = 0)
xx <- difb[dis2[1:critn], ]
xord <- order(xx[, 1])
xx <- xx[xord, ]
temp <- chull(xx)
lines(xx[temp, ])
lines(xx[c(temp[1], temp[length(temp)]), ])
}
if (J == 3) {
scatterplot3d(difb[dis2[1:critn], ], xlab = xlab,
ylab = ylab, zlab = zlab, tick.marks = TRUE)
}
}
list(p.value = sig, est1 = est1, est2 = est2, dif = dif,
n1 = n1, n2 = n2)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.