[R] Kendall's tau-c

[Marc Schwartz]

On Tue, 2006-09-12 at 15:02 +0200, Knut Wenzig wrote:
> Hello,
> 
> I can't find a package which calculates Kendall's tau-c. There is the
> package Kendall, but it only calcuates Kendall's tau-b.
> 
> Here is the example from
> ttp://www2.chass.ncsu.edu/garson/pa765/assocordinal.htm.
> 
> cityriots <- data.frame(citysize=c(1,1,2,2,3,3),
> riotsize=c(1,2,1,2,1,2), weight=c(4,2,2,3,0,4))
> cityriots <- data.frame(lapply(cityriots,function(x)
> rep(x,cityriots$weight)))
> xtabs(~ riotsize+citysize,cityriots)
> 
> tau-c should be .57.
> 
> Do you have a hint?
> 
> Best regards
> 
> Knut Wenzig


Here is some code:

# Calculate CONcordant Pairs in a table
# cycle through x[r, c] and multiply by
# sum(x elements below and to the right of x[r, c])
# x = table
concordant <- function(x)
{
  x <- matrix(as.numeric(x), dim(x))
  
  # get sum(matrix values > r AND > c)
  # for each matrix[r, c]
  mat.lr <- function(r, c)
  { 
    lr <- x[(r.x > r) & (c.x > c)]
    sum(lr)
  }

  # get row and column index for each
  # matrix element
  r.x <- row(x)
  c.x <- col(x)

  # return the sum of each matrix[r, c] * sums
  # using mapply to sequence thru each matrix[r, c]
  sum(x * mapply(mat.lr, r = r.x, c = c.x))
}

# Calculate DIScordant Pairs in a table
# cycle through x[r, c] and multiply by
# sum(x elements below and to the left of x[r, c])
# x = table
discordant <- function(x)
{
  x <- matrix(as.numeric(x), dim(x))
  
  # get sum(matrix values > r AND < c)
  # for each matrix[r, c]
  mat.ll <- function(r, c)
  { 
    ll <- x[(r.x > r) & (c.x < c)]
    sum(ll)
  }

  # get row and column index for each
  # matrix element
  r.x <- row(x)
  c.x <- col(x)

  # return the sum of each matrix[r, c] * sums
  # using mapply to sequence thru each matrix[r, c]
  sum(x * mapply(mat.ll, r = r.x, c = c.x))
}


# Calculate Kendall-Stuart Tau-c
# x = table
calc.KSTc <- function(x)
{
  x <- matrix(as.numeric(x), dim(x))
  
  c <- concordant(x)
  d <- discordant(x)
  m <- min(dim(x))
  n <- sum(x)

  KSTc <- (m * 2 * (c - d)) / ((n ^ 2) * (m - 1))

  KSTc
}


> calc.KSTc(with(cityriots, table(riotsize, citysize)))
[1] 0.5688889


The above code, along with other such measures, will eventually find its
way into the CrossTable() function in the gmodels CRAN package when time
permits (which seems to be in short supply of late...)

HTH,

Marc Schwartz