気象情報解析特論　第11回　相関係数と回帰係数の検定と推定

担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from scipy import signal
from scipy.stats import t
import math

下準備：東京の気温と海面水温のデータセットの読み込み

今回は，東京と宇都宮の気温のデータを用います。

In [2]:
# 東京の気温のデータセットの読み込み
tokyo_temp = np.genfromtxt("Tokyo_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
tokyo = tokyo_temp[:, 2]

# 今回は1882年から2019年の140年分のデータを用いる
tokyo = tokyo[(1982 <= y)*(y <= 2019)]
m = m[(1982 <= y)*(y <= 2019)]
y = y[(1982 <= y)*(y <= 2019)]

# 気候値の計算（基準は140年の各月の平均）
tokyo_clim = np.zeros((12))
for mm in range(1, 13): 
    tokyo_clim[mm-1] = np.nanmean(tokyo[m==mm], 0)

# 偏差の計算
tokyoa = np.zeros((tokyo.shape))
for yy in range(1880, 2020):
    for mm in range(1, 13):
        tokyoa[(y==yy)*(m==mm)] = tokyo[(y==yy)*(m==mm)] - tokyo_clim[mm-1]

# 今回はデトレンドしたものとトレンド入りを両方使うので別の名前をつける
tokyoa_detrended = signal.detrend(tokyoa)

loadfile = 'detrended_ssta_OISST.npz' # デトレンド済みの海面水温偏差の入力ファイル名を定義
ssta_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
ssta = ssta_dataset['ssta'] # 海面水温偏差を変数sstaに保存
lon2 = ssta_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2 = ssta_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = ssta_dataset['y'] # 年(year)を変数yに保存
m = ssta_dataset['m'] # 月(month)を変数mに保存
# 今回は1982年から2019年のデータを用いる
ssta = ssta[:, :, (1982 <= y)*(y <= 2019)]
m = m[(1982 <= y)*(y <= 2019)]
y = y[(1982 <= y)*(y <= 2019)] # y自身のサイズが変わってしまうので，yは一番最後に書き換えないとダメ

# 領域平均をとる関数
def aave(west, east, south, north, var = ssta):
    var = var[(west<=lon2[:, 1])*(lon2[:, 1]<=east),:, :]
    var = var[:, (south<=lat2[1, :])*(lat2[1, :]<=north), :]
    aave_var = np.nanmean(np.nanmean(var, 0), 0)
    return aave_var

# 月別の時系列を2つ描画する関数
def plot_2_mon_time(time_series1, time_series2, lower = -5, upper = 5, init_year=1982, fin_year=2019):
    mon = np.arange(init_year, fin_year+1, 1/12)
    plt.plot(mon, time_series1)
    plt.plot(mon, time_series2, 'r')
    plt.plot(mon, 0*time_series1, 'k')
    plt.xlim(init_year, fin_year)
    plt.ylim(lower, upper)
    plt.show()
    
# 移動平均の関数
def running_mean(time_series, wn):
    b = np.ones(wn)/wn
    time_series_r = np.convolve(time_series, b, mode="same")
    n_conv = math.ceil(wn/2) #2で割って切り上げ
    time_series_r[0] = time_series_r[0] * wn/n_conv
    for i in range(1, n_conv):
        time_series_r[i] = time_series_r[i] * wn/(i+n_conv)
        time_series_r[-i] = time_series_r[-i] * wn/(i + n_conv - (wn % 2))
    return time_series_r

相関が良いのは偶然？

前回は，平均値の差のt検定を勉強しました。今回は，同様にまず相関係数の検定を勉強します。

早速ですが，日本東方沖の海面水温（35°N-45°N, 140°E-170°E）と，東京の気温にはどのくらい相関があるのか，調べてみましょう。

In [3]:
tokyoa_detrended_lp = running_mean(tokyoa_detrended, wn = 5)
tokyoa_detrended_lp = tokyoa_detrended_lp/np.std(tokyoa_detrended_lp)# 振幅を合わせるために正規化

nino34 = aave(140, 170, 35, 45)
nino34_lp = running_mean(nino34, wn = 5)
nino34_lp = nino34_lp/np.std(nino34_lp) # 振幅を合わせるためにこっちも正規化

print(np.corrcoef(tokyoa_detrended_lp, nino34_lp)[0, 1])
plot_2_mon_time(tokyoa_detrended_lp, nino34_lp)

0.5867021076914174

0.586という，非常に微妙な相関係数が出てきました。日本の気温と日本近海の海面水温なのだから，共通の原因で相関を持つのではないかと直観的には思いたくなります。一方で，東京の気温はヒートアイランドに影響されるかもしれないし，海には海流によって運ばれる熱もあるし...。

相関係数の検定

そんなときには，相関係数の検定をします。手順は，前回の「平均値の差の検定」のときとほぼ同じですが，使うべき公式が違います。

I. 有意水準：0.05（信頼度95%）

II. 帰無仮説「東京の気温と日本東方沖の海面水温の間の母相関係数はゼロである」

III. 対立仮説「東京の気温と日本東方沖の海面水温の間の母相関係数はゼロでない」

IV. 検定統計量tは，次の通り。

$$t = \frac{r} {\sqrt {\frac{1-r^2}{N^*-2}} }$$

ただし，$r$はサンプル相関係数，$N^*$は実効的サンプルサイズです。

相関係数の検定を行う場合は，実効的サンプルサイズの推定には以下の公式を用います（Bretherton et al., 1999)。

