[R-sig-Geo] Mann-kendall

[Forrest Stevens]

You're welcome!  I'm glad you found it useful. I agree that the
MannKendall package and functions could use slightly more informative
warnings, though many packages won't give you any such warnings at
all. At least M-K clues you in that you may be making a statistical
error rather than letting one run on completely uninformed. Good luck
getting it all running!

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 7:58 PM, Nuno Sá <nunocesardesa at gmail.com> wrote:
> Hello!
>
> ah! that "return(unlist(MannKendall(x)))" might have been why i didn't
> manage to do it with the calc function,  will try it out, thank you!
>
> I have a landsat image of a large study area, about 1 092 005 cells by 7
> "measurements", hence my avoidance of pixel by pixel but I am going to try
> yours. (I also know I shouldn't be using this test for significance of the
> MK here but this is just a test set that I am using for developing the
> algorithm). Thank you also for the suggestion on how to generate the data.
>
> I'm going to leave the kendall you sugested running through the night. I
> have the following warning though:
>
> WARNING: Error exit, tauk2, IFAULT =  10
>
> I've been reading and it might be related to the low number of n for my
> testing case.
> Found online in a response from the package maintainer A.I. McLeod:
>
> "Thank you though for your comments.  So I will improve the documentation
> for Kendall by terminating the program with an error message when n<=3
> (this case is of no interest to me) and warning message when n<12 that the
> p-values may be inaccurate."
>
>
> Thank you for everything!
> Nuno
>
>
> On 9 October 2014 00:26, Forrest Stevens <forrest at ufl.edu> wrote:
>>
>> 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
>
>
>
>
> --
>
> Nuno César de Sá
> +351 91 961 90 37
>