環境情報論　第8回 回帰係数と相関係数

担当教員: 神山 翼 (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

東京と宇都宮の気温を読み込んで，解析の「下準備」（前回までの復習）

東京の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Tokyo_temp.csv

宇都宮の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Utsunomiya_temp.csv

（出典は第3回に記述）

In [2]:
tokyo_temp = np.genfromtxt("Tokyo_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
utsu_temp = np.genfromtxt("Utsunomiya_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
tokyo = tokyo_temp[:, 2]
utsu = utsu_temp[:, 2]

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

# 気候値の計算
tokyoc= np.zeros((12))
utsuc= np.zeros((12))
for mm in range(1, 13): 
    tokyoc[mm-1] = np.nanmean(tokyo[m==mm], 0)
    utsuc[mm-1] = np.nanmean(utsu[m==mm], 0)

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

# 今回は温暖化には興味がないのでデトレンド
tokyoa = signal.detrend(tokyoa)
utsua = signal.detrend(utsua)

気温偏差の描画

In [3]:
mon = np.arange(1990, 2020, 1/12)
plt.plot(mon,tokyoa)
plt.title('Temperature Anomalies in Tokyo')
Out[3]:
Text(0.5, 1.0, 'Temperature Anomalies in Tokyo')
In [4]:
plt.plot(mon,utsua)
plt.title('Temperature Anomalies in Utsunomiya')
Out[4]:
Text(0.5, 1.0, 'Temperature Anomalies in Utsunomiya')

二つの地点の気温偏差は，なんだかとても似ていますね。間違って同じデータをダウンロードしている？心配なので，重ねて描いてみましょう。

In [5]:
plt.plot(mon, tokyoa, 'r') #東京は赤で描画
plt.plot(mon, utsua)
Out[5]:
[<matplotlib.lines.Line2D at 0x1a1c504910>]

非常に似ているけど，ちょっとずれていますね。それにしても似ています。

気候値を引くところで，季節変化（「夏は暑く，冬は寒い」）はすでに除去していますので，似ているのは季節のせいではありません。ただ，東京が寒い年は，宇都宮も寒いでしょうから，まぁそんなこともあるでしょう。

素朴な疑問として，東京が1℃上がると，宇都宮は何℃上がっているのでしょうか？

また，東京と宇都宮の気温の変動は，どのくらい比例関係にあるといえるのでしょうか？

今日は，そのような「二つの時系列の関係」を定量化するためのお話です。

回帰係数

まず，東京の気温偏差を$x$, 宇都宮の気温偏差を$y$として，$xy$平面上に散布図を描いてみましょう。

In [6]:
import matplotlib as mpl

def draw_scatter(x, y, xname, yname):
    plt.scatter(x, y)
    plt.axhline(y=0, xmin=-10, xmax=10, color='k')
    plt.axvline(x=0, ymin=-10, ymax=10, color='k')
    plt.xlim(-5, 5)
    plt.ylim(-5, 5)
    plt.xlabel(xname)
    plt.ylabel(yname)
    plt.gca().set_aspect('equal', adjustable='box') # 縦横比を1:1にする

draw_scatter(tokyoa, utsua, 'Tokyo (℃)', 'Utsunomiya (℃)')

トレンドを計算するところで既に出てきていますが，点（$x_i$, $y_i$）($i = 1, 2, 3, ..., N$)という点列が与えられたとき（$N$はデータ数），「最小二乗法」という方法によって，点全体の散らばりを最もよく近似する直線$y = ax + b$（回帰直線といいます）における傾き$a$と切片$b$を求めるには，[a, b] = np.polyfit(x, y, 1)という関数を用います。

In [7]:
[a, b] = np.polyfit(tokyoa, utsua, 1)

直線$y=ax+b$を重ねて，赤色で描画してみましょう。

In [8]:
draw_scatter(tokyoa, utsua, 'Tokyo (℃)', 'Utsunomiya (℃)')
x = np.array([-10, 10]) # グラフの端点である5より大きければなんでもよい
plt.plot(x, a*x + b, 'r') # 点(-10, a*(-10)+b)と点(10, a*(10)+b)を繋げる直線を引く
Out[8]:
[<matplotlib.lines.Line2D at 0x109665490>]

特に，この傾き$a$を回帰係数と呼びます。ここでは横軸も縦軸も℃ですので，回帰係数の単位は「℃/℃」または無単位です。

今回の回帰係数の値は

In [9]:
a
Out[9]:
0.9620619831036824

ですので，東京の気温が1℃上がるとき，宇都宮は0.96℃上がっているという解釈になります。

もちろんここまでは「このデータと解析に絶対の信頼を置くとすれば」の話ですので，実際にはこの推定には不確かさがあります。

ただ，ここではそういう細かいことは忘れて数字をそのままの意味で読み取ると，宇都宮の方が東京よりもちょっとだけ，年々の気温変化が緩やかであると言えます。

ちなみに切片$b$については，今回の場合

In [10]:
b
Out[10]:
-8.649122905631364e-17

となっており，これは数学的には厳密にゼロです。

切片$b$は，「一年を通して東京より宇都宮の方が寒い」のような効果を表す数ですので，これは気候値を除去した時点で差っ引かれているからです。

偏差同士の回帰直線は，常に$y=ax$の形をしているということは，頭の片隅に入れておいてください。

相関係数

回帰係数は，最もよく近似する直線の傾きを与えたにすぎず，「どれだけよく直線にのっているか」を表してはいません。どんなに点がバラバラでも，わずかな偏りを見つけてエイヤッと引けてしまうのが回帰直線です。

それに対し，「どれだけよく直線にのっているか」「どれだけ比例関係にあるといえるか」「エイヤッと回帰直線を引くことがどの程度許されるのか」を表す係数を相関係数といい，慣例的には$r$で表します。

点（$x_i$, $y_i$）($i = 1, 2, 3, ..., N$)という点列が与えられたとき（$N$はデータ数），それらの相関係数（$y=ax$にどれだけ点がよく乗っているか）を求めるにはnp.corrcoef(x, y)[1, 0]という関数を使います。

In [11]:
np.corrcoef(tokyoa, utsua)[1, 0]
Out[11]:
0.9385239274981371

相関係数の取りうる範囲は，$-1 \leq r \leq 1$です。なぜなら，数学的には相関係数は偏差を並べた$N$次元ベクトル$\vec{x} := (x_1, x_2, ..., x_N)$と$\vec{y} := (y_1, y_2, ..., y_N)$のなす角$\theta$のコサインで定義されるからです。

$r := cos \theta = \frac{\vec{x}\cdot\vec{y}}{|\vec{x}||\vec{y}|}$

In [12]:
np.dot(tokyoa, utsua)/np.sqrt(np.dot(tokyoa, tokyoa)*np.dot(utsua, utsua))
# np.dot(x, y)はxとyの内積，np.sqrtは平方根
Out[12]:
0.938523927498137

確かにnp.corrcoefで計算したのと同じ値になっていますね。

数学的定義から順を追って考えれば分かると思いますが，点列（$x_i$, $y_i$）の相関係数$r$について

・$r=1$となるケース（$\vec{y} = a\vec{x}, a>0$）は，点列（$x_i$, $y_i$）が完璧に回帰直線$y=ax$にのり，かつ回帰係数$a$が正であるとき

・$0<r<1$となるケースは，回帰直線$y=ax$にパラパラのり，かつ回帰係数$a$が正であるとき（正の相関があるという）

・$r=0$となるケース（$\vec{x} \cdot \vec{y} = 0$）は，$x$と$y$が全く比例関係にないとき（無相関である，直交しているなどという）

・$-1<r<0$となるケースは，回帰直線$y=ax$にパラパラのり，かつ回帰係数$a$が負であるとき（負の相関があるという）

・$r=-1$となるケース（$\vec{y} = a\vec{x}, a<0$）は，点列（$x_i$, $y_i$）が完璧に回帰直線$y=ax$にのり，かつ回帰係数$a$が負であるとき

という風に場合分けされます。

東京と宇都宮の相関係数は0.94ですので，「東京が暖かいときは宇都宮も暖かい」というかなり強い正の相関があるということになります。

回帰係数と相関係数の例

それでは，いくつか例を見てみましょう。

( i ) $a=2, r = 1$のとき（完璧な比例関係）

In [13]:
np.corrcoef(tokyoa, 2*tokyoa)[1, 0]
Out[13]:
0.9999999999999998
In [14]:
draw_scatter(tokyoa, 2*tokyoa, 'Tokyo (℃)', '2*Tokyo (℃)')

( ii ) $a \simeq -1, r \simeq -0.9$のとき（強い負の相関=$x$が大きいほど$y$が小さい）

In [15]:
np.corrcoef(tokyoa, -utsua)[1, 0]
Out[15]:
-0.9385239274981371
In [16]:
draw_scatter(tokyoa, -utsua, 'Tokyo (℃)', '-Utsunomiya (℃)')

( iii ) $a \simeq -0.33, r \simeq -0.9$のとき（iiと相関は同じだが回帰の絶対値が小さいケース）

In [17]:
np.corrcoef(tokyoa, -utsua/3)[1, 0]
Out[17]:
-0.9385239274981366
In [18]:
draw_scatter(tokyoa, -utsua/3, 'Tokyo (℃)', '-(1/3)*Utsunomiya (℃)')

( iii ) $a \simeq -1, r\simeq -0.5$のとき（iiと回帰は同じだが相関が弱いケース）

In [19]:
# 標準正規分布に従う乱数を二種類発生させる
noiseA = (np.random.randn(360))
noiseB = (np.random.randn(360))

np.corrcoef(tokyoa + noiseA, -tokyoa + noiseB)[1, 0]
Out[19]:
-0.5239986991095672
In [20]:
draw_scatter(tokyoa + noiseA, -tokyoa + noiseB, 'Tokyo (℃) + noiseA', '-Tokyo (℃) + noiseB')

( iv ) $r \simeq 0$のとき（$x$と$y$が無関係ならば無相関）

In [21]:
np.corrcoef(noiseA, noiseB) [1, 0]
Out[21]:
0.007959812303424206
In [22]:
draw_scatter(noiseA, noiseB, 'noise', 'noise')

( v ) $r \simeq 0$のとき（無相関でも無関係とは限らない Part1）

In [23]:
np.corrcoef(noiseA, noiseA**2-2)[1, 0]
Out[23]:
-0.09366670446995293
In [24]:
draw_scatter(noiseA, noiseA**2-2, 'noise', 'noise^2-2')

$x$と$y$がそれぞれ東京と宇都宮の気温だと思うことにして，時系列にプロットしてみましょう。

In [25]:
plt.plot(mon, noiseA, 'r') # xを赤でぬる
plt.plot(mon, noiseA**2-2)
plt.xlim([2005, 2010]) #正直見た目だけではxとyの関係はわからない
Out[25]:
(2005, 2010)

( vi ) $r \simeq 0$のとき（無相関でも無関係とは限らない Part2）

In [26]:
np.corrcoef(3*np.cos(2*np.pi*noiseA), 3*np.sin(2*np.pi*noiseA))[1, 0]
Out[26]:
0.00983275880650397
In [27]:
draw_scatter(3*np.cos(2*np.pi*noiseA), 3*np.sin(2*np.pi*noiseA), '3*cos(2*PI*noise)', '3*sin(2*PI*noise)')
In [28]:
plt.plot(mon, 3*np.cos(2*np.pi*noiseA), 'r')
plt.plot(mon, 3*np.sin(2*np.pi*noiseA))
plt.xlim([2000, 2005]) #正直見た目だけではxとyの関係はわからない
Out[28]:
(2000, 2005)

いかがでしたか？（まとめサイト風）

教訓は，相関係数は「$x$と$y$が無関係か」を確かめる指標ではなく，「$x$と$y$は比例関係にあるか」の指標にすぎないということです。

課題8（締切：次回の授業開始前まで）

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

A問題. 回帰係数・相関係数とは何か，自分の言葉で簡単に説明してください。答案は短くて良いです。

B問題. 気象に限らなくて良いので，身近な例で回帰係数と相関係数について説明してください。

C問題. あなたの趣味などの自由なデータについて，上記の例のように散布図を描き，回帰係数や相関係数を計算してみてください。

D問題. アパレル業界で働くミカさんとカズハさんは，東京の湿度が夏物ワンピースの売り上げに影響するのかどうかに興味を持ち，湿度と夏物ワンピース売上額の月別データを数年分用意しました。ミカさんは，湿度と夏物ワンピース売上額のデータについてそのまま相関係数を計算しました。カズハさんは，湿度と夏物ワンピース売上額のデータについて，まずそれぞれの気候値とトレンドを除去して，偏差のデータにしてから相関係数を計算しました。

  1. どちらの相関係数が高く出るでしょうか。

  2. どちらの解析が本目的に対して，より妥当な解析でしょうか。

あなたの考えを，理由とともに自分の言葉で説明してください。また，（衣服に限らなくていいので）適当な「季節変動する経済活動」に関するデータをインターネット上で探し，上記の状況を模したデータ解析の結果を示してください（気象データのダウンロード方法は第3回C問題参照）。

ヒント：ミカさんの解析を，湿度の代わりに南極海氷面積で行うとどうなるでしょうか。

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