# Comparing Kendall's Tau values

I have a single ordinal outcome called RATING, and a series of continuous covariates (P1, P2, P3, etc). I have used Kendall's Tau statistic to assess the concordance/correlation between RATING and each covariate as per below.

#load Kendall package
library(Kendall)

#create example data
RATING<-sample(1:4,200,replace=TRUE)
P1<-RATING+rnorm(200,0,2.5)
P2<-RATING+rnorm(200,0,3)
P3<-RATING+rnorm(200,0,3.5)

#Calculate Kendall's Tau for each predictor
Kendall(RATING,P1)
Kendall(RATING,P2)
Kendall(RATING,P3)


The resulting Tau statistics will differ in size and significance obviously. My question though is whether there a way to determine which are significantly different from others? For example, can I use the results of these tests to say P1 is significantly more strongly correlated with RATING than P3, or that P3 is significantly less correlated with RATING than all other covariates?

I've found discussion of a method to compare Pearson's R values here How to compare the strength of two Pearson correlations? but I believe it assumes the compared R values are from independent samples, its not clear if it can be applied to Kendall's tau values, and its not clear whether its suitable for multiple comparisons.

Any suggestions greatly appreciated

Here's a rough and ready approach using bootstrapping, borrowing a method from Guillame Rousselet.

First I'll replicate your data, modifying P1 slightly such that it will have a strong relationship with the Rating just to illustrate. I'm using the the cor function from the default stats package with the method = "kendall" argument throughout.

set.seed(6000)
Rating <- sample(1:4, 200, replace=TRUE)
df <- data.frame( Rating,
P1 = Rating + rnorm(200, 0, .8),
P2 = Rating + rnorm(200, 0, 3.5),
P3 = Rating + rnorm(200, 0, 3))

#Check the 3 Kendall Tau values
sapply(2:4, function(x) cor(df[ , 1], df[ , x], method = "kendall"))
# [1] 0.6535829 0.1738252 0.2347800


So far so good.

Now we'll want to iteratively resample observations in the dataset, calculate the three overlapping Kendall's $$\tau$$ values, and take the differences between each pair ($$\tau_i-\tau_j)$$. Determining whether the bootstrapped distributions of these differences overlap zero provides your test.

library(boot)
library(broom)

#Create Bootstrap Function
compKendall <- function(data, indices){
df <- data[indices, ] #Resample over rows

#Calculate all 3 overlapping correlations
ks <- sapply(2:4, function(x) cor(df[ , 1], df[ , x], method = "kendall"))

#Arrange pairwise combinations of correlations
kPairs <- combn(ks, 2)

#Return the differences
return( kPairs[1 , ] - kPairs[2 , ] )
}

#Run the bootstrap with 1000 resamples
kDiff <- boot(df, statistic = compKendall, R = 1000)

#Use broom to get a neat dataframe of differences w/ 95% CI's
broom::tidy(kDiff, conf.int = TRUE, conf.method = "perc")

##  A tibble: 3 x 5
#   statistic     bias std.error conf.low conf.high
#      <dbl>    <dbl>     <dbl>    <dbl>     <dbl>
# 1    0.480   0.00102    0.0533    0.378    0.591
# 2    0.419  -0.00308    0.0525    0.312    0.518
# 3   -0.0610 -0.00410    0.0694   -0.205    0.0720


The dataframe shows the estimated differences between Kendall's $$\tau$$ values, organised as $$\tau$$ number 1 - 2, 1 - 3, and 2 - 3.

As expected, the first correlation ($$\tau_1 = .65$$) appears larger than the other two - CI's for the differences exceed zero. However we couldn't conclude there is a significant difference between the second ($$\tau_2 = .17$$) and third ($$\tau_3 = .23$$) overlapping correlations.

• Thanks so much awhug. I will try this out. I don't see anything in the Rousselet link about adjustments for cases where there are many covariates (eg P1, P2, ... P27) and so multiple comparisons. Do you have any thoughts on the necessity of adjusting this method for that circumstance? – Terry B Jul 25 at 7:14
• That is a good question - multiplicity adjustments for CI's are hard to come by. Perhaps you could treat each difference estimate divided by its SE as the test statistic Z. You might save the broom::tidy dataframe as results and determine the p-value for each by 2 * pnorm( -abs(results$statistic / results$std.error) ). From there you could apply a p.adjust() of your choice to the results. But I don't feel 100% confident recommending this without checking the assumptions of the bootstrap. – awhug Jul 29 at 6:59
• Missed that last comment when it came in, but thankyou, I will investigate – Terry B Nov 9 at 12:06