pbad2way:

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

Usage

pbad2way(J, K, x, est = onestep, conall = TRUE, alpha = 0.05, nboot = 2000, grp = NA, op = FALSE, pro.dis = FALSE, MM = FALSE, ...)

Arguments

J
K
x
est
conall
alpha
nboot
grp
op
pro.dis
MM
...

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 = onestep, conall = TRUE, alpha = 0.05, 
    nboot = 2000, grp = NA, op = FALSE, pro.dis = FALSE, MM = FALSE, 
    ...) 
{
    JK <- J * K
    if (is.matrix(x)) 
        x <- listm(x)
    if (!is.na(grp[1])) {
        yy <- x
        for (j in 1:length(grp)) x[[j]] <- yy[[grp[j]]]
    }
    if (!is.list(x)) 
        stop("Data must be stored in list mode or a matrix.")
    for (j in 1:JK) {
        xx <- x[[j]]
        x[[j]] <- xx[!is.na(xx)]
    }
    if (!conall) {
        ij <- matrix(c(rep(1, J)), 1, J)
        ik <- matrix(c(rep(1, K)), 1, K)
        jm1 <- J - 1
        cj <- diag(1, jm1, J)
        for (i in 1:jm1) cj[i, i + 1] <- 0 - 1
        km1 <- K - 1
        ck <- diag(1, km1, K)
        for (i in 1:km1) ck[i, i + 1] <- 0 - 1
        conA <- t(kron(cj, ik))
        conB <- t(kron(ij, ck))
        conAB <- t(kron(cj, ck))
        conAB <- t(kron(abs(cj), ck))
    }
    if (conall) {
        temp <- con2way(J, K)
        conA <- temp$conA
        conB <- temp$conB
        conAB <- temp$conAB
    }
    ncon <- max(nrow(conA), nrow(conB), nrow(conAB))
    if (JK != length(x)) 
        warning("The number of groups does not match the number of contrast coefficients.")
    if (!is.na(grp[1])) {
        xx <- list()
        for (i in 1:length(grp)) xx[[i]] <- x[[grp[i]]]
        x <- xx
    }
    mvec <- NA
    for (j in 1:JK) {
        temp <- x[[j]]
        temp <- temp[!is.na(temp)]
        x[[j]] <- temp
        mvec[j] <- est(temp, ...)
    }
    bvec <- matrix(NA, nrow = JK, ncol = nboot)
    set.seed(2)
    print("Taking bootstrap samples. Please wait.")
    for (j in 1:JK) {
        data <- matrix(sample(x[[j]], size = length(x[[j]]) * 
            nboot, replace = TRUE), nrow = nboot)
        bvec[j, ] <- apply(data, 1, est, ...)
    }
    bconA <- t(conA) %*% bvec
    tvecA <- t(conA) %*% mvec
    tvecA <- tvecA[, 1]
    tempcenA <- apply(bconA, 1, mean)
    veczA <- rep(0, ncol(conA))
    bconA <- t(bconA)
    smatA <- var(bconA - tempcenA + tvecA)
    bconA <- rbind(bconA, veczA)
    if (!pro.dis) {
        if (!op) 
            dv <- mahalanobis(bconA, tvecA, smatA)
        if (op) {
            dv <- out(bconA)$dis
        }
    }
    if (pro.dis) 
        dv = pdis(bconA, MM = MM)
    bplus <- nboot + 1
    sig.levelA <- 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
    bconB <- t(conB) %*% bvec
    tvecB <- t(conB) %*% mvec
    tvecB <- tvecB[, 1]
    tempcenB <- apply(bconB, 1, mean)
    veczB <- rep(0, ncol(conB))
    bconB <- t(bconB)
    smatB <- var(bconB - tempcenB + tvecB)
    bconB <- rbind(bconB, veczB)
    if (!pro.dis) {
        if (!op) 
            dv <- mahalanobis(bconB, tvecB, smatB)
        if (op) {
            dv <- out(bconA)$dis
        }
    }
    if (pro.dis) 
        dv = pdis(bconB, MM = MM)
    sig.levelB <- 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
    bconAB <- t(conAB) %*% bvec
    tvecAB <- t(conAB) %*% mvec
    tvecAB <- tvecAB[, 1]
    tempcenAB <- apply(bconAB, 1, mean)
    veczAB <- rep(0, ncol(conAB))
    bconAB <- t(bconAB)
    smatAB <- var(bconAB - tempcenAB + tvecAB)
    bconAB <- rbind(bconAB, veczAB)
    if (!pro.dis) {
        if (!op) 
            dv <- mahalanobis(bconAB, tvecAB, smatAB)
        if (op) {
            dv <- out(bconAB)$dis
        }
    }
    if (pro.dis) 
        dv = pdis(bconAB, MM = MM)
    sig.levelAB <- 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
    list(sig.levelA = sig.levelA, sig.levelB = sig.levelB, sig.levelAB = sig.levelAB, 
        conA = conA, conB = conB, conAB = conAB)
  }

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