anctspb:

In musto101/wilcox_R: What the package does (short line)

Usage

anctspb(x1, y1, x2, y2, pts = NULL, regfun = tshdreg, fr1 = 1, fr2 = 1, alpha = 0.05, plotit = TRUE, xout = FALSE, outfun = outpro, nboot = 500, SEED = TRUE, xlab = "X", ylab = "Y", ...)

Arguments

x1
y1
x2
y2
pts
regfun
fr1
fr2
alpha
plotit
xout
outfun
nboot
SEED
xlab
ylab
...

Examples

##---- 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, pts = NULL, regfun = tshdreg, fr1 = 1, 
    fr2 = 1, alpha = 0.05, plotit = TRUE, xout = FALSE, outfun = outpro, 
    nboot = 500, SEED = TRUE, xlab = "X", ylab = "Y", ...) 
{
    if (SEED) 
        set.seed(2)
    FLAG = pts
    xy = elimna(cbind(x1, y1))
    if (ncol(xy) > 2) 
        stop("Only one covariate is allowed")
    x1 = xy[, 1]
    y1 = xy[, 2]
    xy = elimna(cbind(x2, y2))
    if (ncol(xy) > 2) 
        stop("Only one covariate is allowed")
    x2 = xy[, 1]
    y2 = xy[, 2]
    if (is.null(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, length(pts), 9)
        dimnames(mat) <- list(NULL, c("X", "Est1", "Est2", "DIF", 
            "TEST", "se", "ci.low", "ci.hi", "p.value"))
        pts = x1[isub]
    }
    mat <- matrix(NA, length(pts), 7)
    dimnames(mat) <- list(NULL, c("X", "Est1", "Est2", "DIF", 
        "ci.low", "ci.hi", "p.value"))
    mat[, 1] <- pts
    bvec1 = matrix(NA, nrow = nboot, ncol = length(pts))
    bvec2 = matrix(NA, nrow = nboot, ncol = length(pts))
    x1 = as.matrix(x1)
    x2 = as.matrix(x2)
    p1 = ncol(x1) + 1
    data <- matrix(sample(length(y1), size = length(y1) * nboot, 
        replace = TRUE), nrow = nboot)
    for (ib in 1:nboot) {
        bvec1[ib, ] = regYsub(x1[data[ib, ], ], y1[data[ib, ]], 
            pts, p1 = p1, regfun = regfun, ...)
    }
    data <- matrix(sample(length(y2), size = length(y2) * nboot, 
        replace = TRUE), nrow = nboot)
    for (ib in 1:nboot) {
        bvec2[ib, ] = regYsub(x2[data[ib, ], ], y2[data[ib, ]], 
            pts, p1 = p1, regfun = regfun, ...)
    }
    dif = bvec1 < bvec2
    L = apply(dif, 2, mean)
    E = bvec1 == bvec2
    T = apply(E, 2, mean)
    pvec = L + 0.5 * T
    pop = 1 - pvec
    pb = cbind(pvec, pop)
    pv = 2 * apply(pb, 1, min)
    ilow <- round((alpha/2) * nboot)
    ihi <- nboot - ilow
    ilow <- ilow + 1
    ciL = NA
    ciU = NA
    difb = bvec1 - bvec2
    for (i in 1:length(pts)) {
        bs = sort(difb[, i])
        ciL[i] = bs[ilow]
        ciU[i] = bs[ihi]
    }
    est1 = regYhat(x1, y1, xr = pts, xout = xout, outfun = outfun, 
        regfun = regfun, ...)
    est2 = regYhat(x2, y2, xr = pts, xout = xout, outfun = outfun, 
        regfun = regfun, ...)
    mat[, 2] = est1
    mat[, 3] = est2
    est = est1 - est2
    mat[, 4] = est
    mat[, 5] = ciL
    mat[, 6] = ciU
    mat[, 7] = pv
    if (!is.null(FLAG[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)])
        }
        est1 = regYhat(x1, y1, regfun = regfun, xr = pts, xout = xout, 
            outfun = outfun, ...)
        est2 = regYhat(x2, y2, regfun = regfun, xr = pts, xout = xout, 
            outfun = outfun, ...)
        mat[, 2] = est1
        mat[, 3] = est2
        est = est1 - est2
        mat[, 4] = est
        mat[, 7] = pv
    }
    if (plotit) {
        if (xout) {
            flag <- outfun(x1, ...)$keep
            x1 <- x1[flag]
            y1 <- y1[flag]
            flag <- outfun(x2, ...)$keep
            x2 <- x2[flag]
            y2 <- y2[flag]
        }
        plot(c(x1, x2), c(y1, y2), type = "n", xlab = xlab, ylab = ylab)
        points(x1, y1, pch = "o")
        points(x2, y2, pch = "+")
        abline(regfun(x1, y1)$coef)
        abline(regfun(x2, y2)$coef, lty = 2)
    }
    list(output = mat)
  }

musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.