データサイエンス 2026 / 06 海面水温との比較

06. 海面水温との比較（大気海洋相互作用に関連して）

SOI と El Niño 監視海域の SST anomaly を比較し、相関係数、回帰式、決定係数、ラグ相関を通して大気海洋相互作用を理解する。

主題
SOI と NINO 海域 SST anomaly の関係
入力データ
SOI_timeseries_updated.txt / sstoi.indices
見るポイント
相関、回帰、R²、共通期間、ラグ相関

1. この回でやること

前回までに SOI を計算し、移動平均や自己相関からその時間変動を見てきた。今回は、SOI を 大気側の指標、El Niño 監視海域の SST anomaly を 海洋側の指標 とみなし、両者の関係を調べる。

  1. SOI 時系列と SST anomaly を読み込む
  2. 共通期間だけを取り出す
  3. NINO1+2, NINO3, NINO4, NINO3.4 の各海域について相関係数を求める
  4. 最も関係の強い海域について回帰式と決定係数を求める
  5. 散布図と回帰直線を描く
  6. ラグ相関を計算し、時間差のある関係を調べる
この回の狙いは、ENSO を「大気」と「海洋」が結びついた変動として理解することである。

2. 相関係数・回帰式・決定係数とは何か

2.1 相関係数

相関係数は、2つの変数がどれくらい一緒に増減するかを表す量である。値は -1 から 1 の範囲をとる。

r = cov(X, Y) / (σX σY)

ここで cov(X, Y) は共分散、σ は標準偏差である。共分散だけだと単位に依存するので、標準偏差で割って無次元化している。

2.2 最小二乗法による回帰式

SST anomaly を説明変数 X、SOI を目的変数 Y として、

Y = aX + b

の形の直線を求める。ここで a は傾き、b は切片である。最小二乗法では、観測値と回帰直線の差の二乗和が最小になるように a と b を決める。

2.3 決定係数

決定係数 R² は、その回帰式がどれだけ当てはまっているかを表す指標である。

R² = 1 - SSR / TSS

ここで SSR は残差平方和、TSS は全平方和である。単回帰では、相関係数の 2 乗にかなり近い意味を持つ。

量意味
相関係数 r2つの変数の増減の向きと強さ
回帰式 Y = aX + bX から Y を近似する直線
決定係数 R²その直線がどれだけ当てはまるか
相関が高いことと、因果関係が直接あることは同じではない。ただし、SOI と SST のように物理的背景がある場合は、相関の意味を考察する価値が高い。

3. 使うデータと注意点

今回は 2 つのファイルを使う。

注意すべきことは、SOI のデータは 2025 年 3 月までであり、SST 側はそれ以降まで含む場合があることである。したがって、解析には必ず 共通期間だけ を使う。

解析では、見た目で同じ時期だと思い込まず、日付列を作ってから機械的に共通期間へそろえるのが安全である。
STEP 1

4. SOI と SST データを読む

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

soi_df = pd.read_csv("SOI_timeseries_updated.txt", sep=r"\s+")
soi_df = soi_df.dropna(subset=["SOI"]).reset_index(drop=True)

soi_df["DATE"] = pd.to_datetime(dict(
    year=soi_df["YEAR"].astype(int),
    month=soi_df["MONTH"].astype(int),
    day=1
))

sst_df = pd.read_csv("sstoi.indices", sep=______, engine="python")
print(sst_df.columns)

sst_df["DATE"] = pd.to_datetime(dict(
    year=sst_df["YR"].astype(int),
    month=sst_df["MON"].astype(int),
    day=1
))

このブロックで大事なこと

  • sep=r"\s+"：空白区切りのテキストを読む
  • dropna(subset=["SOI"])：SOI が欠損の行を落とす
  • pd.to_datetime(...)：年・月から日付を作る
  • print(sst_df.columns)：列名を先に確認する。思い込みで書かない
STEP 2

5. 共通期間にそろえる

SST ファイルでは anomaly 列が ANOM, ANOM.1, ANOM.2, ANOM.3 のように並んでいる。これを分かりやすい列名に直してから、SOI と日付で結合する。

sst_use = sst_df[[
    "DATE",
    "ANOM",
    "ANOM.1",
    "ANOM.2",
    "ANOM.3"
]].copy()

