R: multispatialCCMをためしてみる：Taglibro de H：SSブログ

	ブログをはじめるログイン

R: multispatialCCMをためしてみる　[統計] [編集]

多地点でとった時系列データから因果関係の推定をおこなうRパッケージ"multispatialCCM"をためしてみました。

文献

Multispatial Convergent Cross Mappingの論文
Clark AT, Ye H, Isbell F, et al. (2015) Spatial convergent cross mapping to detect causal relationships from short time series. Ecology 96:1174–1181. doi: 10.1890/14-1479.1
その元になったConvergent Cross Mappingの論文
Sugihara G, May R, Ye H, et al. (2012) Detecting Causality in Complex Ecosystems. Science (80- ) 338:496–500. doi: 10.1126/science.1227079

まずはマニュアルにある例題をなぞってみます。

library(multispatialCCM)

## テスト用データ生成
ccm_data_out <- make_ccm_data()

データを表示してみます。

> ccm_data_out
$Accm
  [1]         NA 0.43901652 0.79590473 0.45745491 0.83460391 0.45385562 0.84711147 0.46431136 0.85405916
 [10] 0.44028015 0.80988939         NA 0.77394628 0.50868099 0.75072820 0.54043692 0.71517458 0.64343894
 [19] 0.67257285 0.66059170 0.59307618 0.76662347         NA 0.81329423 0.24778903 0.65303490 0.65149461
 [28] 0.77784334 0.42598562 0.82328657 0.17839953 0.55973942 0.79396914         NA 0.61901154 0.76976034
 [37] 0.39202567 0.72334355 0.40143345 0.77342292 0.32859612 0.73835989 0.45910903 0.70411524         NA
 [46] 0.07408985 0.25883116 0.71985429 0.65495765 0.71721333 0.59018707 0.75129321 0.53257335 0.78749496
 [55] 0.58023205         NA 0.52757398 0.90314512 0.16201049 0.45751593 0.89654287 0.14312029 0.42308294
 [64] 0.91370467 0.22389105 0.56685459         NA 0.33905268 0.73830754 0.41803639 0.76309231 0.37411698
 [73] 0.77916190 0.32660413 0.73105269 0.38595567 0.83359802         NA 0.82682676 0.40456035 0.81540750
 [82] 0.41965262 0.81463214 0.39700595 0.83158657 0.43483506 0.86150160 0.39205463         NA 0.80021704
 [91] 0.40485886 0.81042904 0.46213196 0.86347764 0.35071113 0.81890450 0.54824255 0.83420078 0.41293884
[100]         NA 0.50315667 0.82310454 0.38496252 0.80126249 0.52895966 0.80645166 0.48950919 0.77597976
[109] 0.48378555 0.81250842         NA 0.77108294 0.51237159 0.76134275 0.57784527 0.79057760 0.46140948
[118] 0.79948424 0.33500715 0.77095052 0.49105923         NA 0.82556749 0.48095290 0.80857228 0.46419525
[127] 0.82178492 0.47398297 0.79040277 0.47869930 0.79937962 0.48158445         NA 0.67433002 0.58197186
[136] 0.72390650 0.51124741 0.81188930 0.24630789 0.62026584 0.66830317 0.60177587 0.69187134         NA
[145] 0.37106585 0.82783432 0.49098194 0.82730979 0.40712843 0.82023510 0.47244781 0.84475828 0.38969284
[154] 0.80092820         NA 0.73744867 0.47758898 0.83066887 0.29786583 0.74315775 0.52140186 0.83230047
[163] 0.27490998 0.69124363 0.53008126         NA 0.69420245 0.65841239 0.59679112 0.80473088 0.26779278
[172] 0.69992847 0.56164770 0.83440522 0.22971113 0.61302620         NA 0.37020055 0.75735869 0.40216174
[181] 0.82121564 0.29879303 0.67701271 0.58641414 0.68690893 0.54920165 0.76971380         NA 0.73939765
[190] 0.62742746 0.71950126 0.67072640 0.61824435 0.79385562 0.31657031 0.77840676 0.38995558 0.82817371
[199]         NA 0.80341495 0.43547307 0.80058175 0.53485761 0.82555841 0.41682826 0.75587140 0.51089785
[208] 0.79190020 0.48775917         NA 0.76485978 0.52061322 0.78038398 0.50492346 0.77698532 0.53808179
[217] 0.78969978 0.50552758 0.79326317 0.52442889

