1 | regYci2G(x1, y1, x2, y2, regfun = tsreg, pts = NULL, ALL = FALSE, npts = 25, plotit = TRUE, SCAT = TRUE, pch1 = "*", pch2 = "+", nboot = 100, ADJ = TRUE, xout = FALSE, outfun = out, SEED = TRUE, p.crit = 0.015, alpha = 0.05, crit = NULL, null.value = 0, plotPV = FALSE, scale = TRUE, span = 0.75, xlab = "X", xlab1 = "X1", xlab2 = "X2", ylab = "p-values", ylab2 = "Y", theta = 50, phi = 25, MC = FALSE, nreps = 1000, pch = "*", ...)
|
x1 |
|
y1 |
|
x2 |
|
y2 |
|
regfun |
|
pts |
|
ALL |
|
npts |
|
plotit |
|
SCAT |
|
pch1 |
|
pch2 |
|
nboot |
|
ADJ |
|
xout |
|
outfun |
|
SEED |
|
p.crit |
|
alpha |
|
crit |
|
null.value |
|
plotPV |
|
scale |
|
span |
|
xlab |
|
xlab1 |
|
xlab2 |
|
ylab |
|
ylab2 |
|
theta |
|
phi |
|
MC |
|
nreps |
|
pch |
|
... |
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 | ##---- 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 = tsreg, pts = NULL, ALL = FALSE,
npts = 25, plotit = TRUE, SCAT = TRUE, pch1 = "*", pch2 = "+",
nboot = 100, ADJ = TRUE, xout = FALSE, outfun = out, SEED = TRUE,
p.crit = 0.015, alpha = 0.05, crit = NULL, null.value = 0,
plotPV = FALSE, scale = TRUE, span = 0.75, xlab = "X", xlab1 = "X1",
xlab2 = "X2", ylab = "p-values", ylab2 = "Y", theta = 50,
phi = 25, MC = FALSE, nreps = 1000, pch = "*", ...)
{
xy = elimna(cbind(x1, y1))
x1 <- as.matrix(x1)
p = ncol(x1)
if (p > 1)
stop("Current version allows one covariate only")
p1 = p + 1
vals = NA
x1 <- xy[, 1:p]
y1 <- xy[, p1]
x1 <- as.matrix(x1)
xy = elimna(cbind(x2, y2))
x2 <- as.matrix(x2)
p = ncol(x2)
p1 = p + 1
vals = NA
x2 <- xy[, 1:p]
y2 <- xy[, p1]
x2 <- as.matrix(x2)
n1 = length(y1)
n2 = length(y2)
n = min(c(n1, n2))
if (xout) {
m <- cbind(x1, y1)
flag <- outfun(x1, plotit = FALSE, ...)$keep
m <- m[flag, ]
n1 = nrow(m)
x1 <- m[, 1:p]
y1 <- m[, p1]
x1 = as.matrix(x1)
m <- cbind(x2, y2)
flag <- outfun(x2, plotit = FALSE, ...)$keep
m <- m[flag, ]
n2 = nrow(m)
n = min(c(n1, n2))
x2 <- m[, 1:p]
y2 <- m[, p1]
x2 = as.matrix(x2)
}
if (is.null(pts)) {
xall = unique(c(x1, x2))
if (ALL)
pts = xall
if (!ALL)
pts = seq(min(xall), max(xall), length.out = npts)
}
if (ADJ) {
if (n < 10)
stop("Should have a sample size of at least 10")
if (alpha == 0.05) {
alpha = p.crit
crit = qnorm(1 - alpha/2)
}
if (!ADJ)
p.crit = alpha
if (n < 20 & max(c(n1, n2)) < 100)
crit = NULL
if (p > 1)
crit = NULL
}
if (is.null(crit) & !ADJ)
crit = qnorm(1 - alpha/2)
if (is.null(crit) & ADJ) {
if (SEED)
set.seed(2)
padj = regYciCV2G(n1, n2, nboot = nreps, regfun = regfun,
MC = MC, SEED = FALSE, ALL = ALL, null.value = null.value,
pts = pts, alpha = alpha, ...)$crit.est
crit = qnorm(1 - padj/2)
p.crit = padj
}
sqsd1 = regYvar(x1, y1, regfun = regfun, pts = pts, nboot = nboot,
SEED = SEED)
sqsd2 = regYvar(x2, y2, regfun = regfun, pts = pts, nboot = nboot,
SEED = SEED)
sd = sqrt(sqsd1 + sqsd2)
est1 = regYhat(x1, y1, regfun = regfun, xr = pts, ...)
est2 = regYhat(x2, y2, regfun = regfun, xr = pts, ...)
pv = 2 * (1 - pnorm(abs(est1 - est2 - null.value)/sd))
est = cbind(pts, est1 - est2, est1 - est2 - crit * sd, est1 -
est2 + crit * sd, pv)
dimnames(est) = list(NULL, c("X", "Est.Dif", "Lower.ci",
"Upper.ci", "p.value"))
if (plotit) {
plotPV = FALSE
plot(c(x1, x2), c(y1, y2), type = "n", xlab = xlab, ylab = ylab2)
reg1 = regfun(x1, y1, ...)$coef
reg2 = regfun(x2, y2, ...)$coef
if (SCAT) {
points(x1, y1, pch = pch1)
points(x2, y2, pch = pch2)
}
abline(reg1)
abline(reg2, lty = 2)
}
if (plotPV) {
if (ncol(x1) > 2)
stop("Can plot only with one or two independent variables")
if (ncol(x1) == 1)
plot(pts, pv, xlab = xlab, ylab = ylab, pch = pch)
if (ncol(x2) == 2)
lplot(pts, pv, xlab = xlab1, ylab = xlab2, zlab = ylab,
span = span, ticktype = "detail", scale = scale,
theta = theta, phi = phi)
}
list(output = est, p.crit = p.crit, crit.value = crit, num.sig = sum(est[,
5] <= p.crit))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.