bwrank:

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

Usage

bwrank(J, K, x, grp = c(1:p), p = J * K)

Arguments

J
K
x
grp
p

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, grp = c(1:p), p = J * K) 
{
    if (is.data.frame(x)) 
        data = as.matrix(x)
    if (is.matrix(x)) 
        x <- listm(x)
    x = x[grp]
    xx <- list()
    nvec <- NA
    alldat <- NA
    klow <- 1 - K
    kup <- 0
    iall = 0
    for (j in 1:J) {
        klow <- klow + K
        kup <- kup + K
        mtemp = elimna(matl(x[klow:kup]))
        for (k in 1:K) {
            iall = iall + 1
            xx[[iall]] = mtemp[, k]
        }
    }
    for (j in 1:p) {
        alldat <- c(alldat, xx[[j]])
        nvec[j] <- length(xx[[j]])
    }
    nmat <- matrix(nvec, J, K, byrow = T)
    for (j in 1:J) {
        if (var(nmat[j, ]) != 0) {
            warning("Number of observations among dependent groups for level", 
                paste(j), " of Factor A are unequal")
            print("Matrix of sample sizes:")
            print(nmat)
        }
    }
    if (sum(is.na(alldat[2:length(alldat)]) > 0)) 
        stop("Missing values not allowed")
    rval <- rank(alldat[2:length(alldat)])
    rdd <- mean(rval)
    xr <- list()
    il <- 1 - nvec[1]
    iu <- 0
    for (j in 1:p) {
        il <- il + nvec[j]
        iu <- iu + nvec[j]
        xr[[j]] <- rval[il:iu]
    }
    v <- matrix(0, p, p)
    Ja <- matrix(1, J, J)
    Ia <- diag(1, J)
    Pa <- Ia - Ja/J
    Jb <- matrix(1, K, K)
    Ib <- diag(1, K)
    Pb <- Ib - Jb/K
    cona <- kron(Pa, Ib)
    conb <- kron(Ia, Pb)
    conab <- kron(Pa, Pb)
    for (k in 1:K) {
        temp <- x[[k]]
        bigm <- matrix(temp, ncol = 1)
        jk <- k
        for (j in 2:J) {
            jk <- jk + K
            tempc <- matrix(x[[jk]], ncol = 1)
            bigm <- rbind(bigm, tempc)
            temp <- c(temp, x[[jk]])
        }
    }
    N <- length(temp)
    rbbd <- NA
    for (k in 1:K) {
        bigm <- xr[[k]]
        jk <- k
        for (j in 2:J) {
            jk <- jk + K
            bigm <- c(bigm, xr[[jk]])
        }
    }
    rbjk <- matrix(NA, nrow = J, ncol = K)
    ic <- 0
    for (j in 1:J) {
        for (k in 1:K) {
            ic <- ic + 1
            rbjk[j, k] <- mean(xr[[ic]])
        }
    }
    for (k in 1:K) rbbd[k] <- mean(rbjk[, k])
    rbj <- 1
    sigv <- 0
    njsam <- 0
    icc <- 1 - K
    ivec <- c(1:K) - K
    for (j in 1:J) {
        icc <- icc + K
        ivec <- ivec + K
        temp <- xr[ivec[1]:ivec[K]]
        temp <- matl(temp)
        tempv <- apply(temp, 1, mean)
        njsam[j] <- nvec[icc]
        rbj[j] <- mean(rbjk[j, ])
        sigv[j] <- var(tempv)
    }
    nv <- sum(njsam)
    phat <- (rbjk - 0.5)/(nv * K)
    sv2 <- sum(sigv/njsam)
    uv <- sum((sigv/njsam)^2)
    dv <- sum((sigv/njsam)^2/(njsam - 1))
    testA <- J * var(rbj)/sv2
    klow <- 1 - K
    kup <- 0
    sv <- matrix(0, nrow = K, ncol = K)
    rvk <- NA
    for (j in 1:J) {
        klow <- klow + K
        kup <- kup + K
        sel <- c(klow:kup)
        m <- matl(xr[klow:kup])
        m <- elimna(m)
        xx <- listm(m)
        xx <- listm(m)
        vsub <- nv * var(m)/(nv * K * nv * K * njsam[j])
        v[sel, sel] <- vsub
        sv <- sv + vsub
    }
    sv <- sv/J^2
    testB <- nv/(nv * K * nv * K * sum(diag(Pb %*% sv))) * sum((rbbd - 
        mean(rbbd))^2)
    testAB <- 0
    for (j in 1:J) {
        for (k in 1:K) {
            testAB <- testAB + (rbjk[j, k] - rbj[j] - rbbd[k] + 
                rdd)^2
        }
    }
    testAB <- nv * testAB/(nv * K * nv * K * sum(diag(conab %*% 
        v)))
    nu1B <- (sum(diag(Pb %*% sv)))^2/sum((diag(Pb %*% sv %*% 
        Pb %*% sv)))
    nu1A <- (J - 1)^2/(1 + J * (J - 2) * uv/sv2^2)
    nu2A <- sv2^2/dv
    nu1AB <- (sum(diag(conab %*% v)))^2/sum(diag(conab %*% v %*% 
        conab %*% v))
    sig.A <- 1 - pf(testA, nu1A, nu2A)
    sig.B <- 1 - pf(testB, nu1B, 1e+06)
    sig.AB <- 1 - pf(testAB, nu1AB, 1e+06)
    list(test.A = testA, sig.A = sig.A, test.B = testB, sig.B = sig.B, 
        test.AB = testAB, sig.AB = sig.AB, avg.ranks = rbjk, 
        rel.effects = phat)
  }

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