$Bccm
  [1]        NA 0.8446986 0.4163134 0.8012023 0.4628400 0.8666743 0.4201954 0.8401107 0.4743030 0.8638343
 [11] 0.4250372        NA 0.6395852 0.8506177 0.4704501 0.9225798 0.3048221 0.7453019 0.6493348 0.8405217
 [21] 0.4603383 0.9407567        NA 0.7895522 0.5499925 0.8645317 0.3693855 0.8760876 0.4410223 0.8865948
 [31] 0.3154533 0.7923424 0.6086976        NA 0.2844921 0.8004781 0.6572119 0.8386878 0.4835386 0.9746074
 [41] 0.1883518 0.5610565 0.9273107 0.2513766        NA 0.2772417 0.7131462 0.7033820 0.7564373 0.6534527
 [51] 0.8528067 0.5211974 0.9448948 0.2670747 0.7019178        NA 0.1724594 0.5327164 0.9817420 0.1488442
 [61] 0.4785245 0.9895916 0.1338863 0.4565686 0.9496481 0.1701330        NA 0.2927145 0.7713924 0.6407255
 [71] 0.8854483 0.3396591 0.8798452 0.4131358 0.9329670 0.2030786 0.6421861        NA 0.4577212 0.8555631
 [81] 0.4439446 0.8295234 0.4369317 0.8487085 0.4505406 0.8218474 0.4748893 0.8877870        NA 0.5750282
 [91] 0.8730500 0.3627040 0.8176051 0.5060148 0.9398667 0.3413530 0.7841011 0.5732435 0.9040860        NA
[101] 0.7946518 0.5805083 0.8117651 0.4197869 0.8842215 0.4003038 0.8628444 0.4109841 0.8984029 0.4060367
[111]        NA 0.4571004 0.9530335 0.2699988 0.7499394 0.7316465 0.7551385 0.7632637 0.7740296 0.7036927
[121] 0.7874270        NA 0.4621049 0.8404663 0.4728596 0.8164959 0.4641774 0.8282664 0.4857008 0.8258110
[131] 0.4585424 0.8281014        NA 0.7784526 0.6729299 0.8468216 0.4958559 0.9176932 0.2618412 0.6436158
[141] 0.8161512 0.6084921 0.8690876        NA 0.8794025 0.3844722 0.8630498 0.5102798 0.8466225 0.3586957
[151] 0.8469254 0.5003982 0.8943609 0.3722045        NA 0.8187618 0.5382421 0.8605279 0.3858501 0.8348715
[161] 0.4890104 0.8601642 0.3451875 0.8198689 0.5597063        NA 0.4436194 0.8965343 0.3272831 0.7917398
[171] 0.6072002 0.8710965 0.3708173 0.8347499 0.4110317 0.9051043        NA 0.4965309 0.8951191 0.2703243
[181] 0.7592161 0.7515354 0.7431187 0.7153286 0.7946363 0.6661621 0.8011591        NA 0.4175404 0.8671194
[191] 0.4319680 0.9074746 0.3815152 0.8207540 0.4502697 0.8857891 0.3624559 0.8505024        NA 0.5309351
[201] 0.9428647 0.3199439 0.7978699 0.5593886 0.8892320 0.3759853 0.8605281 0.4745956 0.8922223        NA
[211] 0.3862380 0.7963640 0.4993282 0.8676377 0.3888409 0.8591025 0.4622835 0.8576267 0.3921987 0.8328497

$time_ccm
  [1] NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1
 [36]  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3
 [71]  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5
[106]  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7
[141]  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9
[176] 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA  1  2  3  4  5  6  7  8  9 10 NA
[211]  1  2  3  4  5  6  7  8  9 10

NAが、時系列データの地点間の区切りになっています。

データをグラフでみてみます。

plot(Accm ~ time_ccm, data = ccm_data_out, type = "l", col = 2,
     ylim = c(min(c(ccm_data_out$Accm, ccm_data_out$Bccm), na.rm = TRUE),
              max(c(ccm_data_out$Accm, ccm_data_out$Bccm), na.rm = TRUE)),
     las = 1)
lines(Bccm ~ time_ccm, data = ccm_data_out, col = 4)

AccmとBccmとをプロットしてみます。1期ずらしたものも表示してみます。

plot(Bccm ~ Accm, data = ccm_data_out, type = "p",
     xlim = c(0, 1), ylim = c(0, 1), las = 1, asp = 1)
plot(Bccm[-1] ~ Accm[-length(Accm)], data = ccm_data_out, type = "p",
     xlim = c(0, 1), ylim = c(0, 1), las = 1, asp = 1,
     xlab = "Accm[t]", ylab = "Bccm[t+1]")
plot(Bccm[-length(Bccm)] ~ Accm[-1], data = ccm_data_out, type = "p",
     xlim = c(0, 1), ylim = c(0, 1), las = 1, asp = 1,
     xlab = "Accm[t+1]", ylab = "Bccm[t]")

これだけみても傾向はあまりはっきりとはしないようです。

ここから解析です。まず、データをデータフレームからとりだしておきます。

Accm <- ccm_data_out$Accm
Bccm <- ccm_data_out$Bccm

