0
$\begingroup$

I'm trying to perform a Mann-Kendall Test of some data. I'm using the code in the following link (https://github.com/mps9506/Mann-Kendall-Trend/blob/master/mk_test.py) modified a little bit bellow for the result to be in an array and it only gives back a p-value (p) and a tau value (z).

def mk_test(x, alpha=0.05):
    """
    This function is derived from code originally posted by Sat Kumar Tomer
    (satkumartomer@gmail.com)
    See also: http://vsp.pnnl.gov/help/Vsample/Design_Trend_Mann_Kendall.htm
    The purpose of the Mann-Kendall (MK) test (Mann 1945, Kendall 1975, Gilbert
    1987) is to statistically assess if there is a monotonic upward or downward
    trend of the variable of interest over time. A monotonic upward (downward)
    trend means that the variable consistently increases (decreases) through
    time, but the trend may or may not be linear. The MK test can be used in
    place of a parametric linear regression analysis, which can be used to test
    if the slope of the estimated linear regression line is different from
    zero. The regression analysis requires that the residuals from the fitted
    regression line be normally distributed; an assumption not required by the
    MK test, that is, the MK test is a non-parametric (distribution-free) test.
    Hirsch, Slack and Smith (1982, page 107) indicate that the MK test is best
    viewed as an exploratory analysis and is most appropriately used to
    identify stations where changes are significant or of large magnitude and
    to quantify these findings.
    Input:
        x:   a vector of data
        alpha: significance level (0.05 default)
    Output:
        trend: tells the trend (increasing, decreasing or no trend)
        h: True (if trend is present) or False (if trend is absence)
        p: p value of the significance test
        z: normalized test statistics`

    Examples
    --------
      >>> x = np.random.rand(100)
      >>> trend,h,p,z = mk_test(x,0.05)
    """
    n = len(x)

    # calculate S
    s = 0
    for k in range(n-1):
        for j in range(k+1, n):
            s += np.sign(x[j] - x[k])

    # calculate the unique data
    unique_x = np.unique(x)
    g = len(unique_x)

    # calculate the var(s)
    if n == g:  # there is no tie
        var_s = (n*(n-1)*(2*n+5))/18
    else:  # there are some ties in data
        tp = np.zeros(unique_x.shape)
        for i in range(len(unique_x)):
            tp[i] = sum(x == unique_x[i])
        var_s = (n*(n-1)*(2*n+5) - np.sum(tp*(tp-1)*(2*tp+5)))/18

    if s > 0:
        z = (s - 1)/np.sqrt(var_s)
        #result = (s - 1)/np.sqrt(var_s)
    elif s == 0:
         z = 0
        #result = 0
    elif s < 0:
        z = (s + 1)/np.sqrt(var_s)
        #result = (s + 1)/np.sqrt(var_s)

    # calculate the p_value
    p = 2*(1-norm.cdf(abs(z)))  # two tail test
    result= np.append(p,z)
    h = abs(z) > norm.ppf(1-alpha/2)

    return np.array(result)

I am then using the following code to perform the test.

out = np.empty((0))
for i in range(145):
    for j in range(192):
        out1 = mk_test(yrmax[:,i,j], alpha=0.05)
        out = np.append(out, out1, axis=0)

I think something is going wrong somewhere while performing the test because I would expect to get z values between -1 and 1, but I'm getting back some values that are greater than 1. Is the coding going wrong or am I misunderstanding what z is and it's not actually tau and hence why I'm getting values I'm not expecting?

$\endgroup$
  • $\begingroup$ $z$ is not $\tau$. $\endgroup$ – Glen_b -Reinstate Monica Nov 1 '17 at 7:15
  • $\begingroup$ Attn close voters: Reopened as the underlying issue is plainly statistical - the answer is an explanation of the difference between $\tau$ and $z(\tau)$, for goodness sake. $\endgroup$ – Glen_b -Reinstate Monica Nov 1 '17 at 22:40
3
$\begingroup$

The $\tau$ statistic is (in effect) the number of pairs of points showing an increase minus the number of pairs of points showing a decrease all divided by the number of pairs of points showing an increase or decrease (which if there's no ties, is just the number of pairs of points).

Incorporating adjustment for ties in the series (but assuming one value per time instant, so no ties there), $\tau=S/D$, where $D = \sqrt{[{n\choose 2}-\sum_j {t_j\choose 2}\}]{n\choose 2}}$ where $t_j$ is the number of ties at the $j$-th value which is non-unique (and if there are no ties, $D={n\choose 2}$ as mentioned above)

So (out of all pairs of points) it's the proportion of increasing pairs minus the proportion of decreasing pairs.

It's a number between -1 and 1.

A z score would divide by the standard deviation of the statistic (say the difference in counts divide by the standard deviation of the difference in counts or $\tau$ divided by its standard deviation, given the population Kendall correlation is 0), and as a result it can be much larger than 1 in absolute value -- in large samples the upper bound approaches $1.5\sqrt{n}$.

The larger $n$ is the smaller the standard deviation in $\tau$ under the null, so a given $\tau$ may correspond to a small-sized $Z$ at a small sample size and a large-sized $Z$ at a large sample size. That is, the relationship between $Z$ and $\tau$ is changes with sample size.

If the null hypotheses is true the $Z$ value will have mean 0 and standard deviation $1$ and in large samples will be approximately standard normal. In those circumstances it will nearly always be between say -2.5 and 2.5 or so, but when the null is false it may be quite far from 0.

With the Mann-Kendall test it seems to be standard to incorporate the continuity correction directly into the z-score (instead of using the correction when calculating the p-value); this just moves the difference in counts one* closer toward 0 before standardizing.

* why 1, rather than 1/2? Because the counts two sets of counts are perfectly negatively correlated so the difference in counts moves by two (one goes up the other goes down); the continuity correction is then half the smallest change.

$\endgroup$
  • $\begingroup$ Thank you for the clarification! I've never taken a stats class, so I'm still learning some stuff. Do you have a recommendation coding wise what to do to get tau? $\endgroup$ – Alex Morrison Nov 1 '17 at 16:28
  • $\begingroup$ Okay, let's assume a maximum of one data point at any given instant of time. See in the code where it has calculated S and var_s and then calculates z? In there where it's computing var_s, compute $D$ as in my edited answer -- where my $t_j$ is tp in the code. Then $\tau=S/D$. You should get the same value from anything that calculates Kendall's tau though, and I bet there are already a bunch of ways to just call something for that. $\endgroup$ – Glen_b -Reinstate Monica Nov 1 '17 at 22:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.