bptd:

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

Usage

bptd(x, tr = 0.2, alpha = 0.05, con = 0, nboot = 599)

Arguments

x
tr
alpha
con
nboot

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, tr = 0.2, alpha = 0.05, con = 0, nboot = 599) 
{
    if (is.data.frame(x)) 
        x = as.matrix(x)
    if (!is.list(x) && !is.matrix(x)) 
        stop("Data must be stored in a matrix or in list mode.")
    if (is.list(x)) {
        if (is.matrix(con)) {
            if (length(x) != nrow(con)) 
                stop("The number of rows in con is not equal to the number of groups.")
        }
    }
    if (is.list(x)) {
        mat <- matrix(0, length(x[[1]]), length(x))
        for (j in 1:length(x)) mat[, j] <- x[[j]]
    }
    if (is.matrix(x)) 
        mat = x
    J <- ncol(mat)
    Jm <- J - 1
    if (sum(con^2) == 0) {
        d <- (J^2 - J)/2
        con <- matrix(0, J, d)
        id <- 0
        for (j in 1:Jm) {
            jp <- j + 1
            for (k in jp:J) {
                id <- id + 1
                con[j, id] <- 1
                con[k, id] <- 0 - 1
            }
        }
    }
    if (is.matrix(x)) {
        if (ncol(x) != nrow(con)) 
            stop("The number of rows in con is not equal to the number of groups.")
        mat <- x
    }
    if (sum(is.na(mat) >= 1)) 
        stop("Missing values are not allowed.")
    J <- ncol(mat)
    connum <- ncol(con)
    bvec <- matrix(0, connum, nboot)
    set.seed(2)
    xcen <- matrix(0, nrow(mat), ncol(mat))
    xbars <- matrix(0, nboot, ncol(mat))
    psihat <- matrix(0, connum, nboot)
    print("Taking bootstrap samples. Please wait.")
    data <- matrix(sample(nrow(xcen), size = nrow(mat) * nboot, 
        replace = TRUE), nrow = nboot)
    for (j in 1:J) {
        xcen[, j] <- mat[, j] - mean(mat[, j], tr)
        xbars[, j] <- apply(data, 1, bptdmean, xcen[, j], tr)
    }
    for (ic in 1:connum) {
        paste("Working on contrast number", ic)
        bvec[ic, ] <- apply(data, 1, bptdsub, xcen, tr, con[, 
            ic])
        psihat[ic, ] <- apply(xbars, 1, bptdpsi, con[, ic])
    }
    bvec <- psihat/sqrt(bvec)
    bvec <- abs(bvec)
    icrit <- round((1 - alpha) * nboot)
    critvec <- apply(bvec, 2, max)
    critvec <- sort(critvec)
    crit <- critvec[icrit]
    psihat <- matrix(0, connum, 4)
    dimnames(psihat) <- list(NULL, c("con.num", "psihat", "ci.lower", 
        "ci.upper"))
    test <- matrix(NA, connum, 3)
    dimnames(test) <- list(NULL, c("con.num", "test", "se"))
    isub <- c(1:nrow(mat))
    tmeans <- apply(mat, 2, mean, trim = tr)
    sqse <- 1
    psi <- 1
    for (ic in 1:ncol(con)) {
        sqse[ic] <- bptdsub(isub, mat, tr, con[, ic])
        psi[ic] <- sum(con[, ic] * tmeans)
        psihat[ic, 1] <- ic
        psihat[ic, 2] <- psi[ic]
        psihat[ic, 3] <- psi[ic] - crit * sqrt(sqse[ic])
        psihat[ic, 4] <- psi[ic] + crit * sqrt(sqse[ic])
        test[ic, 1] <- ic
        test[ic, 2] <- psi[ic]/sqrt(sqse[ic])
        test[ic, 3] <- sqrt(sqse[ic])
    }
    list(test = test, psihat = psihat, crit = crit, con = con)
  }

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