rlmr-rastercausality

Granger causal tests on raster stacks in R

I am attempting to do a pixel-wise granger causal test on two raster stacks with 60 rasters each. The example below has only 20 rasters:

library(raster)
library(lmtest)

r <- raster(ncol=10, nrow=10)
r[]=1:ncell(r)
S <- stack(r,r,r,r,r,r,r,r,r,r,r,r)
R <- stack(r,r,r,r,r,r,r,r,r,r,r,r)
FNO2<-stack(S,R)

The original function using th "lmtest" package is:

D<- grangertest(degg ~ dchick, order=4)

Here's a modification I did to run the original grangertest function on raster stacks?

funG <- function(x) { if (is.na(x[1])){ NA } else { grangertest(x[13:24] ~ x[1:12],order=1 )}}
granger<-calc(FNO2,funG)

Where FNO2 is the stack of both raster stacks. I get the error below:

Error in `colnames<-`(`*tmp*`, value = c("x", "y", "x_1", "y_1")) : 
length of 'dimnames' [2] not equal to array extent

How do I modify this function for rasters please?

Solution

  • You need to inspect what grangertest returns:

    library(lmtest)
x <- grangertest(egg ~ chicken, order = 3, data = ChickEgg)

class(x)
#[1] "anova"      "data.frame"

x
#Granger causality test
#
#Model 1: egg ~ Lags(egg, 1:3) + Lags(chicken, 1:3)
#Model 2: egg ~ Lags(egg, 1:3)
#  Res.Df Df      F Pr(>F)
#1     44                 
#2     47 -3 0.5916 0.6238

    That is not something we can stick in a RasterLayer

    str(x)
#Classes ‘anova’ and 'data.frame':       2 obs. of  4 variables:
# $ Res.Df: num  44 47
# $ Df    : num  NA -3
# $ F     : num  NA 0.592
# $ Pr(>F): num  NA 0.624
# - attr(*, "heading")= chr  "Granger causality test\n" "Model 1: egg ~ Lags(egg, 1:3) + Lags(chicken, 1:3)\nModel 2: egg ~ Lags(egg, 1:3)"
>

    I am not sure what you are after, but if it were the p value, perhaps

    x$'Pr(>F)'[2]
#[1] 0.6237862

    Then we can do change the function to something like:

    funG <- function(x) { if (is.na(x[1])){ 
               NA 
          } else { 
               grangertest(x[13:24] ~ x[1:12],order=1 )$'Pr(>F)'[2]
          }
     }

    Example RasterStack:

    library(raster)
r <- raster(ncol=10, nrow=10)
set.seed(9)
FNO2 <- stack(sapply(1:24, function(i) setValues(r, runif(ncell(r)))))

    Test the function:

    d <- as.vector(FNO2[1])
funG(d)
#[1] 0.03248222

    Now:

    granger<-calc(FNO2,funG)

granger
#class       : RasterLayer 
#dimensions  : 10, 10, 100  (nrow, ncol, ncell)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#data source : in memory
#names       : layer 
#values      : 0.007425836, 0.9895796  (min, max)