2

corと を使用して、2 つのラスター ブリック間のピクセル単位の相関と有意性 (p 値) を取得しようとしていますcor.test。私のデータはここにあります：

どちらもかなり小さく、全体で 2MB 未満です。

StackOverflow と r-sig-geo に関する以前の議論から、次の 2 つのコード (両方とも Robert Hijmans から) を見つけました。

#load data
require(raster)
sa <- brick("rain.grd")
sb <- brick("pet.grd")

# code 1
funcal <- function(xy) {
    xy <- na.omit(matrix(xy, ncol=2))
    if (ncol(xy) < 2) {
        NA
    } else {
        cor(xy[, 1], xy[, 2])
    }
}
s <- stack(sa, sb)
calc(s, funcal)
# code 2
stackcor <- function(s1, s2, method='spearman') { 
        mycor <- function(v) { 
                x <- v[1:split] 
                y <- v[(split+1):(2*split)] 
                cor(x, y, method=method) 
        } 
        s <- stack(s1, s2) 
        split <- nlayers(s)/2 
        calc(s, fun=mycor ) 
}

どちらのコードも、cor相関グリッドを生成する関数で期待どおりに機能します。ただし、 p 値を抽出するために を にcor置き換えてみました。cor.test

# using the first code
funcal <- function(xy) {
    xy <- na.omit(matrix(xy, ncol=2))
    if (ncol(xy) < 2) {
        NA
    } else {
        cor.test(xy[, 1], xy[, 2],method="pearson")$p.value
    }
}
s <- stack(sa, sb)
calc(s, funcal)

次のエラーが発生しました（RStudioのトラックバックを使用）：

Error in cor.test.default(xy[, 1], xy[, 2], method = "pearson") : 
  not enough finite observations 
8 stop("not enough finite observations") 
7 cor.test.default(xy[, 1], xy[, 2], method = "pearson") 
6 cor.test(xy[, 1], xy[, 2], method = "pearson") 
5 FUN(newX[, i], ...) 
4 apply(v, 1, fun) 
3 .local(x, fun, ...) 
2 calc(meistack, brick.cor.pval) 
1 calc(meistack, brick.cor.pval)

以前の r-sig-geo ディスカッションで、ユーザーがこのエラーについて質問しましたが、回答がありませんでした。corそのため、もう一度質問したところ、行列を入力できて入力できないと指摘されたがcor.test、データを数値ベクトルに変換した後でも、問い合わせに対する 1 つの応答を受け取りました。

cor.pval <- function(xy) { # Pearson product moment correlation
  xy <- na.omit(as.matrix(xy, ncol=2))
  x <- xy[,1:11]
  y <- xy[,12:22]
  #  if (ncol(xy) < 2) {
   # NA
  #} else {
    cor.test(x[1,], y[1,])$p.value
  #}
}
calc(s, cor.pval)

次のエラーに直面しています：

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

誰かがこれを手伝ってくれるかどうか疑問に思っていましたか？

私の sessionInfo() は以下の通りです:

R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] car_2.0-18          nnet_7.3-7          MASS_7.3-29         rgdal_0.8-11       
 [5] plyr_1.8            rasterVis_0.21      hexbin_1.26.2       latticeExtra_0.6-26
 [9] RColorBrewer_1.0-5  lattice_0.20-23     raster_2.1-49       maptools_0.8-27    
[13] sp_1.0-13          

loaded via a namespace (and not attached):
[1] foreign_0.8-55 tools_3.0.1    zoo_1.7-10

ありがとう！

4

2 に答える 2