R 中栅格堆栈的 Granger 因果检验

Granger causal tests on raster stacks in R

我正在尝试对两个栅格堆栈进行像素级格兰杰因果检验,每个栅格堆栈有 60 个栅格。下面的示例只有 20 个栅格:

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)

使用 "lmtest" 包的原始函数是:

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

这是我对 运行 栅格堆栈上的原始 grangertest 函数所做的修改?

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

其中 FNO2 是两个光栅堆栈的堆栈。我收到以下错误:

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

请问如何为栅格修改此函数?

你需要检查什么 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

那不是我们可以在 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)"
> 

我不确定你在追求什么,但如果是 p 值,也许

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

然后我们可以将函数更改为:

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

RasterStack 示例:

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

测试功能:

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

现在:

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)