sst_use = sst_use.rename(columns={
    "ANOM": "NINO1+2_ANOM",
    "ANOM.1": "NINO3_ANOM",
    "ANOM.2": "NINO4_ANOM",
    "ANOM.3": "NINO3.4_ANOM"
})

df = pd.merge(soi_df, sst_use, on="DATE", how=______)

print(df["DATE"].min())
print(df["DATE"].max())
print(len(df))

ここでのポイント

  • rename(columns={...})：列名を分かりやすくする
  • pd.merge(..., on="DATE")：日付で 2 つの表を結合する
  • how="inner"：両方に存在する日付だけを残す
STEP 3

6. 各海域の相関係数・回帰式・決定係数

まずは 1 つの関数にまとめて、どの海域が最も SOI と強く結びついているかを調べる。

def calc_stats(x, y):
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    mask = np.isfinite(x) & np.isfinite(y)
    x = x[mask]
    y = y[mask]

    r = np.corrcoef(x, y)[0, 1]
    a, b = np.polyfit(x, y, 1)

    y_pred = a * x + b
    ss_res = np.sum((y - y_pred) ** 2)
    ss_tot = np.sum((y - np.mean(y)) ** 2)
    r2 = 1 - ss_res / ss_tot

    return r, a, b, r2, x, y, y_pred
regions = ["NINO1+2_ANOM", "NINO3_ANOM", "NINO4_ANOM", "NINO3.4_ANOM"]
results = []

for region in regions:
    r, a, b, r2, x_valid, y_valid, y_pred = calc_stats(df[region], df["SOI"])
    results.append({
        "Region": region,
        "Correlation_r": r,
        "Slope_a": a,
        "Intercept_b": b,
        "R2": r2
    })

result_df = pd.DataFrame(results)
result_df["abs_r"] = result_df["Correlation_r"].abs()
result_df = result_df.sort_values("abs_r", ascending=______)
print(result_df)
ENSO 監視では一般に NINO3.4 が強い関係を示すことが多い。結果がどうなるか、自分で確かめることが大事である。
STEP 4

7. 散布図と回帰直線

最も相関の強い海域として NINO3.4 を選び、SOI との散布図を描く。回帰直線も一緒に描き、図中に相関係数と決定係数を入れる。

target = "NINO3.4_ANOM"
r, a, b, r2, x_valid, y_valid, y_pred = calc_stats(df[target], df["SOI"])

plt.figure(figsize=(6, 6))
plt.scatter(x_valid, y_valid, alpha=0.5)

xx = np.linspace(np.min(x_valid), np.max(x_valid), 100)
yy = a * xx + b
plt.plot(xx, yy, linewidth=2)

plt.xlabel("NINO3.4 SST anomaly")
plt.ylabel("SOI")
plt.title("NINO3.4 anomaly vs SOI")
plt.grid(True, alpha=0.3)

plt.text(
    0.05, 0.95,
    f"r = {r:.2f}\nR² = {r2:.2f}\nSOI = {a:.2f} × SST + {b:.2f}",
    transform=plt.gca().transAxes,
    verticalalignment="top"
)

plt.tight_layout()
plt.show()

実際の結果を以下に示す。

図：NINO3.4 と SOI の関係

図：NINO3.4 海面水温偏差と SOI の散布図および回帰直線

考えてみよう
  • なぜ相関係数は負になりやすいのか
  • 回帰直線の傾きが負とはどういう意味か
  • R² が 1 に近くないのはなぜか
STEP 5

8. ラグ相関（フィルタなし）

今度は、SST と SOI を時間的にずらしながら相関を求める。これがラグ相関である。ここでは lag < 0 なら SOI が先行、lag > 0 なら SST が先行と読む。

def cross_correlation(x, y, max_lag):
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    lags = np.arange(-max_lag, max_lag + 1)
    cc = np.full(len(lags), np.nan)

    for i, lag in enumerate(lags):
        if lag < 0:
            x1 = x[-lag:]
            y1 = y[:len(y) + lag]
        elif lag > 0:
            x1 = x[:len(x) - lag]
            y1 = y[lag:]
        else:
            x1 = x
            y1 = y

        mask = np.isfinite(x1) & np.isfinite(y1)
        if np.sum(mask) > 2:
            cc[i] = np.corrcoef(x1[mask], y1[mask])[0, 1]

    return lags, cc
