1 |
x1 |
|
y1 |
|
x2 |
|
y2 |
|
regfun |
|
xout |
|
outfun |
|
nboot |
|
alpha |
|
SEED |
|
plotit |
|
xlab |
|
ylab |
|
... |
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 | ##---- 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, regfun = lsfit, xout = FALSE, outfun = outpro,
nboot = 200, alpha = 0.05, SEED = TRUE, plotit = FALSE, xlab = "X",
ylab = "Y", ...)
{
if (SEED)
set.seed(2)
X = elimna(cbind(x1, y1, x2, y2))
x1 = as.matrix(x1)
x2 = as.matrix(x2)
p = ncol(x1)
p1 = p + 1
p2 = p + 2
p3 = p1 + p
p4 = p3 + 1
x1 = X[, 1:p]
y1 = X[, p1]
x2 = X[, p2:p3]
y2 = X[, p4]
n = length(y1)
if (xout) {
flag1 = outfun(x1)$out.id
flag2 = outfun(x2)$out.id
flag = unique(c(flag1, flag2))
if (length(flag) > 0)
X = X[-flag, ]
x1 = X[, 1:p]
y1 = X[, p1]
x2 = X[, p2:p3]
y2 = X[, p4]
}
nk = length(y1)
x1 = as.matrix(x1)
x2 = as.matrix(x2)
data <- matrix(sample(length(y1), size = length(y1) * nboot,
replace = TRUE), nrow = nboot)
bvec1 <- apply(data, 1, regboot, x1, y1, regfun, ...)
bvec2 <- apply(data, 1, regboot, x2, y2, regfun, ...)
dif = t(bvec1) - t(bvec2)
est1 = lsfit(x1, y1)$coef
est2 = lsfit(x2, y2)$coef
estdif = est1 - est2
se = apply(dif, 2, sd)
pvec = NA
test = NA
test = estdif/se
df = nk - 1
pvec = 2 * (1 - pt(abs(test), df))
if (plotit) {
reg2plot(x1, y1, x2, y2, xlab = xlab, ylab = ylab)
}
lvec = "Intercept"
ci = matrix(NA, nrow = p1, ncol = 3)
ci[, 1] = c(0:p)
ci[, 2] = estdif + qt(alpha/2, df) * se
ci[, 3] = estdif - qt(alpha/2, df) * se
dimnames(ci) = list(NULL, c("Param", "ci.low", "ci.hi"))
for (j in 2:p1) lvec = c(lvec, paste("slope", j - 1))
pvec = array(pvec, dimnames = lvec)
list(n = n, n.keep = nk, est.dif = estdif, est.1 = est1,
est.2 = est2, test.stat = test, standard.error = se,
p.values = pvec, conf.intervals = ci)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.