4  エルニーニョ南方振動

エルニーニョ・南方振動（ENSO: El Niño–Southern Oscillation）は熱帯太平洋の大気海洋結合現象で、その影響は熱帯だけでなく、大気の波動を通じて世界各地に及ぶ。

ENSOを監視するために、次の領域で平均された海面水温が指標して用いられている。

Code 
library(terra)

shift.lon <- function(v, dlon=0) {
  e <- crop(v, ext(dlon, 180, -90, 90))
  w <- crop(v, ext(-180, dlon, -90, 90))
  w <- shift(w, 360)
  rbind(e, w)
}

l50 <- shift.lon(vect("ne_50m_land/ne_50m_land.shp"))
ocean <- vect(ext(120, 300, -20, 20), crs(l50))
land <- crop(l50, ocean)
nino12 <- vect(ext(270, 280, -10, 0), crs(l50))
nino3 <- vect(ext(210, 270, -5, 5), crs(l50))
nino34 <- vect(ext(190, 240, -5, 5), crs(l50))
nino4 <- vect(ext(160, 210, -5, 5), crs(l50))
ninow <- vect(ext(130, 150, 0, 15), crs(l50))
par(mar = c(5, 4, 2, 1))
plot(ocean, col = "lightblue", border = NA,
     main = "ENSO monitoring regions",
     pax = list(
       xat = seq(120, 300, by = 30),
       yat = seq(-20, 20, by =10),
       cex.axis = 1.2
     ),
     xlab = "longitude", ylab = "",
     cex.lab = 1.2)
plot(land, border = "brown", col = "bisque", add = TRUE)
plot(nino12, border = "black", add = TRUE)
text(275, -5, "1+2")
plot(nino3, border = "blue", col = rgb(0, 0, 1, 0.3), add = TRUE)
text(250, 0, "Niño 3")
plot(nino4, border = "red", col = rgb(1, 0, 0, 0.3), add = TRUE)
text(175, 0, "Niño 4")
plot(nino34, border = "purple", lw = 2, add = TRUE)
text(220, 0, "Niño 3.4")
plot(ninow, border = "green", col = rgb(0, 1, 0, 0.3), lw = 2, add = TRUE)
text(140, 7.5, "west")

ENSO監視指数は、気象庁から入手できる。

このページのデータは、空白区切りで提供されている。 Rで空白区切りのデータを読むには、read.table()を使う。 Niño 3指数を読んでみよう。

df3 <- read.table("https://www.data.jma.go.jp/tcc/tcc/products/elnino/index/sstindex/base_period_9120/Nino_3/sst")
head(df3)
      JAN  FEB  MAR  APR  MAY  JUN  JUL  AUG  SEP  OCT  NOV  DEC
1949 25.1 26.5 26.7 27.3 27.1 26.1 25.0 24.7 24.2 24.0 24.1 24.1
1950 24.5 25.3 26.3 26.4 26.2 25.7 24.7 24.5 23.8 23.7 23.8 24.0
1951 24.5 25.6 26.4 27.1 27.1 26.8 26.5 26.1 25.8 26.3 26.3 25.9
1952 25.9 26.4 26.9 27.3 26.6 25.5 24.8 24.2 24.2 24.1 24.4 24.6
1953 25.6 26.6 27.3 27.8 27.5 27.0 26.0 25.4 25.5 25.2 25.1 25.0
1954 24.8 25.8 26.5 26.4 25.7 25.1 24.2 23.9 23.8 23.6 23.8 23.8
tail(df3)
      JAN  FEB  MAR  APR  MAY  JUN  JUL  AUG  SEP  OCT  NOV  DEC
2021 25.1 26.0 26.7 26.8 26.7 26.6 25.6 24.8 24.3 24.2 24.2 23.9
2022 24.3 25.3 26.5 26.6 26.0 25.9 25.4 24.6 23.9 23.9 24.1 24.3
2023 25.1 26.3 27.7 28.1 28.2 27.9 27.6 27.3 27.1 27.2 27.3 27.4
2024 27.7 28.1 28.4 28.3 27.1 26.4 25.7 24.7 24.7 24.6 24.9 24.7
2025 25.5 26.6 27.9 27.7 26.9 26.4 25.8 24.9 24.4 24.5 24.5 24.5
2026 25.2 26.5 27.5 28.3 99.9 99.9 99.9 99.9 99.9 99.9 99.9 99.9