x_raw = df["NINO3.4_ANOM"].to_numpy()
y_raw = df["SOI"].to_numpy()

mask_raw = np.isfinite(x_raw) & np.isfinite(y_raw)
x_raw = x_raw[mask_raw]
y_raw = y_raw[mask_raw]

lags_raw, cc_raw = cross_correlation(x_raw, y_raw, max_lag=24)
imin_raw = np.nanargmin(cc_raw)
lag_min_raw = lags_raw[imin_raw]
corr_min_raw = cc_raw[imin_raw]
STEP 6

9. ラグ相関（13か月移動平均後）

生データのままだと月ごとのゆらぎが大きいので、13か月移動平均をかけた後でもラグ相関を求める。ここでは 121か月移動平均は使わず、短周期のノイズを弱めるための 13か月平滑化だけを使う。

def moving_average(x, window):
    return pd.Series(x).rolling(window=window, center=True).mean().to_numpy()

x_smooth = moving_average(df["NINO3.4_ANOM"].to_numpy(), 13)
y_smooth = moving_average(df["SOI"].to_numpy(), ____)

mask_smooth = np.isfinite(x_smooth) & np.isfinite(y_smooth)
x_smooth = x_smooth[mask_smooth]
y_smooth = y_smooth[mask_smooth]

lags_smooth, cc_smooth = cross_correlation(x_smooth, y_smooth, max_lag=24)
imin_smooth = np.nanargmin(cc_smooth)
lag_min_smooth = lags_smooth[imin_smooth]
corr_min_smooth = cc_smooth[imin_smooth]
強いフィルタをかけすぎると、月単位のラグ情報そのものがぼやける。今回は、短周期のばらつきを弱める最小限の平滑化として 13 か月移動平均だけを使う。
STEP 7

10. フィルタ前後の比較

最後に、フィルタなしと 13か月移動平均後のラグ相関を 1 枚に重ねて描く。

plt.figure(figsize=(8, 6))
plt.plot(lags_raw, cc_raw, linestyle="--", color="tab:blue", label="Raw")
plt.plot(lags_smooth, cc_smooth, color="tab:orange", linewidth=2, label="13-month smoothed")

plt.axhline(0, color="k", linewidth=0.8)
plt.axvline(0, color="k", linewidth=0.8)

plt.scatter(lag_min_raw, corr_min_raw, color="tab:blue", s=70, zorder=3)
plt.scatter(lag_min_smooth, corr_min_smooth, color="tab:orange", s=70, zorder=3)

plt.annotate(
    f"raw\nlag={lag_min_raw}\nr={corr_min_raw:.2f}",
    xy=(lag_min_raw, corr_min_raw),
    xytext=(lag_min_raw - 8, corr_min_raw - 0.12),
    arrowprops=dict(arrowstyle="->", color="tab:blue"),
    color="tab:blue",
    fontsize=12
)

plt.annotate(
    f"smoothed\nlag={lag_min_smooth}\nr={corr_min_smooth:.2f}",
    xy=(lag_min_smooth, corr_min_smooth),
    xytext=(10, corr_min_smooth + 0.40),
    arrowprops=dict(arrowstyle="->", color="tab:orange"),
    color="tab:orange",
    fontsize=12
)

plt.xlabel("Lag (month)")
plt.ylabel("Correlation")
plt.title("Cross-correlation (Raw vs Smoothed)")
plt.legend(loc="lower left")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

計算結果を図で確認する。

図：ラグ相関（Raw vs Smoothed）

図：フィルタ前後のラグ相関比較。平滑化により関係が明瞭になる

比較図では、フィルタ後のほうが谷がなめらかになり、lag = 0 付近の関係が見やすくなる。これは関係が新しく生まれたのではなく、もともとの関係が見えやすくなったと考えるべきである。

11. まとめ

12. 解答の表示

まずは上の穴埋めコードを自分で考えること。確認したい場合だけ、パスワードを入力して解答と完全スクリプトを表示する。