whitney

whitney <- function(x, groups, weights = NULL)  {

  f <- ! is.na(x) & ! is.na(groups)
  x <- x[f]
  groups <- groups[f]
  weights <- weights[f]
  lv <- levels(groups)

  sort.x <- sort(x, index.return=TRUE)
  ranks <- 1:length(x)
  unique.x <- unique(sort.x$x)

  if (! is.null(weights)) {
    require(Hmisc)
    ranks <- wtd.rank(x, weights=weights)
  } else {
    for (i in unique.x) {
      f <- sort.x$x == i
      rf <- ranks[f]
      ranks[f] <- rep(mean(rf), length(rf))
    }
    new.index <- sort(sort.x$ix, index.return=TRUE)$ix
    ranks <- ranks[new.index]
  }


  g1 <- ranks[groups == lv[1]]
  g2 <- ranks[groups == lv[2]]

  if (! is.null(weights)) {
    wt1 <- weights[groups == lv[1]]
    wt2 <- weights[groups == lv[2]]

    n1 <- sum(wt1)
    n2 <- sum(wt2)
    rank1 <- sum(wt1*g1)
    rank2 <- sum(wt2*g2)
  } else {
    n1 <- length(g1)
    n2 <- length(g2)
    rank1 <- sum(g1)
    rank2 <- sum(g2)
  }

  u1 <- n1*n2 + (n1*(n1+1))/2 - rank1
  u2 <- n1*n2 + (n2*(n2+1))/2 - rank2
  u <- min(u1,u2)

  if (u > (n1*n2)/2) {
    u <- n1*n2 - u
    w <- rank1
  } else {
    w <- rank2
  }

  mu <- (n1*n2)/2
  sigma <- sqrt((n1*n2*(n1+n2+1))/12)

  z <- (u - mu)/sigma

  if (n1 < 50 & n2 < 50) {
    p <- pwilcox(u, n1, n2)*2
  } else {
    p <- pnorm(-abs(z))*2
  }

  result <- list(u = u, w = w, z = z, "p-value" = p)
  result$names <- levels(groups)
  result$ranksum <- c(rank1, rank2)
  result$n <- c(n1, n2)

  class(result) <- "whitney"

  return(result)
}

print.whitney <- function(x) {
  tab <- data.frame("Grupos" = x$names,
                    "n" = x$n,
                    "Suma de rangos" = x$ranksum)
  print(tab)
  cat("\n")
  p <- format(c(x$u, x$w, x$z))
  cat(paste("U de Mann-Whitney: ", p[1], "\n"))
  cat(paste("    W de Wilcoxon: ", p[2], "\n"))
  cat(paste("                Z: ", p[3], "\n"))
  cat("\n")
  cat(paste("p-value: ", x[["p-value"]], "\n"))
  invisible(x)
}