1 |
x1 |
|
y1 |
|
x2 |
|
y2 |
|
fr1 |
|
fr2 |
|
alpha |
|
plotit |
|
pts |
|
sm |
|
pr |
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 112 113 114 115 116 117 118 119 120 121 | ##---- 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, y1, x2, y2, fr1 = 1, fr2 = 1, alpha = 0.05, plotit = TRUE,
pts = NA, sm = FALSE, pr = T)
{
if (pr) {
print("NOTE: Confidence intervals are adjusted to control the probability")
print("of at least one Type I error.")
print("But p-values are not")
}
if (is.na(pts[1])) {
npt <- 5
isub <- c(1:5)
test <- c(1:5)
xorder <- order(x1)
y1 <- y1[xorder]
x1 <- x1[xorder]
xorder <- order(x2)
y2 <- y2[xorder]
x2 <- x2[xorder]
n1 <- 1
n2 <- 1
vecn <- 1
for (i in 1:length(x1)) n1[i] <- length(y1[near(x1, x1[i],
fr1)])
for (i in 1:length(x1)) n2[i] <- length(y2[near(x2, x1[i],
fr2)])
for (i in 1:length(x1)) vecn[i] <- min(n1[i], n2[i])
sub <- c(1:length(x1))
isub[1] <- min(sub[vecn >= 12])
isub[5] <- max(sub[vecn >= 12])
isub[3] <- floor((isub[1] + isub[5])/2)
isub[2] <- floor((isub[1] + isub[3])/2)
isub[4] <- floor((isub[3] + isub[5])/2)
mat <- matrix(NA, 5, 9)
dimnames(mat) <- list(NULL, c("X", "n1", "n2", "DIF",
"TEST", "se", "ci.low", "ci.hi", "p.value"))
critv <- NA
if (alpha == 0.05)
critv <- smmcrit(500, 5)
if (alpha == 0.01)
critv <- smmcrit01(500, 5)
if (is.na(critv))
critv <- smmval(rep(999, 5), alpha = alpha)
for (i in 1:5) {
g1 <- y1[near(x1, x1[isub[i]], fr1)]
g2 <- y2[near(x2, x1[isub[i]], fr2)]
g1 <- g1[!is.na(g1)]
g2 <- g2[!is.na(g2)]
test <- msmed(g1, g2)
mat[i, 1] <- x1[isub[i]]
mat[i, 2] <- length(g1)
mat[i, 3] <- length(g2)
mat[i, 4] <- median(g1) - median(g2)
mat[i, 5] <- test$test[3]
mat[i, 6] <- test$test[5]
cilow <- mat[i, 4] - critv * mat[i, 6]
cihi <- mat[i, 4] + critv * mat[i, 6]
mat[i, 7] <- cilow
mat[i, 8] <- cihi
mat[i, 9] <- test$test[6]
}
}
if (!is.na(pts[1])) {
if (length(pts) >= 29)
stop("At most 28 points can be compared")
n1 <- 1
n2 <- 1
vecn <- 1
for (i in 1:length(pts)) {
n1[i] <- length(y1[near(x1, pts[i], fr1)])
n2[i] <- length(y2[near(x2, pts[i], fr2)])
}
mat <- matrix(NA, length(pts), 9)
dimnames(mat) <- list(NULL, c("X", "n1", "n2", "DIF",
"TEST", "se", "ci.low", "ci.hi", "p.value"))
critv <- NA
if (length(pts) >= 2) {
if (alpha == 0.05)
critv <- smmcrit(500, length(pts))
if (alpha == 0.01)
critv <- smmcrit01(500, length(pts))
if (is.na(critv))
critv <- smmval(rep(999, length(pts)), alpha = alpha)
}
if (length(pts) == 1)
critv <- qnorm(1 - alpha/2)
for (i in 1:length(pts)) {
g1 <- y1[near(x1, pts[i], fr1)]
g2 <- y2[near(x2, pts[i], fr2)]
g1 <- g1[!is.na(g1)]
g2 <- g2[!is.na(g2)]
test <- msmed(g1, g2)
mat[i, 1] <- pts[i]
mat[i, 2] <- length(g1)
mat[i, 3] <- length(g2)
if (length(g1) <= 5)
print(paste("Warning, there are", length(g1),
" points corresponding to the design point X=",
pts[i]))
if (length(g2) <= 5)
print(paste("Warning, there are", length(g2),
" points corresponding to the design point X=",
pts[i]))
mat[i, 4] <- median(g1) - median(g2)
mat[i, 5] <- test$test[3]
mat[i, 6] <- test$test[5]
cilow <- mat[i, 4] - critv * mat[i, 6]
cihi <- mat[i, 4] + critv * mat[i, 6]
mat[i, 7] <- cilow
mat[i, 8] <- cihi
mat[i, 9] <- test$test[6]
}
}
if (plotit)
runmean2g(x1, y1, x2, y2, fr = fr1, est = median, sm = sm)
list(output = mat, crit = critv)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.