最適なembedding dimensionを決定します。論文によると、embedding dimensionとは、LOO CVで予測能力がもっともたかくなるようなtime step数とのことです。

# 最適なEを求める
maxE <- 5 # Maximum E to test
# 出力を格納する行列
Emat <- matrix(nrow = maxE - 1, ncol=2)
colnames(Emat) <- c("A", "B")

for(E in 2:maxE) {
  #Uses defaults of looking forward one prediction step (predstep)
  #And using time lag intervals of one time step (tau)
  Emat[E - 1, "A"] <- SSR_pred_boot(A = Accm, E = E,
                                    predstep = 1, tau = 1)$rho
  Emat[E - 1, "B"] <- SSR_pred_boot(A = Bccm, E = E,
                                    predstep = 1, tau = 1)$rho
}

# グラフをみて、rhoを最大化するEを決める
matplot(2:maxE, Emat, type = "l", col = 1:2, lty = 1:2,
          xlab = "E", ylab = "rho", lwd = 2, las = 1)
legend("bottomleft", c("A", "B"),
       lty = 1:2, col = 1:2, lwd = 2, bty = "n")

ということで、Aについては2、Bについては3ということになりました。

E_A <- 2
E_B <- 3

次に、非線形性と確率的ノイズのチェックをします。論文によると、CCMはcoupled nonlinear systemに適用されるのでシステムが完全にランダムなものではないこと、またノイズが大きすぎないことを確認する必要があるとのことです。

signal_A_out <- SSR_check_signal(A = Accm, E = E_A, tau = 1,
                                 predsteplist = 1:10)
signal_B_out <- SSR_check_signal(A = Bccm, E = E_B, tau = 1,
                                 predsteplist = 1:10)

結果です。

> signal_A_out
$predatout
   predstep        rho
1         1 0.41092654
2         2 0.37104873
3         3 0.23067061
4         4 0.18685019
5         5 0.03532828
6         6 0.09244104
7         7 0.18237867
8         8 0.12596805
9         9 0.02908599
10       10 0.49743669

$rho_pre_slope
   Estimate    Pr(>|t|) 
-0.01269674  0.50750092 

$rho_predmaxCI
[1] 0.3849831 0.5953362

> signal_B_out
$predatout
   predstep         rho
1         1  0.55893684
2         2  0.62654258
3         3  0.34631536
4         4  0.26467229
5         5  0.23766947
6         6  0.16090055
7         7  0.13669368
8         8  0.08082902
9         9  0.52535764
10       10 -0.02419999

$rho_pre_slope
   Estimate    Pr(>|t|) 
-0.04693735  0.04188685 

$rho_predmaxCI
[1] -0.1690604  0.1216839

非線形なシステムであれば、predstepの増加に対してrhoが減少するはず（ただし逆はなりたたない）とのことです。一方、完全にランダムなシステムではrhoは一貫してひくい値となるということです。

いよいよ解析実行です。100のブートストラップをおこなっています。

# AがBの原因か?
CCM_boot_A <- CCM_boot(Accm, Bccm, E_A, tau = 1, iterations = 100)
# BがAの原因か?
CCM_boot_B <- CCM_boot(Bccm, Accm, E_B, tau = 1, iterations = 100)

#Test for significant causal signal
#See R function for details
CCM_significance_test <- ccmtest(CCM_boot_A, CCM_boot_B)

> CCM_significance_test
pval_a_cause_b pval_b_cause_a 
           0.2            0.0

BがAの原因になっているという仮説は有意でしたが、AがBの原因になっているという仮説は有意ではありませんでした。

結果をグラフで表示します。

plotxlimits <- range(c(CCM_boot_A$Lobs, CCM_boot_B$Lobs))

#Plot "A causes B"
plot(CCM_boot_A$Lobs, CCM_boot_A$rho, type = "l", col = 1, lwd = 2,
     xlim = c(plotxlimits[1], plotxlimits[2]), ylim = c(0, 1),
     xlab = "L", ylab = "rho")
#Add +/- 1 standard error
matlines(CCM_boot_A$Lobs,
         cbind(CCM_boot_A$rho - CCM_boot_A$sdevrho,
               CCM_boot_A$rho + CCM_boot_A$sdevrho),
         lty = 3, col = 1)

#Plot "B causes A"
lines(CCM_boot_B$Lobs, CCM_boot_B$rho, type = "l", col = 2, lty = 2, lwd = 2)
#Add +/- 1 standard error
matlines(CCM_boot_B$Lobs,
         cbind(CCM_boot_B$rho - CCM_boot_B$sdevrho,
               CCM_boot_B$rho + CCM_boot_B$sdevrho),
         lty = 3, col = 2)

legend("topleft",
       c("A causes B", "B causes A"),
       lty = c(1, 2), col = c(1, 2), lwd = 2, bty = "n")

BがAの原因になっているというモデルの方がrhoの値がたかくなっています。