$$N^* = N \frac{1-r_1r_2}{1+r_1r_2}$$

ただし，$N^*$は「実効的サンプル数」，$N$は「元のサンプル数」，$r_1$と$r_2$はそれぞれの時系列の「ラグ1自己相関」です。

V. 棄却域は，次の通り。

$$|t| > t_{975, N^*-2}$$

ただし$t_{975, N^*-2}$は，統計的自由度$N^*-2$のt分布関数における97.5%の閾値（前回同様，両側に2.5%ずつ外れ値を仮定する両側検定）。

では実際にやってみます。

In [4]:
## 検定統計量を計算するのに必要な量を計算する関数。引数は時系列。
def calc_t_corr(x_1, x_2):
    # サンプル相関係数
    r = np.corrcoef(x_1, x_2)[0, 1]
    # 実効的サンプルサイズ
    N = x_1.shape[0]
    r_1 = np.corrcoef(x_1[0:-1], x_1[1:])[0, 1]
    r_2 = np.corrcoef(x_2[0:-1], x_2[1:])[0, 1]
    N_star = N * (1-r_1*r_2)/(1+r_1*r_2)
    # 検定統計量
    t_sample = r / (np.sqrt((1 - r**2)/(N_star - 2)))
    # 中を見たければコメントアウトを外す
    # print('r, N, N_star, t_sample are')
    # print([r, N, N_star, t_sample])
    return r, t_sample, N_star

r, t_sample, N_star = calc_t_corr(tokyoa_detrended_lp, nino34_lp)
print('\n t =', t_sample)

# 有意水準と統計的自由度から，棄却域の計算
t_975 = t.ppf(0.975, N_star - 2)
print('\n t_975 =', t_975)

 t = 4.277526082537193

 t_975 = 2.030402564059088

検定統計量$t=4.278$が棄却域に入っているため，帰無仮説を棄却し，対立仮説「東京の気温と日本東方沖の海面水温の間の母相関係数はゼロでない」を採択します。よって，東京の気温と日本東方沖の海面水温の間の相関係数0.59は，stochastisticに大きくなっているわけではなく，なんらかの意味があるのだろうと考えるのです。

相関係数を検定する際の注意点

ただし，相関係数を検定する際に注意しなければいけない点が二つあります。

アポステリオリ(a posteriori)とアプリオリ(a priori)

検定において有意であるとされれば，時系列の相関に何らかの意味があると言って良いのかと言うと，実は違います。「95%の信頼度において有意な相関係数の大きさ」を知ってから，その条件に合うような時系列を自然界から見つけて来るのは簡単だからです。

たとえば火星大接近に興味を持ち，どうしても地球-火星間の距離と降水量の間に何か関係を見出したいと思ったアマチュア天文家が，地上100カ所の観測所の降水量データを持って来て，地球-火星間の距離との相関係数を計算したとします。そうすると，95%の信頼度で検定している以上，100地点のうちに約5カ所ほどは，地球-火星間の距離と有意な相関係数を持つような降水パターンを示す観測所を見つけて来ることが出来ます。例えばその場所の一つがたまたま東京の大手町だったとして，そのアマチュア天文家が「大手町の雨は火星と関係あり！」という論文を書くことは，当然ながら許されません。

このように，先に有意な相関係数の大きさを「カンニング」してから，それに見合うような気候変数を探して来るのは アポステリオリな議論(a posteriori reasoning) と言って，少なくともこのような文脈における相関係数の検定では御法度です。ここで，アポステリオリとは「結果から原因に向かう」という意味です。

