spmcpa:

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

Usage

spmcpa(J, K, x, est = tmean, JK = J * K, grp = c(1:JK), con = 0, avg = FALSE, alpha = 0.05, nboot = NA, pr = TRUE, ...)

Arguments

J
K
x
est
JK
grp
con
avg
alpha
nboot
pr
...

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 (J, K, x, est = tmean, JK = J * K, grp = c(1:JK), con = 0, 
    avg = FALSE, alpha = 0.05, nboot = NA, pr = TRUE, ...) 
{
    if (is.matrix(x)) {
        y <- list()
        for (j in 1:ncol(x)) y[[j]] <- x[, j]
        x <- y
    }
    if (pr) 
        print("As of Sept. 2005, est defaults to tmean")
    JK <- J * K
    data <- list()
    for (j in 1:length(x)) {
        data[[j]] <- x[[grp[j]]]
    }
    x <- data
    jp <- 1 - K
    kv <- 0
    kv2 <- 0
    for (j in 1:J) {
        jp <- jp + K
        xmat <- matrix(NA, ncol = K, nrow = length(x[[jp]]))
        for (k in 1:K) {
            kv <- kv + 1
            xmat[, k] <- x[[kv]]
        }
        xmat <- elimna(xmat)
        for (k in 1:K) {
            kv2 <- kv2 + 1
            x[[kv2]] <- xmat[, k]
        }
    }
    xx <- x
    set.seed(2)
    nvec <- NA
    jp <- 1 - K
    for (j in 1:J) {
        jp <- jp + K
        nvec[j] <- length(x[[jp]])
    }
    if (avg) {
        d <- (J^2 - J)/2
        con <- matrix(0, J, d)
        id <- 0
        Jm <- J - 1
        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 (!avg) {
        MJK <- K * (J^2 - J)/2
        JK <- J * K
        MJ <- (J^2 - J)/2
        cont <- matrix(0, nrow = J, ncol = MJ)
        ic <- 0
        for (j in 1:J) {
            for (jj in 1:J) {
                if (j < jj) {
                  ic <- ic + 1
                  cont[j, ic] <- 1
                  cont[jj, ic] <- 0 - 1
                }
            }
        }
        tempv <- matrix(0, nrow = K - 1, ncol = MJ)
        con1 <- rbind(cont[1, ], tempv)
        for (j in 2:J) {
            con2 <- rbind(cont[j, ], tempv)
            con1 <- rbind(con1, con2)
        }
        con <- con1
        if (K > 1) {
            for (k in 2:K) {
                con1 <- push(con1)
                con <- cbind(con, con1)
            }
        }
    }
    d <- ncol(con)
    if (is.na(nboot)) {
        if (d <= 4) 
            nboot <- 1000
        if (d > 4) 
            nboot <- 5000
    }
    bloc <- matrix(NA, nrow = J, ncol = nboot)
    print("Taking bootstrap samples. Please wait.")
    mvec <- NA
    ik <- 0
    for (j in 1:J) {
        paste("Working on level ", j, " of Factor A")
        x <- matrix(NA, nrow = nvec[j], ncol = K)
        for (k in 1:K) {
            ik <- ik + 1
            x[, k] <- xx[[ik]]
            if (!avg) 
                mvec[ik] <- est(xx[[ik]], ...)
        }
        tempv <- apply(x, 2, est, ...)
        data <- matrix(sample(nvec[j], size = nvec[j] * nboot, 
            replace = TRUE), nrow = nboot)
        bvec <- matrix(NA, ncol = K, nrow = nboot)
        mat <- listm(x)
        for (k in 1:K) {
            temp <- x[, k]
            bvec[, k] <- apply(data, 1, rmanogsub, temp, est, 
                ...)
        }
        if (avg) {
            mvec[j] <- mean(tempv)
            bloc[j, ] <- apply(bvec, 1, mean)
        }
        if (!avg) {
            if (j == 1) 
                bloc <- bvec
            if (j > 1) 
                bloc <- cbind(bloc, bvec)
        }
    }
    if (avg) 
        bloc <- t(bloc)
    connum <- d
    psihat <- matrix(0, connum, nboot)
    test <- 1
    for (ic in 1:connum) {
        psihat[ic, ] <- apply(bloc, 1, bptdpsi, con[, ic])
        test[ic] <- (sum(psihat[ic, ] > 0) + 0.5 * sum(psihat[ic, 
            ] == 0))/nboot
        test[ic] <- min(test[ic], 1 - test[ic])
    }
    ncon <- ncol(con)
    if (alpha == 0.05) {
        dvec <- c(0.025, 0.025, 0.0169, 0.0127, 0.0102, 0.00851, 
            0.0073, 0.00639, 0.00568, 0.00511)
        if (ncon > 10) {
            avec <- 0.05/c(11:ncon)
            dvec <- c(dvec, avec)
        }
    }
    if (alpha == 0.01) {
        dvec <- c(0.005, 0.005, 0.00334, 0.00251, 0.00201, 0.00167, 
            0.00143, 0.00126, 0.00112, 0.00101)
        if (ncon > 10) {
            avec <- 0.01/c(11:ncon)
            dvec <- c(dvec, avec)
        }
    }
    if (alpha != 0.05 && alpha != 0.01) {
        dvec <- alpha/c(1:ncon)
        dvec[1] <- alpha/2
    }
    temp2 <- order(0 - test)
    ncon <- ncol(con)
    zvec <- dvec[1:ncon]
    sigvec <- (test[temp2] >= zvec)
    output <- matrix(0, connum, 6)
    dimnames(output) <- list(NULL, c("con.num", "psihat", "p.value", 
        "p.sig", "ci.lower", "ci.upper"))
    tmeans <- mvec
    psi <- 1
    output[temp2, 4] <- zvec
    for (ic in 1:ncol(con)) {
        output[ic, 2] <- sum(con[, ic] * tmeans)
        output[ic, 1] <- ic
        output[ic, 3] <- test[ic]
        temp <- sort(psihat[ic, ])
        temp3 <- round(output[ic, 4] * nboot) + 1
        icl <- round(dvec[ncon] * nboot) + 1
        icu <- nboot - (icl - 1)
        output[ic, 5] <- temp[icl]
        output[ic, 6] <- temp[icu]
    }
    output[, 3] <- 2 * output[, 3]
    output[, 4] <- 2 * output[, 4]
    num.sig <- sum(output[, 3] <= output[, 4])
    list(output = output, con = con, num.sig = num.sig)
  }

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