Mann-Kendall

S and sign

When I first started this deep dive into Mann-Kendall I got a little hung up on the sign function aka signum. It’s a little abstract at first (IMO) but is really just a counting algorithm. At it’s core, it really is just figuring out whether a value is \(\lt 0\), \(\gt 0\) or \(= 0\). The sign function is within the \(S\) statistic. The \(S\) statistic is approximately normally distributed (remember above and solving for \(u_c\)) with a mean of 0 and a variance of \(\sqrt{Var(S)}\). Below is the equation for \(S\),

For example, let’s say we had some observations over time (important because they are ordered) and we wanted to know whether the point after (\(X(t')\)) is bigger or smaller than (\(X(t)\)) we could just take \(X(t')-X(t)\) and if it’s negative it means it’s smaller and if it’s positive it means it’s bigger (also 0 if they’re equal). But let’s say we do this as a rolling subtraction, e.g. \(x_1-x_1, x_2-(x_1, x_2),x_3-(x_1, x_2,x_3),\dots, x_n-(x_1,\dots,x_n)\). Now if we did this for each point as we progress to the end we would have a lot of positive and negative values and maybe some equal results right? In the code below you can see that this is what is going on, that is we are subtracting but rolling in time from the beginning to end.

x <- yaak_peaks$Peak
n <- length(x)
S <- 0.0   
for(j in 1:n) {
    S2 <- list(data.frame(rolling_sub = x[j] - x[1:j], name = paste0(j)))
    S <- append(S, S2)
  } 
head(S)

## [[1]]
## [1] 0
## 
## [[2]]
##   rolling_sub name
## 1           0    1
## 
## [[3]]
##   rolling_sub name
## 1       -4170    2
## 2           0    2
## 
## [[4]]
##   rolling_sub name
## 1       -5830    3
## 2       -1660    3
## 3           0    3
## 
## [[5]]
##   rolling_sub name
## 1       -5290    4
## 2       -1120    4
## 3         540    4
## 4           0    4
## 
## [[6]]
##   rolling_sub name
## 1       -4420    5
## 2        -250    5
## 3        1410    5
## 4         870    5
## 5           0    5

Let’s take a closer look in the graph below. The data.frame below S is just taking the peaks and subtracting \(X(t')-X(t)\) in this rolling pattern which inevitably is getting larger and larger samples as time goes on, e.g. \(x_n-(x_1,\dots,x_n)\). As you can see, values with more blue tend to be higher than values with more black. This will matter in the next steps.

Ok, the next step would be to add them up right? Well, sure but sign does this in an pretty cool and effective way that doesn’t let outliers drive the final results. Remember, the sign function is a counting method and all it does is determine whether the value is -, + or 0 and then gives either a -1, 1 or 0 as a result! By doing this the sign() makes it easy to interpret all of the -, + and 0’s.

\[ sign(a)=\begin{cases}\frac{a}{\left|a\right|}, & \text{if} \ x \neq 0 \\ 0, & \text{if} \ x = 0 \end{cases}=\large\begin{cases}1, & \text{if} \ a \gt 0 \\ 0, & \text{if} \ a = 0 \\ -1, & \text{if} \ a \lt 0 \end{cases} \]

So if we do this to the graph above we can see what each iteration expresses in terms of lower, higher or equal to in point graph. The point graph is a intuitive graphic for me in this problem because S is essentially a reflection of ‘the more +, - or 0’s you have as time goes on, the more this will determine your trend’.

sign_it <- S_df %>% 
  group_by(name) %>% 
  nest() %>%
  mutate(sign_result = map(data, ~sign(.$rolling_sub))) %>% 
  select(name, sign_result) %>% 
  unnest(sign_result) %>% 
  ungroup()

So the next step in solving for \(S\) is to sum up all of those points! This will give us a final sum and an reflection of a trend, e.g. more negative means downward and more positive means upward and closer to zero means no trend!

sum(sign_it$sign_result)

## [1] -373

Now if we look at the Kendall::MannKendall() function you’ll see we get the same answer. How cool is that!?

#Mann-Kendall function S score

Kendall::MannKendall(yaak_peaks$Peak)$S[1]

## [1] -373

\(\sqrt{Var(S)}\)

Now let’s move on to the variance part of the \(S\) equation. The \(VAR(S)\) will help us finish the \(S\) and also give us an idea of the what’s going on with the variation in the data. To solve for \(VAR(S)\) we need to use the equation below,

\[ VAR(S)=\frac{1}{18}\left[n(n-1)(2n+5)-\sum_{p-1}^{g}t_{p}(t_{p}-1)(2t_{p}+5)\right] \]

where \(n\) is equal to the sample size and \(t_p\) is equal to repeated values in the sample. Most of the time in hydrology we will not get repeated values when we sample but in this example we did so we’ll need to account for those repeats using \(\sum_{p-1}^{g}t_{p}(t_{p}-1)(2t_{p}+5)\) above, e.g. peaks that are the same. There are also some rules for sample size where \(n \le 10\) will need to use a table of probabilities and not perform \(VAR(S)\) but most of the time you’ll have over 10 samples (i hope 😉).

yaak_peaks %>% count(Peak) %>% filter(n>1)

## # A tibble: 5 x 2
##    Peak     n
##   <dbl> <int>
## 1  4600     2
## 2  6490     3
## 3  6600     2
## 4 10200     2
## 5 12100     2

So let’s solve for \(VAR(S)\)! We’ll need the length of the vector to get \(n\) and we’ll need to perform the summation of the repeated values.

n <- length(yaak_peaks$Peak)
repeats <- yaak_peaks %>% count(Peak) %>% filter(n>1)

repeat_res <- vector()
for(i in 1:length(repeats$n)){
  rr <- repeats$n[i]*(repeats$n[i]-1)*(2*repeats$n[i] + 5)
  repeat_res <- append(repeat_res, rr)
}

var_s <- (1/18)*(n*(n-1)*(2*n+5) - sum(repeat_res))
var_s

## [1] 31192.33

And now if we look at the Kendall::MannKendall() result we can see that they are the same!

Kendall::MannKendall(yaak_peaks$Peak)$varS[1]

## [1] 31192.33

With that we can now complete the equation for \(S\)!