ではどうすれば良いのかというと，例えば「火星から大手町に降り注ぐ未知の物質が存在することが物理学の法則で予言された」などというような，火星と大手町の降水に有意な結果を得てもおかしくないという「統計以外に基づく期待」が必要なのです。これを，アプリオリな期待(a priori expectation) といいます。アプリオリというのは，アポステリオリの反義語で，「原因から結果に向かう」という意味です。

もしアプリオリな期待があれば，「なぜ中野坂上ではなく大手町なのか」という理由付けが出来ます。「大手町には他と違って〇〇というアプリオリな期待がある。その物理学によって選ばれし観測地点が，統計学的にも100地点から偶然に選ばれた5地点のうちの一つと一致したのだから，これは偶然ではないだろう」と言えるからこそ，火星と大手町の関係に意味が生まれるのです。

原因が先か，結果が先か，ここではそれが重要な意味をなしています。同じ「スキー場で恋に落ちた」という事象を考えたときでも，たまたまなんとなく気に入っていた相手にスキー場でバッタリ会ったら運命を感じるかもしれないですが，スキー場で出会った相手だけに絞って好きになるのでは「スキー場で恋に落ちる」のは当然であり，運命を感じることはありません。 ここで前者が運命を感じたのは，「アプリオリに期待していた相手」にスキー場で会えたからです。偶然性をよりどころとしている統計検定において，その「運命」をねじまげてしまっているのがアポステリオリな議論です。

ただし，アポステリオリな考え方自体が全て悪いというわけではないことは注意しましょう。実際，気候力学におけるほとんどの発見は，まずデータを見ておかしな現象を見つけてから，その現象を説得力を持って説明出来る物理を探すという「アポステリオリ的な思考」でなされています。その意味で，アポステリオリ自体が問題なのではなく， アポステリオリな現象に対してアプリオリな期待が必要な統計検定を行い，有意という結果を得てもなおアプリオリな期待を見つけられないまま，ただ闇雲に検定結果のみを信じようとするのが問題なのです。

発展的な注：実際，1972年に発見されたマッデン・ジュリアン振動(MJO)という熱帯大気の現象は，「アポステリオリ専用の検定」まで持ち出して完全にアポステリオリに発見され，そしてそれを説明する物理的なメカニズムは分かっていません。ここで当然，「アポステリオリ専用の検定」は，普段用いられているようなアプリオリ統計よりも，有意であるという結果を得るのが厳しくなっています。逆に，完全にアプリオリに発見された気候の現象は非常に少なく，19世紀後半に観測された「大気太陰潮汐」や，松野(1966)による「混合ロスビー重力波」「赤道ケルビン波」などが例として挙げられます。大気太陰潮汐については，ニュートンが万有引力を発見した際にまとめたPrincipia Mathematica (Newton, 1687)の中で，「月の引力が大気を引っ張り気圧を変化させるような現象があってもおかしくない」ことを示唆し，それを確かめる地上観測が2世紀後になされました。混合ロスビー重力波も，数式をいじっていたら導かれた現象が，その後すぐに柳井-丸山波やウォレス-カウスキー波として観測に引っかかりました。

相関関係は因果関係を意味しない（Correlation is not causation)

たとえば，南極の氷の量とアイスの売り上げには強い相関があることがアプリオリに期待されます。しかし，「南極の氷が増えた！よし，アイス買いに行こう！」という購買行動を取る人がいたとしたら常軌を逸しています。単に，北半球が夏になると，「南半球は冬になる」のと「アイスが売れる」のとが同時に起こるだけです。南極の氷の量とアイスの売り上げには，因果関係はないのです。

このように，共通要因CによってAとBが両方影響を受けることによって，AとB同士に因果関係がほとんどないにも関わらず相関係数が大きくなってしまう現象を疑似相関ということもあります。ただし，因果関係がないだけで，相関関係はあるので，「疑似」相関という名前はあまり良くない気がしています。

世の中には，面白い疑似相関をまとめたサイトがいっぱいあります（笑）。たとえば https://life-analyze24.com/gijisoukan-example/ あたりでしょうか。

相関係数の推定（フィッシャーのZ変換）

前回，モンテカルロ法によって「統計的推定」を行い，「真の回帰係数（＝母回帰係数）」があるだろう範囲を推定したのを覚えてますでしょうか。同様に，母相関係数がゼロでないと期待されるときには，フィッシャーのZ変換 という方法を使って母相関係数を推定することができます。