各行が1949年からの年、各列が月になっている。 データのない月には 99.9 が入っているので、最後の年を削除する。

df3 <- df3[-nrow(df3), ]

data.frame は列の長さ（行の数、標本数）が同じであるリストなので、 c()を使って、一続きのベクトルにする。 各年毎に月の方向（行方向）に並べたベクトルを作りたいが、 Rはデータフレームでも列優先なので、各月毎に年の方向（列方向）にデータが格納されている。 t() で転置してから、 c() で次元の属性を取り、ベクトルにする。

nino3 <- c(t(df3))

平均と分散を計算してみよう。

\(X_i,\,i=1, \dots,n\) の平均 \(\mu\) は

\[ \mu = \frac{1}{n}\sum_{i=1}^n X_i \]

不偏分散は

\[ \sigma^2 = \frac{1}{n- 1}\sum_{i=1}^n (X_i - \mu)^2 \]

で定義される。

n <- length(nino3)
nino3_mean <- sum(nino3) / n
nino3_var <- sum((nino3 - nino3_mean)^2) / (n - 1)
print(nino3_mean)
[1] 25.77186
print(nino3_var)
[1] 1.639901

平均は mean() 、分散は var() で計算できる。

mean(nino3)
[1] 25.77186
var(nino3)
[1] 1.639901

各月の平均を apply() を使って計算する。 引数は、作用させる変数、次元、函数の順に与える。 行（次元1）が年、列（次元2）が月なので、各月の平均は次のように計算する。

apply(df3, 2, mean)
     JAN      FEB      MAR      APR      MAY      JUN      JUL      AUG 
25.45065 26.23377 27.00130 27.36883 26.93506 26.32857 25.56104 24.92597 
     SEP      OCT      NOV      DEC 
24.73896 24.80519 24.90779 25.00519

colMeans() を使うと効率的に計算できる。

colMeans(df3)
     JAN      FEB      MAR      APR      MAY      JUN      JUL      AUG 
25.45065 26.23377 27.00130 27.36883 26.93506 26.32857 25.56104 24.92597 
     SEP      OCT      NOV      DEC 
24.73896 24.80519 24.90779 25.00519

同様に各年の年平均も計算できる。

head(apply(df3, 1, mean))
    1949     1950     1951     1952     1953     1954 
25.40833 24.90833 26.20000 25.40833 26.16667 24.78333 
tail(apply(df3, 1, mean))
    2020     2021     2022     2023     2024     2025 
25.49167 25.40833 25.06667 27.26667 26.27500 25.80000 
head(rowMeans(df3))
    1949     1950     1951     1952     1953     1954 
25.40833 24.90833 26.20000 25.40833 26.16667 24.78333 
tail(rowMeans(df3))
    2020     2021     2022     2023     2024     2025 
25.49167 25.40833 25.06667 27.26667 26.27500 25.80000

月毎の不偏分散は apply() を使うと次のように計算できる。

apply(df3, 2, var)
      JAN       FEB       MAR       APR       MAY       JUN       JUL       AUG 
1.1256904 0.6454238 0.4101299 0.3545420 0.4862543 0.5996992 0.6953042 0.7706323 
      SEP       OCT       NOV       DEC 
0.8468831 1.0220779 1.2904648 1.3986569

不偏分散は平均を使って次のように書ける。

\[ \sigma^2 = \frac{1}{n-1}\left\{\sum_{i = 1}^{n}X_i^2 - \frac{1}{n}\left(\sum_{i = 1}^{n}X_i\right)^2\right\} = \frac{n}{n-1}\left\{\frac{1}{n}\sum_{i = 1}^{n}X_i^2 - \mu^2\right\} \]

この式を利用すると、高速な colSums() や colMeans() を使って分散を計算できる。

ny <- nrow(df3)
(colSums(df3^2) - colSums(df3)^2 / ny) / (ny - 1)
      JAN       FEB       MAR       APR       MAY       JUN       JUL       AUG 
1.1256904 0.6454238 0.4101299 0.3545420 0.4862543 0.5996992 0.6953042 0.7706323 
      SEP       OCT       NOV       DEC 
0.8468831 1.0220779 1.2904648 1.3986569 
(colMeans(df3^2) - colMeans(df3)^2) * ny / (ny - 1)
      JAN       FEB       MAR       APR       MAY       JUN       JUL       AUG 
1.1256904 0.6454238 0.4101299 0.3545420 0.4862543 0.5996992 0.6953042 0.7706323 
      SEP       OCT       NOV       DEC 
