0

我正在尝试对两个光栅堆栈进行逐像素格兰杰因果检验，每个堆栈有 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

请问如何为栅格修改此功能？

4

1 回答 1

1

您需要检查grangertest返回的内容：

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]
          }
     }

光栅堆栈示例：

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)
于 2016-06-22T00:05:07.870 回答