opregpb:

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

Usage

opregpb(x, y, nboot = 1000, alpha = 0.05, om = TRUE, ADJ = TRUE, nullvec = rep(0, ncol(x) + 1), plotit = TRUE, opdis = 2, gval = sqrt(qchisq(0.95, ncol(x) + 1)))

Arguments

x
y
nboot
alpha
om
ADJ
nullvec
plotit
opdis
gval

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 (x, y, nboot = 1000, alpha = 0.05, om = TRUE, ADJ = TRUE, 
    nullvec = rep(0, ncol(x) + 1), plotit = TRUE, opdis = 2, 
    gval = sqrt(qchisq(0.95, ncol(x) + 1))) 
{
    x <- as.matrix(x)
    m <- cbind(x, y)
    p1 <- ncol(x) + 1
    m <- elimna(m)
    x <- m[, 1:ncol(x)]
    x <- as.matrix(x)
    y <- m[, p1]
    if (nrow(x) != length(y)) 
        stop("Sample size of x differs from sample size of y")
    if (!is.matrix(x)) 
        stop("Data should be stored in a matrix")
    print("Taking bootstrap samples. Please wait.")
    data <- matrix(sample(length(y), size = length(y) * nboot, 
        replace = TRUE), nrow = nboot)
    bvec <- apply(data, 1, regboot, x, y, regfun = opreg)
    bvec <- t(bvec)
    dvec <- alpha/(c(1:ncol(x)))
    test <- NA
    icl0 <- round(alpha * nboot/2)
    icl <- round(alpha * nboot/(2 * ncol(x)))
    icu0 <- nboot - icl0
    icu <- nboot - icl
    output <- matrix(0, p1, 6)
    vlabs = "Intercept"
    for (j in 2:p1) vlabs[j] = paste("Slope", j - 1)
    dimnames(output) <- list(vlabs, c("Param.", "p.value", "p.crit", 
        "ci.lower", "ci.upper", "s.e."))
    pval <- NA
    for (i in 1:p1) {
        output[i, 1] <- i - 1
        se.val <- var(bvec[, i])
        temp <- sort(bvec[, i])
        output[i, 6] <- sqrt(se.val)
        if (i == 1) {
            output[i, 4] <- temp[icl0 + 1]
            output[i, 5] <- temp[icu0]
        }
        if (i > 1) {
            output[i, 4] <- temp[icl + 1]
            output[i, 5] <- temp[icu]
        }
        pval[i] <- sum((temp > nullvec[i]))/length(temp)
        if (pval[i] > 0.5) 
            pval[i] <- 1 - pval[i]
    }
    fac <- 2
    if (ADJ) {
        nval <- length(y)
        if (nval < 20) 
            nval <- 20
        if (nval > 60) 
            nval <- 60
        fac <- 2 - (60 - nval)/40
    }
    pval[1] <- 2 * pval[1]
    pval[2:p1] <- fac * pval[2:p1]
    output[, 2] <- pval
    temp2 <- order(0 - pval[2:p1])
    zvec <- dvec[1:ncol(x)]
    sigvec <- (test[temp2] >= zvec)
    output[temp2 + 1, 3] <- zvec
    output[1, 3] <- NA
    output[, 2] <- pval
    om.pval <- NA
    temp <- opreg(x, y)$coef
    if (om && ncol(x) > 1) {
        temp2 <- rbind(bvec[, 2:p1], nullvec[2:p1])
        if (opdis == 1) 
            dis <- pdis(temp2, pr = FALSE, center = temp[2:p1])
        if (opdis == 2) {
            cmat <- var(bvec[, 2:p1] - apply(bvec[, 2:p1], 2, 
                mean) + temp[2:p1])
            dis <- mahalanobis(temp2, temp[2:p1], cmat)
        }
        om.pval <- sum((dis[nboot + 1] <= dis[1:nboot]))/nboot
    }
    nval <- length(y)
    if (nval < 20) 
        nval <- 20
    if (nval > 60) 
        nval <- 60
    adj.pval <- om.pval/2 + (om.pval - om.pval/2) * (nval - 20)/40
    if (ncol(x) == 2 && plotit) {
        plot(bvec[, 2], bvec[, 3], xlab = "Slope 1", ylab = "Slope 2")
        temp.dis <- order(dis[1:nboot])
        ic <- round((1 - alpha) * nboot)
        xx <- bvec[temp.dis[1:ic], 2:3]
        xord <- order(xx[, 1])
        xx <- xx[xord, ]
        temp <- chull(xx)
        lines(xx[temp, ])
        lines(xx[c(temp[1], temp[length(temp)]), ])
    }
    list(output = output, om.pval = om.pval, adj.om.pval = adj.pval)
  }

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