0.8468831 1.0220779 1.2904648 1.3986569

二つの変数 \(X,\,Y\) の共分散は次のように定義される。

\[ \mathrm{cov}(X, Y) = \frac{1}{n - 1}\sum_{i=1}^{n}(X_i - \mu_X)(Y_i - \mu_Y) \]

1月と7月の共分散を計算してみよう。

cov17 <- sum((df3$JAN - sum(df3$JAN) / ny) * (df3$JUL - sum(df3$JUL) / ny)) / (ny - 1)
print(cov17)
[1] -0.01234279

共分散は cov() や var() でも計算できる。

cov(df3$JAN, df3$JUL)
[1] -0.01234279
var(df3$JAN, df3$JUL)
[1] -0.01234279

相関係数は共分散を二つの変数の標準偏差で割ったものである。

\[ r(X, Y) = \frac{\sum_{i=1}^{n}(X_i - \mu_X)(Y_i - \mu_Y)} { \sqrt{\sum_{i=1}^{n}(X_i - \mu_X)^2} \sqrt{\sum_{i=1}^{n}(Y_i - \mu_Y)^2}} \]

sd1 <- sd(df3$JAN)
sd7 <- sd(df3$JUL)
cov17 / (sd1 * sd7)
[1] -0.01395136

相関係数は cor() で計算できる。

cor(df3$JAN, df3$JUL)
[1] -0.01395136
Note練習
  1. 各月の間の相関係数を計算しよう。
Code 
cor(df3)
            JAN         FEB        MAR       APR       MAY       JUN
JAN  1.00000000  0.95519690 0.84808644 0.6642305 0.3785150 0.1252096
FEB  0.95519690  1.00000000 0.92570217 0.7823039 0.5202152 0.2649108
MAR  0.84808644  0.92570217 1.00000000 0.8979479 0.6566527 0.4098332
APR  0.66423049  0.78230393 0.89794793 1.0000000 0.8697025 0.6437206
MAY  0.37851500  0.52021518 0.65665270 0.8697025 1.0000000 0.8940658
JUN  0.12520959  0.26491081 0.40983318 0.6437206 0.8940658 1.0000000
JUL -0.01395136  0.11375028 0.27532327 0.5113797 0.7719961 0.9447767
AUG -0.10950345  0.01347903 0.17641025 0.4086110 0.6873973 0.8737458
SEP -0.07239307  0.01795208 0.17651292 0.3986948 0.6593087 0.8324020
OCT -0.09605304 -0.01771509 0.12924250 0.3462847 0.6149149 0.7942519
NOV -0.10873771 -0.03287571 0.10850443 0.3168589 0.5808505 0.7653961
DEC -0.11063295 -0.04325644 0.09206688 0.2997562 0.5604379 0.7434683
            JUL         AUG         SEP         OCT         NOV         DEC
JAN -0.01395136 -0.10950345 -0.07239307 -0.09605304 -0.10873771 -0.11063295
FEB  0.11375028  0.01347903  0.01795208 -0.01771509 -0.03287571 -0.04325644
MAR  0.27532327  0.17641025  0.17651292  0.12924250  0.10850443  0.09206688
APR  0.51137974  0.40861096  0.39869484  0.34628475  0.31685891  0.29975616
MAY  0.77199614  0.68739726  0.65930875  0.61491493  0.58085046  0.56043787
JUN  0.94477674  0.87374581  0.83240199  0.79425186  0.76539609  0.74346825
JUL  1.00000000  0.95840399  0.89964779  0.85729773  0.82279643  0.79983536
AUG  0.95840399  1.00000000  0.94747041  0.90051886  0.87722215  0.85420992
SEP  0.89964779  0.94747041  1.00000000  0.96855530  0.94456695  0.92310953
OCT  0.85729773  0.90051886  0.96855530  1.00000000  0.97988242  0.96115110
NOV  0.82279643  0.87722215  0.94456695  0.97988242  1.00000000  0.98455493
DEC  0.79983536  0.85420992  0.92310953  0.96115110  0.98455493  1.00000000
  1. 異なる海域間の相関を計算しよう。
Code 
df4 <- read.table("https://www.data.jma.go.jp/tcc/tcc/products/elnino/index/sstindex/base_period_9120/Nino_4/sst")
df4 <- df4[-nrow(df4), ]
nino4 <- c(t(df4))
cor(nino3, nino4)
[1] 0.4827971