Mann-kendall

Not to dissuade any interesting new algorithm development, and I may
be missing something in the motivation for what you're trying to
accomplish, but have you looked at the Kendall package and tried using
it on your raster object?  I've been pleasantly surprised at its
efficiency for at least moderately sized raster stacks.  I've
implemented something like the following on real MODIS time series and
it works well (simulated data provided here):


#install.packages("Kendall")

require(raster)
require(Kendall)

## Simulate a small 25x25x25 raster with uncorrelated but trended series:
set.seed(2002)
slopes <- matrix(rnorm(625, 0, 1), nrow=25, ncol=25)

r <- raster(slopes*1+matrix(rnorm(625, 0, 1), nrow=25, ncol=25))
s <- stack(r)

for (i in 2:25) {
r <- raster(slopes*i+matrix(rnorm(625, 0, 1), nrow=25, ncol=25))
s <- addLayer(s, r)
}


## Confirm on a series with a high slope:
slopes[1,3]
plot(1:25, s[1,3,])
MannKendall(s[1,3,])

## vs. one that has a very small slope:
slopes[1,1]
plot(1:25, s[1,1,])
MannKendall(s[1,1,])


## Mann-Kendall on one pixel:
MannKendall(s[1,1,])


## Perform Mann-Kendall on all pixels:
a <- calc(s, fun=function(x) { return(unlist(MannKendall(x)))})

## Plot pixel-level Kendall Tau and p-values:
plot(a$tau)
plot(a$sl)



Again, I might be missing something about your goals but hopefully
this helps someone out there with similar problems.

Sincerely,
Forrest

--
Forrest R. Stevens
Ph.D. Candidate, QSE3 IGERT Fellow
Department of Geography
Land Use and Environmental Change Institute
University of Florida
www.clas.ufl.edu/users/forrest


On Wed, Oct 8, 2014 at 6:02 PM, Nuno S? <nunocesardesa at gmail.com> wrote:

Hello

Using the integration by part summation was quite fast (I owe a beer to

the

ones who've developed the raster package!!)

So, let variable Z.combi be the Z value of the Mann-kendall.

And PHI(x) be the CDF of the Standar normal distribution w/ mean 0 and

var

1.

According to Neeti the p-value is given by:

p = 2 * (1 - PHI(|Z|)

Now, on the code, using the integral expansion:

-----------------------

library(phangorn) #to have double factorials

phi.aprox <- 0.5+exp((-1/2)*abs(Z.combi)^2)/sqrt(2*pi)
phi.integ <- 0

for (n in 1:100){

  phi.integ <- phi.integ + Z.combi^(2*n+1)/dfactorial(2*n+1)
}

phi.aprox <- phi-aprox*phi.integ

p.value.rst.neeti <- 2*(1-phi.aprox)

-------------------

The result falls within 0 to 1 which is a good indication. But all

depends

on everything I've asked before.

Thank you!



On 8 October 2014 22:23, Nuno S? <nunocesardesa at gmail.com> wrote:

Hello!

I've been trying for two days to implement the Mann-Kendall test. It is
available in IDRISI and it would be cool that it would be available

here.

I believe I made part of it correctly but I am stuck in two steps.

1. Calculating the ties/NA coefficients. (Algorithm problem..)

2. Estimating the p-value on the last step (Lack of statistical

knowledge

problem maybe)

So here it goes what I have so far. Bear with me that I have tried to

keep

my variables names with a certain "instant" meaning.

---------------------- Problem here, is identifying the coefficients for
the ties/NA -------------------

dummy.stack <- raster(stack.ndvi,layer=1) #required to initiate stack
within loop
mtx <- matrix(0,ncol=nlayers(stack.ndvi),nrow=(nlayers(stack.ndvi)-1))
#this allows the construction of a matrix with tie values signalled to 1

for (i in 1:(nlayers(stack.ndvi)-1)){

  print("----- i ----")
  print(i)
  for (j in (i+1):nlayers(stack.ndvi)){

    print(j)

    dummy.stack <-

stack(dummy.stack,sign(raster(stack.ndvi,layer=i)-raster(stack.ndvi,layer=j)))

    tie.break <-
sign(raster(stack.ndvi,layer=i)-raster(stack.ndvi,layer=j)) == 0
    cc <- maxValue(tie.break)
    if( cc == 1){
      mtx[i,j-1]=1
    }
  }
}

S.stack <- dropLayer(dummy.stack,1) #removing dummy

S.stati <- raster(S.stack,layer=1) #for the somatory - initiating in the
second term (can be substituted by in-built overlay function)

for (i in 2:nlayers(S.stack)){

  S.stati <- S.stati + raster(S.stack,layer=i)
}

-----------------------------------------------------------------------

#Ignoring the existence of ties/na (which I've read is acceptable for
large enough datasets)
#then:

  var.mk.sqr <- sqrt(n*(n-1)*(2*n+5)/18)

  S.c1 <- ((S.stati > 0)*S.stati - 1 )/ var.mk.sqr
  S.c2 <- (S.stati ==0)*S.stati* 0
  S.c3 <- ((S.stati < 0)*S.stati + 1 )/ var.mk.sqr

  Z.combi <- overlay(S.c1,S.c2,S.c3,fun="sum")

#"The P value (probability value p) of the MK statistic S of sample data
can
be estimated using the normal cumulative distribution function:"

(Yue,2002

and Neeti, 2011)
# So, how to implement the CDF using raster data? Am I forced to do

series

expansions? would this not be to slow with raster? Or is there a quicker
way to iterate pixel by pixel?

Sorry for the long email and thank you in advance for any
comments/corrections/recommendations!

Best regards to all,
Nuno S?


--

Nuno C?sar de S?
+351 91 961 90 37

--

Nuno C?sar de S?
+351 91 961 90 37

        [[alternative HTML version deleted]]