[R-sig-Geo] Mann-kendall

[Forrest Stevens]

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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo