cmanova:

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

Usage

cmanova(J, K, x, grp = c(1:JK), JK = J * K)

Arguments

J
K
x
grp
JK

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:JK), JK = J * K) 
{
    x = elimna(x)
    if (is.matrix(x)) 
        x <- listm(x)
    xx <- list()
    nvec <- NA
    jk <- 0
    for (j in 1:J) {
        for (k in 1:K) {
            jk <- jk + 1
            xx[[jk]] <- x[[grp[jk]]]
            if (k == 1) 
                nvec[j] <- length(xx[[jk]])
        }
    }
    N <- sum(nvec)
    RVALL <- matrix(0, nrow = N, K)
    x <- xx
    jk <- 0
    rmean <- matrix(NA, nrow = J, ncol = K)
    for (j in 1:J) {
        RV <- matrix(0, nrow = nvec[j], ncol = K)
        jk <- jk + 1
        temp1 <- matrix(x[[jk]], ncol = 1)
        for (k in 2:K) {
            jk <- jk + 1
            temp1 <- cbind(temp1, x[[jk]])
        }
        X <- temp1
        if (j == 1) 
            XALL <- X
        if (j > 1) 
            XALL <- rbind(XALL, X)
        n <- nvec[j]
        for (i in 1:n) {
            for (ii in 1:n) {
                temp3 <- sqrt(sum((X[i, ] - X[ii, ])^2))
                if (temp3 != 0) 
                  RV[i, ] <- RV[i, ] + (X[i, ] - X[ii, ])/temp3
            }
            RV[i, ] <- RV[i, ]/nvec[j]
            if (j == 1 && i == 1) 
                sighat <- RV[i, ] %*% t(RV[i, ])
            if (j > 1 || i > 1) 
                sighat <- sighat + RV[i, ] %*% t(RV[i, ])
        }
    }
    for (i in 1:N) {
        for (ii in 1:N) {
            temp3 <- sqrt(sum((XALL[i, ] - XALL[ii, ])^2))
            if (temp3 != 0) 
                RVALL[i, ] <- RVALL[i, ] + (XALL[i, ] - XALL[ii, 
                  ])/temp3
        }
        RVALL[i, ] <- RVALL[i, ]/N
    }
    bot <- 1 - nvec[1]
    top <- 0
    for (j in 1:J) {
        bot <- bot + nvec[j]
        top <- top + nvec[j]
        flag <- c(bot:top)
        rmean[j, ] <- apply(RVALL[flag, ], 2, mean)
    }
    sighat <- sighat/(N - J)
    shatinv <- solve(sighat)
    KW <- 0
    for (j in 1:J) {
        KW <- KW + nvec[j] * t(rmean[j, ]) %*% shatinv %*% rmean[j, 
            ]
    }
    df <- K * (J - 1)
    sig.level <- 1 - pchisq(KW, df)
    list(test.stat = KW[1, 1], df = df, p.value = sig.level)
  }

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