rfanova:

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

Usage

rfanova(x, grp = 0)

Arguments

x
grp

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, grp = 0) 
{
    library(MASS)
    if (!is.list(x)) 
        x <- listm(x)
    chk = tlist(x)
    if (chk != 0) 
        print("Warning: tied values detected")
    xall <- NA
    if (sum(grp) == 0) {
        J <- length(x)
        grp <- c(1:J)
    }
    if (sum(grp) > 0) 
        J <- length(grp)
    nval <- 1
    nrat <- 1
    nmax <- 0
    rbar <- 1
    mrbar <- 0
    for (j in grp) {
        temp <- x[[j]]
        temp <- temp[!is.na(temp)]
        nrat[j] <- (length(temp) - 1)/length(temp)
        nval[j] <- length(temp)
        if (j == grp[1]) 
            xall <- temp
        if (j != grp[1]) 
            xall <- c(xall, temp)
        if (length(temp) > nmax) 
            nmax <- length(temp)
    }
    pv <- array(NA, c(J, nmax, J))
    tv <- matrix(NA, J, nmax)
    rv <- matrix(0, J, nmax)
    for (i in 1:J) {
        data <- x[[i]]
        data <- data[!is.na(data)]
        for (j in 1:length(data)) {
            tempr <- data[j] - xall
            rv[i, j] <- length(tempr[tempr >= 0])
            for (l in 1:J) {
                templ <- x[[l]]
                templ <- templ[!is.na(templ)]
                temp <- data[j] - templ
                pv[i, j, l] <- length(temp[temp >= 0])
            }
            tv[i, j] <- sum(pv[i, j, ]) - pv[i, j, i]
        }
        rbar[i] <- sum(rv[i, ])/nval[i]
        mrbar <- mrbar + sum(rv[i, ])
    }
    amat <- matrix(0, J, J)
    for (i in 1:J) {
        temptv <- tv[i, ]
        temptv <- temptv[!is.na(temptv)]
        amat[i, i] <- (length(temptv) - 1) * var(temptv)
        for (l in 1:J) {
            tempp <- pv[l, , i]
            tempp <- tempp[!is.na(tempp)]
            if (l != i) {
                amat[i, i] <- amat[i, i] + (length(tempp) - 1) * 
                  var(tempp)
            }
        }
        for (j in 1:J) {
            if (j > i) {
                for (l in 1:J) {
                  temp1 <- pv[l, , i]
                  temp2 <- pv[l, , j]
                  temp1 <- temp1[!is.na(temp1)]
                  temp2 <- temp2[!is.na(temp2)]
                  if (i != l && l != j) 
                    amat[i, j] <- amat[i, j] + (length(temp1) - 
                      1) * var(temp1, temp2)
                }
                temp1 <- pv[i, , j]
                temp2 <- tv[i, ]
                temp1 <- temp1[!is.na(temp1)]
                temp2 <- temp2[!is.na(temp2)]
                amat[i, j] <- amat[i, j] - (length(temp1) - 1) * 
                  var(temp1, temp2)
                temp1 <- pv[j, , i]
                temp2 <- tv[j, ]
                temp1 <- temp1[!is.na(temp1)]
                temp2 <- temp2[!is.na(temp2)]
                amat[i, j] <- amat[i, j] - (length(temp1) - 1) * 
                  var(temp1, temp2)
            }
            amat[j, i] <- amat[i, j]
        }
    }
    N <- sum(nval)
    amat <- amat/N^3
    amati <- ginv(amat)
    uvec <- 1
    mrbar <- mrbar/N
    for (i in 1:J) uvec[i] <- nval[i] * (rbar[i] - mrbar)/(N * 
        (N + 1))
    testv <- N * prod(nrat) * uvec %*% amati %*% uvec
    test <- testv[1, 1]
    df <- J - 1
    siglevel <- 1 - pchisq(test, df)
    list(test = test, p.value = siglevel, df = df)
  }

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