いま，サンプル相関係数を$r$とすると，

$$Z = \frac{1}{2}\ln\left(\frac{1+r}{1-r}\right)$$

という検定統計量$Z$は，次のような平均$\mu_z$と標準偏差$\sigma_z$を持つ正規分布に従うことが知られています。

$$\mu_z = \frac{1}{2}\ln\left(\frac{1+\rho}{1-\rho}\right), \ \ \ \sigma_z = \frac{1}{\sqrt{N^*-3}}$$

ただし，$\rho$は母相関係数です。つまり，まずサンプル相関係数を$Z$に変換することで平均$\mu_z$を推定し，その$\mu_z$を$Z$から相関係数に逆変換することによって母相関係数を推定すれば良いのです。

いま，$Z$が正規分布を満たすなら，$\mu_z$は95%（=1.96標準偏差）の信頼度で

$$Z - 1.96 \sigma_Z <\mu_z < Z + 1.96 \sigma_Z \iff \mu_{z1} <\mu_z < \mu_{z2}$$

の範囲に入ります。ただし，ここで計算した端点の値を$\mu_{z1}:= Z - 1.96 \sigma_Z$, $\mu_{z2}:=Z + 1.96 \sigma_Z$としました。

このときこれらの端点は，母相関係数の推定範囲の下端$\rho_1$と上端$\rho_2$との間に

$$\mu_{z1} = \frac{1}{2}\ln\left(\frac{1+\rho_1}{1-\rho_1}\right), \ \ \ \mu_{z2} = \frac{1}{2}\ln\left(\frac{1+\rho_2}{1-\rho_2}\right)$$

すなわち

$$\rho_1 = \tanh(\mu_{z1}), \ \ \ \rho_2 = \tanh(\mu_{z2})$$

という関係が成り立ちますので，母相関係数$\rho$は

$$\tanh(\mu_{z1}) < \rho < \tanh(\mu_{z2})$$

すなわち

$$\tanh(Z - 1.96 \sigma_Z) < \rho < \tanh(Z + 1.96 \sigma_Z)$$

と推定されます。

では，実際にやってみましょう。理論は複雑に見えたかもしれませんが，実践は「公式当てはめゲー」なので簡単です。

In [5]:
# Zとsigma_Zを計算
Z = 1/2*np.log((1+r)/(1-r))
sigma_Z = 1/np.sqrt(N_star-3)

# 下端と上端を計算
lower_bound = np.tanh(Z-1.96*sigma_Z)
upper_bound = np.tanh(Z+1.96*sigma_Z)

# 結果の表示
print(str(lower_bound) + ' < true correlation < ' + str(upper_bound))

0.3237076978648697 < true correlation < 0.7655390302478353

よって，東京の気温と日本東方沖の海面水温の間の母相関係数は，95%の信頼度で0.324から0.766の間に入っていると推定されます！

回帰係数（トレンド）の検定と推定

だいぶ，検定と推定の考え方に慣れてきたでしょうか？最後に総まとめとして，回帰係数（トレンド）の検定と推定を行ってみましょう。

1982年から2019年の間に，東京の気温は有意に上がったでしょうか？

In [6]:
mon = np.arange(1982, 2020, 1/12)
tokyoa = tokyoa[(1982 <= y)*(y <= 2019)]
plt.plot(mon, tokyoa)
[a, b] = np.polyfit(mon, tokyoa, 1) 
plt.plot(mon, a*mon + b, 'r')
plt.show()

回帰係数の検定

I. 有意水準：0.05（信頼度95%）

II. 帰無仮説「1982年から2019年の東京の気温のトレンドはゼロである」

III. 対立仮説「1982年から2019年の東京の気温のトレンドはゼロでない」

IV. 検定統計量tは，次の通り。

$$t = \frac{a}{\left(\frac{\sigma_e}{\sqrt{N}\sigma_x}\right)}$$

ただし，$a$はサンプル回帰係数（回帰直線は$y = ax+b$），$N$は（実際の）サンプルサイズ，$\sigma_x$は横軸の値$x$のサンプル標準偏差，$\sigma_e$は誤差項標準偏差

$$\sigma_e = \sigma_y\sqrt{\frac{N}{N^*-2}(1-r^2)}$$

