ridit.R

ridit <- function(x, category.var, reference.var = NULL, alpha=0.05) {
  if (missing(category.var) || is.null(category.var) || !is.character(category.var))
    stop("Must specify category variable")
  if (!is.factor(x[[category.var]])) stop("category variable must be a factor")
  if (is.null(reference.var)) stop("Must specify reference group")

  ## make local copy
  tmp <- x
  tmp[, c(category.var)] <- list(NULL)
  tmp <- apply(tmp, 2, function(x) x/sum(x))

  ## compute riddit
  do.ridit <- function(props) { ## props should be in order of levels (highest to lowest)
    r <- rep(-1, length(props))
    for (i in 1:length(props)) {
      if (i == length(props)) vals <- 0
      else vals <- props[(i+1):length(props)]
      r[i] <- sum(vals) + 0.5*props[i]
    }
    return(r)
  }
  ridit.vals <- apply(tmp, 2, do.ridit)

  ridit.vals <- cbind(x[[category.var]], data.frame(ridit.vals))
  names(ridit.vals)[1] <- category.var

  ## mean ridit for each group wrt specified reference
  ref.ridit <- ridit.vals[[reference.var]]
  mean.ridits <- apply(tmp, 2, function(col) sum(col*ref.ridit))

  
  ## CI's
  nref <- sum(x[[reference.var]])
  k <- length(mean.ridits)-1
  K <- k*(k+1)/2
  B <- qnorm(1-alpha/K/2)
  cis <- do.call(rbind, lapply(names(mean.ridits), function(colname) {
    n <- sum(x[[colname]])
    r <- mean.ridits[[colname]]
    f <- B*sqrt(nref+n)/(2*sqrt(3*nref*n))
    return(data.frame(low=r-f, high=r+f))
  }))
  cis <- cbind(group=names(mean.ridits), cis)

  ## remove reference var
  ridit.vals[, c(reference.var)] <- list(NULL)
  mean.ridits <- mean.ridits[ names(mean.ridits) != reference.var ]
  cis <- subset(cis, group != reference.var)
  list(category.ridit=ridit.vals, mean.ridit=mean.ridits, ci=cis)
}