1 |
x |
|
y |
|
tr |
|
alpha |
|
plotit |
|
plotfun |
|
op |
|
VL |
|
cor.op |
|
xlab |
|
ylab |
|
PB |
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 | ##---- 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, y, tr = 0.2, alpha = 0.05, plotit = FALSE, plotfun = splot,
op = TRUE, VL = TRUE, cor.op = FALSE, xlab = "Groups", ylab = "",
PB = FALSE)
{
if (tr == 0.5)
stop("Use medpb to compare medians.")
if (tr > 0.5)
stop("Can't have tr>.5")
library(MASS)
x <- x[!is.na(x)]
y <- y[!is.na(y)]
h1 <- length(x) - 2 * floor(tr * length(x))
h2 <- length(y) - 2 * floor(tr * length(y))
q1 <- (length(x) - 1) * winvar(x, tr)/(h1 * (h1 - 1))
q2 <- (length(y) - 1) * winvar(y, tr)/(h2 * (h2 - 1))
df <- (q1 + q2)^2/((q1^2/(h1 - 1)) + (q2^2/(h2 - 1)))
crit <- qt(1 - alpha/2, df)
m1 = mean(x, tr)
m2 = mean(y, tr)
mbar = (m1 + m2)/2
dif = m1 - m2
low <- dif - crit * sqrt(q1 + q2)
up <- dif + crit * sqrt(q1 + q2)
test <- abs(dif/sqrt(q1 + q2))
yuen <- 2 * (1 - pt(test, df))
xx = c(rep(1, length(x)), rep(2, length(y)))
pts = c(x, y)
top = var(c(m1, m2))
if (!PB) {
if (tr == 0)
cterm = 1
if (tr > 0)
cterm = area(dnormvar, qnorm(tr), qnorm(1 - tr)) +
2 * (qnorm(tr)^2) * tr
bot = winvar(pts, tr = tr)/cterm
}
if (PB)
bot = pbvar(pts)/1.06
e.pow = top/bot
if (e.pow > 1) {
x0 = c(rep(1, length(x)), rep(2, length(y)))
y0 = c(x, y)
e.pow = wincor(x0, y0, tr = tr)$cor^2
}
if (plotit) {
plot(xx, pts, xlab = xlab, ylab = ylab)
if (op)
points(c(1, 2), c(m1, m2))
if (VL)
lines(c(1, 2), c(m1, m2))
}
list(ci = c(low, up), p.value = yuen, dif = dif, se = sqrt(q1 +
q2), teststat = test, crit = crit, df = df, Var.Explained = e.pow,
Effect.Size = sqrt(e.pow))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.