です（発展的な注：この公式は普通の教科書に載っている誤差項標準偏差から導出できますが，結構大変でした）。ただし，$r$は$x$と$y$のサンプル相関係数，$\sigma_y$は縦軸の値$y$のサンプル標準偏差です。回帰係数の検定を行う場合は，実効的サンプルサイズ$N^*$の推定には以下の公式を用います（Bretherton et al., 1999)。

$$N^* = N \frac{1-r_{1y}}{1+r_{1y}}$$

ただし，$r_{1y}$は$y$の「ラグ1自己相関」です。

V. 棄却域は，次の通り。

$$|t| > t_{975, N^*-2}$$

ただし$t_{975, N^*-2}$は，統計的自由度$N^*-2$のt分布関数における97.5%の閾値（前回同様，両側に2.5%ずつ外れ値を仮定する両側検定）。

回帰係数の推定

回帰係数の推定範囲は，前述の記法をそのまま用いると

$$a - t_{975, N^*-2}\frac{\sigma_e}{\sqrt{N}\sigma_x} < \alpha < a + t_{975, N^*-2}\frac{\sigma_e}{\sqrt{N}\sigma_x}$$

です。ただし，$\alpha$が母回帰係数です。

発展的な注：この推定範囲において「下端>0または上端<0」が，検定の棄却域と同値です（簡単なので示してみましょう）。

さて，東京の気温のトレンド有意でしょうか？推定範囲はどの程度でしょうか？

In [7]:
## 検定統計量を計算するのに必要な量を計算する関数。引数は時系列。
def calc_t_regr(x, y):
    # サンプル回帰係数・相関係数
    a = np.polyfit(x, y, 1)[0]
    r = np.corrcoef(x, y)[0, 1]
    # 実効的サンプルサイズ
    N = x.shape[0]
    r_1y = np.corrcoef(y[0:-1], y[1:])[0, 1]
    N_star = N * (1-r_1y)/(1+r_1y)
    # サンプル標準偏差
    sigma_x = np.std(x)
    sigma_y = np.std(y)
    # 誤差項標準偏差
    sigma_e = sigma_y*np.sqrt(N/(N_star-2)*(1-r**2))  
    # 検定統計量
    t_sample = a / (sigma_e/(np.sqrt(N)*sigma_x))
    # 中を見たければコメントアウトを外す
    print('a, N, N_star, t_sample are')
    print([a, N, N_star, t_sample, np.sqrt(sigma_e/(np.sqrt(N)*sigma_x))])
    return a, t_sample, N_star, sigma_e, sigma_x

# 統計検定量などの計算
a, t_sample, N_star, sigma_e, sigma_x = calc_t_regr(mon, tokyoa)
print('\n t =', t_sample)

# 有意水準と統計的自由度から，棄却域の計算
t_975 = t.ppf(0.975, N_star - 2)
print('\n t_975 =', t_975)

# トレンドの推定区間
lower_bound = a - t_975*sigma_e/(np.sqrt(N_star)*sigma_x)
upper_bound = a + t_975*sigma_e/(np.sqrt(N_star)*sigma_x)

# 結果の表示
print('\n' + str(lower_bound) + ' < true regression < ' + str(upper_bound))

a, N, N_star, t_sample are
[0.024396895284859343, 456, 210.31330789959716, 3.617702458034104, 0.0821203626516759]

 t = 3.617702458034104

 t_975 = 1.9714173075372092

0.00482065322515203 < true regression < 0.043973137344566655

対立仮説「1982年から2019年の東京の気温のトレンドはゼロでない」を採択します。また，トレンドの推定範囲は，95%の信頼度で0.005 °C/年から0.044 °C/年の間に入っていると推定されます！

ここまでたどり着いたあなたなら，まだ検定や推定をやったことのない変数についても，統計学の教科書で検定統計量を調べれば自分でできそうですね！（？）

課題11（締切：一週間後の月曜10時40分まで）

この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。

A問題. 気象情報解析特論の感想を自由にお願いします。

B問題. 2000-2013年の期間では，東京の気温偏差に有意なトレンドはみとめられますか？

ヒント：この期間は，地球温暖化の停滞（ハイエイタス）と呼ばれている期間です。

C問題. あなたの趣味などに関する好きなデータを用いて，相関係数の検定を行ってください。また，フィッシャーのZ変換を用いて，相関係数の95%信頼区間を推定してください。

D問題. Mann-Kendall検定について調べてください。t検定と比べて，どのようなメリットがありますか？

提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論11，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。