2016-06-21 8 views
0

각각 60 개의 래스터가있는 두 래스터 스택에서 픽셀 단위의 Granger 인과 테스트를 수행하려고합니다.Raster Stack에 대한 Granger 인과 테스트 R

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"패키지 번째 사용하여 원래의 기능은 다음과 같습니다 : 아래의 예는 20 래스터를 가지고

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 

래스터를 위해이 기능을 어떻게 수정합니까?

답변

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 

그러면 우리는 뭔가를 li로 변경할 수 있습니다. KE :

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) 
+0

감사 @RobertH. 이것은 효과가 있었다. –