[SciPy-dev] Possible Error in Kendall's Tau (scipy.stats.stats.kendalltau)

[josef.pktd at gmail.com]

On Wed, Mar 18, 2009 at 9:29 AM, Bruce Southey <bsouthey at gmail.com> wrote:
> Almer S. Tigelaar wrote:
>> Hello,
>>
>> On Wed, 2009-03-18 at 13:11 +0100, Sturla Molden wrote:
>>
>>> So it seems Hollander & Wolfe and Minitab says 0.67, whereas Numerical
>>> Receipes says 1.0. Intuitively a vector correlation should be exactly
>>> correlated with itself, but I am inclined to trust Hollander & Wolfe
>>> more than Numerical Receipes.
>>>
>>
>> Ah, I was under the impression you already checked Hollander & Wolfe.
>> Anyway, it seems my initial interpretation was right then. Repeating the
>> formula here (augmented) for future reference:
>>
>> Kendall's tau-b (tie handling):
>> -------------------------------
>> Given two rankings R1 and R2, Kendall's tau-b is calculated by:
>>         t = (P - Q) / SQRT((P + Q + T) * (P + Q + U))
>> where P is the number of concordant pairs, Q the number of discordant
>> pairs, T the number of ties in R1 and U the number of ties in R2.
>> [Ties are always counted regardless of whether they occur for the same
>> pair in R1 and R2 or different pairs]
>> -------------------------------
>>
>> Some tests I ran today with the R implementation of Kendall's Tau(-a)
>> and the original implementation in SciPy.stats.stats (Kendall's Tau-b)
>> seem to suggests that if we do NOT count ties on the same pair (the
>> current situation in SciPy.stats.stats) effectively Kendall's Tau-b
>> gives the same outcomes as Kendall's Tau-a for about 36 test cases.
>>
>> This seems to suggest that Kendall's Tau-b (tie correction) in SciPy as
>> it is behaves like Kendall's Tau-a (no tie correction), possibly because
>> of leaving out ties on identical pairs in T and U above.
>>
>> I unfortunately do not have the time to mathematically prove (or
>> disprove) the equivalence of Kendall's Tau-a and the current SciPy
>> implementation right now, but I thought I'd be useful to mention these
>> test results.
>>
>>
>
> Hi,
> This link might be useful as it has worked examples:
> http://faculty.chass.ncsu.edu/garson/PA765/assocordinal.htm
>
> I find that the SAS documentation for Proc Freq very useful:
> http://support.sas.com/onlinedoc/913/getDoc/en/statug.hlp/freq_index.htm
> http://support.sas.com/onlinedoc/913/getDoc/en/statug.hlp/freq_sect20.htm
>
> Also, does these implementation depend on the type of array?
> It would be great to have a single function that accepts an array,
> masked array or an object that can be converted into an array.
>
> Finally, this measure assumes ordinal data but there is no type checking
> done in the Scipy function.
>
> Bruce



I'm giving up, this takes too much time:

the sas reference by Bruce hat a reference to a book by Agresti, where
I finally found a clear formal definition and it excludes matching
ties in the in the denominator:

http://books.google.ca/books?id=hpEzw4T0sPUC&dq=Agresti,+Alan+(1996).+Introduction+to+categorical+data+analysis&printsec=frontcover&source=bn&hl=en&ei=hADBSamzBpKWMrTttJoI&sa=X&oi=book_result&resnum=4&ct=result#PPA68,M1

categorical data:
sas and spss are using kendall's tau and the other association
measures for contingency tables, and sas says kendall tau-b is only
appropriate for categorical data.
I think this is not true in general, a saw other applications of
kendall's tau for continuous (normal) random variables, and since in
this case there are no ties all definitions for kendalls tau or gamma
produce identical results, corresponding to the probability and
correlation definition.

the definition that Bruce found:  Kendall's tau according to
Hollander, M.,  and D. A. Wolfe. 1999,
I think, corresponds to kendall's tau-a since there is no tie
correction in the denominator.

general principles:
I like the correlation interpretation that I mentioned,
In a reference on the wikipedia page, they gave a preference ordering
interpretation of kendall's tau, how strongly are two preference
orderings, e.g. rankings related. They refered to it only for strict
preferences, but the same should hold for weak preferences. If you
consider [1,1,2] as preference ranking, than two individuals with the
same weak ranking have exactly the same preferences, and the
correlation should be 1 and not 2/3

I don't know about kendall's tau-c because m seems to be specific to
contingency tables, while all other measures have a more general
interpretation. Similar in view of the general definitions, I don't
understand the talk about square or rectangular tables. But I don't
have a good intuition for contingency tables.

I haven't checked the comparison with R, but I would be very surprised
if stats.kendalltau corresponds to kendall's tau-a in R.

the version of kendall's tau-b that adds matching ties in the
denominator is wrong from my reading.

So, I prefer to leave stats.kendalltau as it is and maybe add on
option for tiehandling, so that also obtain kendall's tau-a and maybe
gamma (in the definition of sas and spss)

Josef

attached is my test script
-------------- next part --------------
import numpy as np
from scipy import stats
from numpy.testing import assert_equal

def kendalltaustat(x, y):
    '''calculate Kendall's tau-b correlation statistic
    this is just the (non-central) correlation of all pairwise rankings

    Notes
    -----
    tau-b = (P - Q) / sqrt(((P + Q + Tx - Txy)(P + Q + Ty - Txy)))
    equivalent to
    tau-b = (P - Q) / sqrt(((n*(n-1) - Tx)*((n*(n-1) - Ty)))
    total number of pairs (n over 2) = (n*(n-1)
    Tx number of ties in x
    Ty number of ties in y
    Txy number of pairs that are tied in x and in y
    Tx - Txy number of pairs only tied in x and not in y
    Ty - Txy number of pairs only tied in y and not in x
    P = sum_pairs( (x1-x2)*(y1-y2)>0)   concordant pairs
    Q = sum_pairs( (x1-x2)*(y1-y2)<0)   discordant pairs
    '''
    # not this is using 2 times the pairs
    # calculate indicators of all pairs
    ppos1 = np.sign((x[:,np.newaxis] - x)).astype(float).ravel()
    ppos2 = np.sign((y[:,np.newaxis] - y)).astype(float).ravel()
    #correlation coefficient without mean correction
    tau = np.dot(ppos1,ppos2) / np.sqrt(np.dot(ppos1,ppos1) * np.dot(ppos2,ppos2))
    return tau

def kendalltaustata(x, y):
    '''calculate Kendall's tau-a correlation statistic
    sum of rank indicator (-1,0,+1) divided by all pairs
    '''
    n = x.shape[0]
    tau = (np.sum(np.sign((x[:,np.newaxis]-x)*(y[:,np.newaxis]-y))))/(1.0*n*(n-1))
    return tau

def gammaassoc(x, y):
    ''' this is gamma correlation statistic
    this is just the (non-central) correlation of all pairwise rankings
    '''
    n = x.shape[0]
    conc = np.sign((x[:,np.newaxis]-x)*(y[:,np.newaxis]-y))  
    tau = np.sum(conc)/float(np.sum(np.abs(conc)))
    return tau

def kendalltaustatba(x, y):
    ''' calculate Kendall's tau-b correlation statistic with different
    tie handling, adding simultanous ties in x and y in denominator

    '''
    n = x.shape[0]
    conc = np.sign((x[:,np.newaxis]-x)*(y[:,np.newaxis]-y))
    strict = np.sum(np.abs(conc))
    tiesx = np.sum((x[:,np.newaxis]-x)==0) - n
    tiesy = np.sum((y[:,np.newaxis]-y)==0) - n
    tau = np.sum(conc)/np.float(np.sqrt((strict + tiesx) * (strict + tiesx)))
    return tau

x1a = np.array([0, 1, 3, 3, 4, 5, 5, 7, 8])
x1b = np.array([1, 3, 3, 4, 5, 5, 7, 8, 9])
x1c = np.array([1, 3, 3, 3, 5, 5, 7, 8, 9])
x1d = np.array([1, 3, 3, 3, 5, 5, 7, 8, 6])
x1 = np.array([1,1,2])
x2a = np.array([1,1,2,2])
x2b = np.array([1,2,3,4])

data = [(x1a, x1a),
        (x1a, x1b),
        (x1a, x1c),
        (x1a, 10-x1c),
        (x1, x1),
        (x2a, x2b),
        (x2b, x2b),
        (x2a, 3-x2a),
        (x2a, 3-x2b),
        ]

data = [(x1a, x1a, 1.0, 0.00017455009626808976), 
        (x1a, x1b, 0.94117647058823528, 0.0004116823127257346),
        (x1a, x1c, 0.93982554701579024, 0.00041964794448139537),
        (x1a, x1d, 0.81855773449762381, 0.0021244496819179748),
        (x1a, 10-x1c, -0.93982554701579024, 0.00041964794448139537),
        (x1, x1, 1.0, 0.11718509694604401),
        (x2a, x2b, 0.81649658092772615, 0.0960923383021903),
        (x2b, x2b, 1.0, 0.041540072431798185),
        (x2a, 3-x2a, -1.0, 0.041540072431798185),
        (x2a, 3-x2b, -0.81649658092772615, 0.0960923383021903),
        ]


for x,y, tc, pc in data:
    t1 = kendalltaustat(x, y)
    ts, ps = stats.kendalltau(x, y)
    ta = kendalltaustata(x, y)
    taa = gammaassoc(x, y)
    tba = kendalltaustatba(x, y)
    print ts, t1, ta, taa, tba
    #print (ts, ps)
    assert_equal(t1, ts)
    assert_equal(ts, tc)
    assert_equal(ps, pc)
for i in range(10):
    x = np.random.randn(20)
    y = np.random.randn(20)
    t1 = kendalltaustat(x, y)
    ts, ps = stats.kendalltau(x, y)
    ta = kendalltaustata(x, y)
    taa = gammaassoc(x, y)
    tba = kendalltaustatba(x, y)
    #print t1, (ts, ps)
    print ts, t1, ta, taa, tba
    assert_equal(t1,ts)


violence = np.array([1,1,1,2,2,2])
rating = np.array([1,2,3,1,2,3])
count = np.array([10, 5, 2, 9, 12, 16])
vi = np.repeat(violence,count)
ra = np.repeat(rating,count)

print 'stats.kendalltau, 0.34932342479397172',stats.kendalltau(vi, ra)
print 'stats.spearmanr, 0.370', stats.spearmanr(vi, ra)  # wrong
vir = stats.rankdata(vi)
rar = stats.rankdata(ra)
print 'pearson, 0.370', np.corrcoef(vi, ra)[0,1]
print 'rankdata correlation, 0.370', np.corrcoef(vir, rar)[0,1] # correct