こーひーをぶんなぐれ！ - Work

この授業の目的

アート系学生が作品が，マーケティングなどの演習のためにデータ処理・分析を行う場合を想定して，統計などで頻繁に使われているR言語の基礎と，データ処理の方法と統計の基礎についてを学ぶ．数値を扱う内容・統計に関する内容はPython3版と対応させており，同じデータ，同じ内容である．

ネットコンテンツの取得やそのパース，文字列処理などについてはPythonの方がライブラリが揃っているためPythonを推奨する．逆にRは数学的なライブラリ（例えばグラフの描画や，あるパラメータを最適化するなど）がビルトインであるため，数学的な処理に関してはRを推奨する．（Pythonは追加パッケージがそれを担う）

（アルゴリズムを擬似コードで書きやすいものはPython，数式で書きやすいものはR，がそれぞれ有利，と言うイメージ）

本コンテンツはJupyter Notebookで動作するRを対象とする．Microsoftが提供するAzure NotebookのRは，PythonのGoogle Colaboratoryと同じくローカルに環境構築が一切いらないのでオススメ．

また，本コンテンツでは，読み込ませるデータの作成（オープンデータを整える）には表計算ソフトを対話的に用いる（Pyhotnの広義ノートの方ではプログラム的に処理している）．これは本Rの講義の方が，Pythonの講義に比べ，プログラミングの技能を獲得していない学生を対象としているからである．

コンテンツ

1. Rの基本，配列array，関数定義
2. Rのベクタとマトリクス
3. CSVのPlot，data.frameの取り扱い
4. Rのdata.frame用の高階関数apply
5. 散布図の描画
6. CSVデータからの相関係数の算出
7. CSVデータからの単純線形回帰分析
8. CSVデータからのt検定
9. CSVデータからの分散分析と多重比較
10. CSVデータからの重回帰分析と因果関係
11. CSVデータからの主成分分析
12. CSVデータからの線形二項分類
13. CSVデータからの因子分析

Azure NotebooksでRを使う方法

Microsoftが運営しているので，Microsoftアカウントを取得しておく．

Azure Notebooksのページに行き，右上のSign InからMSアカウントでサインイン（ログイン）する．

上の「My Projects」を選択．

「New Project」で授業用（もしくは自習用）のプロジェクトを作成する．

プロジェクト名と，プロジェクトID（URLになる．なんでも良い）を入れて，「Create」．

空のプロジェクトができる．要するに「フォルダ」である．

「＋」マークをクリックし，「Notebook」を選択．

ノートブックの名前をつけて，使用する言語「R」を選択して，「New」

新しいノートブックができたので，クリックする．

Rのコードが入力できるようになる．入力したら「Run」でその入力部分を実行できる．

本項でもGoogle Colaboratoryを使用する．（もちろんファイル読み込みの部分以外はJupyter Notebookでもよいし，一気にPythonのソースコードをまとめて書いて実行しても問題はない）

主成分分析

これまで「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」という4つの観測値を扱ってきた．この「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」という問題を考える．

ところが，観測値が4つある時，つまり説明変数が4つある時，2次元グラフはそれぞれの軸の組み合わせて6通り存在する．さらに3つ以上の観測値を軸にして見ることはできないため，グラフでぱっと見ではわからない．

そこで「主成分分析」という「次元削減」手法を用いる．

主成分分析は，説明変数が多数存在する時に，元のデータの特徴を十分に表しつつ，かつ次元が少ない説明変数のセットを「合成」する手法である．

実際のには，上図のようにただ足し合わせるだけではなく，変換行列を作成し，それによって次元削減をするが，その変換行列を作成する際に「元のデータの特徴を十分に表現できる」ことを目標として調整をしていく．

これは3DCGのレンダリングの時の3次元座標から2次元座標の変換行列（カメラの位置，アングル，拡大度等から算出される）を「3Dモデルのディティールが上手く掴める」ように調整するのに，イメージ的には近い．

十分に表現できるような「カメラの位置」「カメラの角度」「画角」「拡大」などをうまく調整し，3Dデータの特徴を十分に表現できるような2次元画像にするのと，主成分分析のための変換行列の作成は，似たようなもんである．

使用するデータと加工

今回はデータとして，相関係数を出した時と同じ「東日本の県庁所在地」の「1977年〜2018年」の「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」のCSVファイルを使う．

さらにそこから，盛岡，甲府，横浜に絞ってスタック形式のデータを作る．

GoogleスプレッドシートにCSVを読み込む．

例によって不要な行と列を削除．

左下の「＋」マークを押して，新規シートを追加．

シートの名前をわかりやすくつける．

最初のシートから盛岡の全データをコピーして，

新作ったシートに貼り付ける．（一番左の列は地点ラベルが入るので，2番目の列から貼り付ける）

盛岡のデータが貼り付けられた状態．

一番左の列名を「観測地点」にして，

1つのセルに「盛岡」と書き込んでコピーし，

データの他の行に貼り付ける

盛岡のデータは完成．

同じように甲府のデータもコピペして地点ラベルをつける．

横浜のデータも同様にする．

列名を，Rでも読み込める短いものにする．

これまでと同じく「ファイル」→「ウェブに公開」

スタック形式データのシートを選択し，CSVとして公開．公開するシートの選択を間違わないように．

import numpy as np
import pandas as pd

df = pd.read_csv('公開したURL')
df.head()

そのままグラフを描画してみる

まず，説明変数4つをx軸とy軸に取り，2次元のグラフを描画している．説明変数の組み合わせは合計6通りあるが，とりあえず3つの組み合わせでグラフを描画する．

# 生データを散布図でプロットしてみる
# 説明変数が4つで2次元に表示するためには，6通りのxとyの組み合わせがあるが，とりあえず3つだけ
import matplotlib.pyplot as plt

# 都市名ラベルの0番目を除いた説明変数4つ分で3つのグラフを描く
for i in range(1, 4, 1):
  for j, a_city in enumerate(['Morioka', 'Kofu', 'Yokohama']):
    # 42年分 * 3都市分なので，都市ごとに色を変えるため，データを分けてscatterを呼び出す
    x = df.iloc[j*42:(j+1)*42, i].astype(np.float32)
    y = df.iloc[j*42:(j+1)*42, i+1].astype(np.float32)
    plt.scatter(x, y, label=a_city)
  if i==1:
    plt.xlabel('Ave. of Temp in Aug.')
    plt.ylabel('Ave. of Max Temp. in Aug')
  elif i==2:
    plt.xlabel('Ave. of Max Temp. in Aug')
    plt.ylabel('Ave. of Min Temp. in Aug')
  else:
    plt.xlabel('Ave. of Min Temp. in Aug')
    plt.ylabel('Count of over 25 deg in Aug')
  plt.legend()
  plt.show()

X軸が「8月の平均気温」，Y軸が「8月の日最高気温の平均」

X軸が「8月の日最高気温の平均」，Y軸が「8月の日最低気温の平均」

X軸が「8月の日最低気温の平均」，Y軸が「25度を超えた日数」

今は都市ごとに色分けをしているため判別しやすいが，横浜はどのグラフでも分離して見える（線形判別分析が簡単にできそう，すなわち判別の「直線」を引くことが可能に見える，「線形分離可能」と呼ぶ．）のに対して，盛岡と甲府はかなり被っているように見える．（盛岡と甲府は同じ直線上に載っているように見える）

主成分分析は「このグラフが分離して見えるような表示になるような「合成変数」を作成すること」が目的である．合成変数は任意の数（次元）で作成することができる．合成変数を2つ（2次元）にしてしまえば，2次元グラフでも簡単に表示できるようになるはず．

主成分分析の実行

scikit-learnの中に入っているPCAクラスが使いやすい．主成分分析は「説明変数に対して行う」ものなので，PCAクラスでfitする時に食わせるデータは「説明変数のみ」である．

説明変数間で単位が違う場合「標準化」or「正規化」という作業をしなければならない（値のスケールをうまい具合に整えてくれるイメージ）．これは同じくscikit-learnの中に含まれているStandardScalerクラスで簡単に行うことができる．標準化するのは同じく説明変数のみであることに注意．

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# 単位が違う説明変数（3つの説明変数の単位は摂氏温度で，1つの単位は日数である）を「標準化」or「正規化」する
scaler = StandardScaler()
df_sc = scaler.fit_transform(df.loc[:, ['ave', 'max.ave', 'min.ave', '25over']])

# 合成変数を2次元で作ると指定
pca_model = PCA(n_components=2)
# データを学習して合成変数への変換行列を作る
pca_model.fit(df_sc)

# 変換行列ができたので，
# 4次元(8月の平均気温，8月の日最高気温の平均，8月の日最低気温の平均，25度を超えた日数）の説明変数を，
# 2次元の合成変数に変換する（合成する）
transformed_variables = pca_model.transform(df.loc[:,  ['ave', 'max.ave', 'min.ave', '25over']])
# 使いやすいようにデータフレーム化しておく
transformed_df = pd.DataFrame(transformed_variables)

# 合成変数2つの散布図を描く
for i, a_city in enumerate( ['Morioka', 'Kofu', 'Yokohama']):
  plt.scatter(transformed_df.iloc[i*42:(i+1)*42, 0], transformed_df.iloc[i*42:(i+1)*42, 1], label=a_city)
plt.legend()
plt.show()

結果は以下のようになる．

このグラフでは，X軸が主成分つまり合成変数0(Principal Component)0 , Y軸が主成分つまり合成変数1である．

合成変数では，盛岡，甲府，横浜が綺麗に分離しており，線形判別分析が容易になっている．この3都市の8月の気温に関しては，主成分分析 - 合成変数が有効であることを示している．

この状態で最初の目的である「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」を達成できたかを考えると，「直線の補助線を引くことで，「この線のこちら側が盛岡で，この線のこちら側が甲府で，この線のこちら側が横浜です」とうように明確に説明できるようになる」ので，目的を達成したと言える．

寄与率

主成分分析でもう一つ確認しなければならないのが「寄与率」である．主成分分析をした後に，モデルを保持している変数のexplained_variance_ratio_という変数の中に寄与率が入っている．

# 寄与率の表示
print(pca_model.explained_variance_ratio_)

寄与率の出力結果は以下の通り．

[0.97995861 0.01258219]

寄与率というのは，それぞれの主成分つまり合成変数が「データをどれだけ表現できているか」を表したものである．寄与率の合計が1.0に十分に近ければ，合成変数のみでデータを十分表現できていると言っても良い．

寄与率の出力は配列になっており，0番目が主成分合成変数0の寄与率，1番目が主成分合成変数1の寄与率である．合計が1に十分に近いので，この2次元だけで元データを十分に表現していると言って良い．つまり「4次元あった変数は2次元に落とし込むことが可能」ということである．

逆に合計寄与率が低かった場合，「この次元数では元のデータを表現できているとは言い難い」となる．この場合，合成変数の次元数を増やす必要がある．1次元ずつ増やしていって，寄与率が1に十分近くなる次元数を探す．次元数を増やしても合計寄与率が上がらなくなった場合はそこで止める．

主成分分析の合成変数間の相関係数

# 合成変数間の相関係数の表示
print(transformed_df.corr())

「合成変数間の相関係数」は，普通にデータフレームの内部の相関係数を出力すれば良い．

              0             1
0  1.000000e+00 -1.182588e-15
1 -1.182588e-15  1.000000e+00

合成変数間の相関係数は非常に小さい（ほぼ0）．というのも，主成分分析の内部計算は「合成変数間の相関係数を低くする」ように調整するからである．逆に言えば「相関が高いものをまとめる」という調整とも言える．

本項でもGoogle Colaboratoryを使用する．（もちろんファイル読み込みの部分以外はJupyter Notebookでもよいし，一気にPythonのソースコードをまとめて書いて実行しても問題はない）

主成分分析

これまで「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」という4つの観測値を扱ってきた．この「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」という問題を考える．

ところが，観測値が4つある時，つまり説明変数が4つある時，2次元グラフはそれぞれの軸の組み合わせて6通り存在する．さらに3つ以上の観測値を軸にして見ることはできないため，グラフでぱっと見ではわからない．

そこで「主成分分析」という「次元削減」手法を用いる．

主成分分析は，説明変数が多数存在する時に，元のデータの特徴を十分に表しつつ，かつ次元が少ない説明変数のセットを「合成」する手法である．

実際のには，上図のようにただ足し合わせるだけではなく，変換行列を作成し，それによって次元削減をするが，その変換行列を作成する際に「元のデータの特徴を十分に表現できる」ことを目標として調整をしていく．

これは3DCGのレンダリングの時の3次元座標から2次元座標の変換行列（カメラの位置，アングル，拡大度等から算出される）を「3Dモデルのディティールが上手く掴める」ように調整するのに，イメージ的には近い．

十分に表現できるような「カメラの位置」「カメラの角度」「画角」「拡大」などをうまく調整し，3Dデータの特徴を十分に表現できるような2次元画像にするのと，主成分分析のための変換行列の作成は，似たようなもんである．

使用するデータ

今回はデータとして，相関係数を出した時と同じ「東日本の県庁所在地」の「1977年〜2018年」の「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」のCSVファイルを使う．

コード例

1つ目のセルはいつもの通り，CSVをGoogle Colaboにアップロードする．

from google.colab import files
uploaded = files.upload()

2つ目のセルでは，相関係数を出した時と同じく，抜き出してDataFrameを構築する．

import numpy as np
import pandas as pd
from pandas import Series, DataFrame

df = pd.read_csv('tohoku-kanto-aug.csv', encoding='SHIFT-JIS', header=2)
# 前と同じ処理をして，抜き出しておく
for a_column in ['青森', '秋田', '盛岡', '山形', '仙台', '福島', '水戸', '宇都宮', 'さいたま', '前橋', '甲府', '千葉', '東京', '横浜']:
  for i in range(1, 11, 3):
    for j in range(0, 2, 1):
      df = df.drop(a_column + '.' + str(i + j), axis=1)

df = df.drop(1)

# DataFrameの中に欠損値がある場合は，0.0を代わりに入れておく
df = df.fillna(0)
df.head()

2つ目のセルの実行結果は以下のようになる．

3つ目のセルでは，その中から3つの都市を抜き出して，重回帰分析の項でやったようなスタック形式のデータを構築する．

全てのデータを変数ごとに配列につないで，順列に並んだデータフレームを構築する．ついでに都市名ラベルの列も作ってくっつけておく．そして全ての配列をvstack関数で繋ぐ．（vertical stack，つまり垂直方向にスタックする，の意味）

vstackしただけだと，想定とは逆に重なってしまうので，transposeして，配列の次元を入れ替えておく．

今回は，盛岡，甲府，横浜を選んだ．

# グラフ表示するためにローマ字ラベルを作っておく
all_city_list = (['Aomori', 'Akita', 'Morioka', 'Yamagata', 'Sendai', 'Fukushima', 'Mito', 'Utsunomiya', 'Saitama', 'Maebashi', 'Kofu', 'Chiba', 'Tokyo', 'Yokohama'])

# 盛岡，甲府，横浜を選ぶ
choosed_cities = (['Morioka', 'Kofu', 'Yokohama'])

# 選んだ都市の
# average_temp: 8月の平均気温
# max_average_temp: 8月の最高気温の月平均
# min_avaerage_temp: 8月の最低気温の月平均
# num_of_days_25over: 
# として，それぞれ42年分を，1つ配列としてくっつける

# 選んだ3都市の42年分のラベルを繋げておく　(42年分*3都市)
city_label = np.array(choosed_cities)
city_label = np.repeat(city_label, 42)

# データをつなげる
# まず最初に空の配列を作っておいて，
average_temp = np.array([])
max_average_temp = np.array([])
min_average_temp = np.array([])
num_of_days_25over = np.array([])

# 選んだ都市の，8月の平均気温，最高気温の月平均，最低気温の月平均，25度を超えた日数を
# それぞれ1つの配列につなぐ
for a_city in choosed_cities:
  i = all_city_list.index(a_city)
  average_temp = np.concatenate((average_temp, df.iloc[1:, 1+i*4]))
  max_average_temp = np.concatenate((max_average_temp, df.iloc[1:, 1+i*4+1]))
  min_average_temp = np.concatenate((min_average_temp, df.iloc[1:, 1+i*4+2]))
  num_of_days_25over = np.concatenate((num_of_days_25over, df.iloc[1:, 1+i*4+3]))

# vstackで重ねて，配列の次元を入れ替えて，DataFrameにしておく
df_X1 = pd.DataFrame(np.vstack([city_label, average_temp, max_average_temp, min_average_temp, num_of_days_25over]).transpose())
# DataFrameの列名をつけておく
df_X1 = df_X1.rename(columns={0: '都市名', 1:'8月の平均気温', 2:'8月の日最高気温の月平均', 3:'8月の日最低気温の月平均', 4:'25度を超えた日数'})
# DataFrameの冒頭を確認
df_X1.head()

3つ目のセルの実行結果は以下の通りである．

4つ目のセルでは，説明変数4つをx軸とy軸に取り，2次元のグラフを描画している．説明変数の組み合わせは合計6通りあるが，とりあえず3つの組み合わせでグラフを描画する．

# 生データを散布図でプロットしてみる
# 説明変数が4つで2次元に表示するためには，6通りのxとyの組み合わせがあるが，とりあえず3つだけ
import matplotlib.pyplot as plt

# 都市名ラベルの0番目を除いた説明変数4つ分で3つのグラフを描く
for i in range(1, 4, 1):
  for j, a_city in enumerate(choosed_cities):
    # 42年分 * 3都市分なので，都市ごとに色を変えるため，データを分けてscatterを呼び出す
    x = df_X1.iloc[j*42:(j+1)*42, i].astype(np.float32)
    y = df_X1.iloc[j*42:(j+1)*42, i+1].astype(np.float32)
    plt.scatter(x, y, label=a_city)
  if i==1:
    plt.xlabel('Ave. of Temp in Aug.')
    plt.ylabel('Ave. of Max Temp. in Aug')
  elif i==2:
    plt.xlabel('Ave. of Max Temp. in Aug')
    plt.ylabel('Ave. of Min Temp. in Aug')
  else:
    plt.xlabel('Ave. of Min Temp. in Aug')
    plt.ylabel('Count of over 25 deg in Aug')
  plt.legend()
  plt.show()

4つ目のセルの実行結果のグラフ3つは以下の通り．

X軸が「8月の平均気温」，Y軸が「8月の日最高気温の平均」

X軸が「8月の日最高気温の平均」，Y軸が「8月の日最低気温の平均」

X軸が「8月の日最低気温の平均」，Y軸が「25度を超えた日数」

今は都市ごとに色分けをしているため判別しやすいが，横浜はどのグラフでも分離して見える（線形判別分析が簡単にできそう，すなわち判別の「直線」を引くことが可能に見える，「線形分離可能」と呼ぶ．）のに対して，盛岡と甲府はかなり被っているように見える．（盛岡と甲府は同じ直線上に載っているように見える）

主成分分析は「このグラフが分離して見えるような表示になるような「合成変数」を作成すること」が目的である．合成変数は任意の数（次元）で作成することができる．合成変数を2つ（2次元）にしてしまえば，2次元グラフでも簡単に表示できるようになるはず．

5つ目のセルでは，主成分分析を行い，グラフを描画している．scikit-learnに簡易的な主成分分析キットがあるので，これを利用する．

説明変数間で単位が違う場合「標準化」or「正規化」という作業をしなければならない（値のスケールをうまい具合に整えてくれるイメージ）．これは同じくscikit-learnの中に含まれているStandardScalerクラスで簡単に行うことができる．標準化するのは同じく説明変数のみであることに注意．

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# 標準化する
scaler = StandardScaler()
df_X2 = scaler.fit_transform(df_X1.iloc[:, 1:])

# 合成変数を2次元で作る
pca_model = PCA(n_components=2)
pca_model.fit(df_X2)
# 変換行列ができたので，多次元の説明変数を2次元の合成変数に変換する
transformed_variables = pca_model.transform(df_X1.iloc[:, 1:])
# 使いやすいようにデータフレーム化しておく
transformed_df = pd.DataFrame(transformed_variables)
# 合成変数2つの散布図を描く
for i, a_city in enumerate(choosed_cities):
  plt.scatter(transformed_df.iloc[i*42:(i+1)*42, 0], transformed_df.iloc[i*42:(i+1)*42, 1], label=a_city)
plt.xlabel('Principal Component 0')
plt.ylabel('Principal Component 1')
plt.legend()
plt.show()

5つ目のセルで描画したグラフは以下の通り．

このグラフでは，X軸が主成分つまり合成変数0(Principal Component)0 , Y軸が主成分つまり合成変数1である．

合成変数では，盛岡，甲府，横浜が綺麗に分離しており，線形判別分析が容易になっている．この3都市の8月の気温に関しては，主成分分析 - 合成変数が有効であることを示している．

この状態で最初の目的である「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」を達成できたかを考えると，「直線の補助線を引くことで，「この線のこちら側が盛岡で，この線のこちら側が甲府で，この線のこちら側が横浜です」とうように明確に説明できるようになる」ので，目的を達成したと言える．

寄与率

主成分分析でもう一つ確認しなければならないのが「寄与率」である．主成分分析をした後に，モデルを保持している変数のexplained_variance_ratio_という変数の中に寄与率が入っている．

# 寄与率の表示
print(pca_model.explained_variance_ratio_)

寄与率の出力結果は以下の通り．

[0.97995861 0.01258219]

寄与率というのは，それぞれの主成分つまり合成変数が「データをどれだけ表現できているか」を表したものである．寄与率の合計が1.0に十分に近ければ，合成変数のみでデータを十分表現できていると言っても良い．

寄与率の出力は配列になっており，0番目が主成分合成変数0の寄与率，1番目が主成分合成変数1の寄与率である．合計が1に十分に近いので，この2次元だけで元データを十分に表現していると言って良い．つまり「4次元あった変数は2次元に落とし込むことが可能」ということである．

逆に合計寄与率が低かった場合，「この次元数では元のデータを表現できているとは言い難い」となる．この場合，合成変数の次元数を増やす必要がある．1次元ずつ増やしていって，寄与率が1に十分近くなる次元数を探す．次元数を増やしても合計寄与率が上がらなくなった場合はそこで止める．

主成分分析の合成変数間の相関係数

# 合成変数間の相関係数の表示
print(transformed_df.corr())

「合成変数間の相関係数」は，普通にデータフレームの内部の相関係数を出力すれば良い．

              0             1
0  1.000000e+00 -1.182588e-15
1 -1.182588e-15  1.000000e+00

合成変数間の相関係数は非常に小さい（ほぼ0）．というのも，主成分分析の内部計算は「合成変数間の相関係数を低くする」ように調整するからである．逆に言えば「相関が高いものをまとめる」という調整とも言える．

主成分分析

これまで「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」という4つの観測値を扱ってきた．この「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」という問題を考える．

ところが，観測値が4つある時，つまり説明変数が4つある時，2次元グラフはそれぞれの軸の組み合わせて6通り存在する．さらに3つ以上の観測値を軸にして見ることはできないため，グラフでぱっと見ではわからない．

そこで「主成分分析」という「次元削減」手法を用いる．

主成分分析は，説明変数が多数存在する時に，元のデータの特徴を十分に表しつつ，かつ次元が少ない説明変数のセットを「合成」する手法である．

実際のには，上図のようにただ足し合わせるだけではなく，変換行列を作成し，それによって次元削減をするが，その変換行列を作成する際に「元のデータの特徴を十分に表現できる」ことを目標として調整をしていく．

これは3DCGのレンダリングの時の3次元座標から2次元座標の変換行列（カメラの位置，アングル，拡大度等から算出される）を「3Dモデルのディティールが上手く掴める」ように調整するのに，イメージ的には近い．

十分に表現できるような「カメラの位置」「カメラの角度」「画角」「拡大」などをうまく調整し，3Dデータの特徴を十分に表現できるような2次元画像にするのと，主成分分析のための変換行列の作成は，似たようなもんである．

使用するデータと加工

今回はデータとして，相関係数を出した時と同じ「東日本の県庁所在地」の「1977年〜2018年」の「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」のCSVファイルを使う．

さらにそこから，盛岡，甲府，横浜に絞ってスタック形式のデータを作る．

GoogleスプレッドシートにCSVを読み込む．

例によって不要な行と列を削除．

左下の「＋」マークを押して，新規シートを追加．

シートの名前をわかりやすくつける．

最初のシートから盛岡の全データをコピーして，

新作ったシートに貼り付ける．（一番左の列は地点ラベルが入るので，2番目の列から貼り付ける）

盛岡のデータが貼り付けられた状態．

一番左の列名を「観測地点」にして，

1つのセルに「盛岡」と書き込んでコピーし，

データの他の行に貼り付ける

盛岡のデータは完成．

同じように甲府のデータもコピペして地点ラベルをつける．

横浜のデータも同様にする．

列名を，Rでも読み込める短いものにする．

新しく作ったスタック形式のシートを右クリックし「別のワークブックにコピー」→「新しいスプレッドシート」（Rは公開されたURLの関係上，単一ブックに単一シートでないと読み込めないようだ）

コピーされた新しいスプレッドシートを開く

新しいスプレッドシートの名前を整えて，

いつものようにCSVでウェブに公開してURLを取得しておく．

url <- '公開したURL'
df <- read.csv(url)
head(df)

いつものようにRで読み込み，冒頭数行を表示．

生のデータの散布図を描く

まず，説明変数4つをx軸とy軸に取り，2次元のグラフを描画している．説明変数の組み合わせは合計6通りあるが，とりあえず3つの組み合わせでグラフを描画する．

# スタック形式のデータを行指定でバラしてから，
library(stats)
morioka.df <- df[c(1:42), ]
kofu.df <- df[c(43:84), ]
yokohama.df <- df[c(85:126), ]


# X軸を平均気温，Y軸を日最高気温の月平均の散布図でプロットする
plot(morioka.df$ave, morioka.df$max.ave, col=2, xlim=c(18, 30), ylim=c(22, 36), ann=FALSE)
par(new=TRUE)
plot(kofu.df$ave, kofu.df$max.ave, col=3, xlim=c(18, 30), ylim=c(22, 36), ann=FALSE)
par(new=TRUE)
plot(yokohama.df$ave, kofu.df$max.ave, col=4, xlim=c(18, 30), ylim=c(22, 36), xlab='Ave. temp. of Aug.', ylab='Ave. of max temp. of aug.')

X軸が「8月の平均気温」，Y軸が「8月の日最高気温の平均」

# X軸を日最高気温の月平均, Y軸を日最低気温の月平均の散布図でプロットする
plot(morioka.df$max.ave, morioka.df$min.ave, col=2, xlim=c(24, 36), ylim=c(16, 26), ann=FALSE)
par(new=TRUE)
plot(kofu.df$max.ave, kofu.df$min.ave, col=3, xlim=c(24, 36), ylim=c(16, 26), ann=FALSE)
par(new=TRUE)
plot(yokohama.df$max.ave, kofu.df$min.ave, col=4, xlim=c(24, 36), ylim=c(16, 26), xlab='Ave of max temp. in Aug.', ylab='Ave. of Min temp. of Aug.')

X軸が「8月の日最高気温の平均」，Y軸が「8月の日最低気温の平均」

# X軸を日最低気温の月平均, Y軸を25度を超えた日数の散布図でプロットする
plot(morioka.df$min.ave, morioka.df$X25over, col=2, xlim=c(16, 26), ylim=c(0, 31), ann=FALSE)
par(new=TRUE)
plot(kofu.df$min.ave, kofu.df$X25over, col=3, xlim=c(16, 26), ylim=c(0, 31), ann=FALSE)
par(new=TRUE)
plot(yokohama.df$min.ave, kofu.df$X25over, col=4, xlim=c(16, 26), ylim=c(0, 31), xlab='Ave of min temp. in Aug.', ylab='Count of over 25 degree in Aug.')

X軸を日最低気温の月平均, Y軸を25度を超えた日数

今は都市ごとに色分けをしているため判別しやすいが，青(col=4)で示された横浜はどのグラフでも分離して見える（線形判別分析が簡単にできそう，すなわち判別の「直線」を引くことが可能に見える，「線形分離可能」と呼ぶ．）のに対して，盛岡と甲府はかなり被っているように見える．（盛岡と甲府は同じ直線上に載っているように見える）

主成分分析は「このグラフが分離して見えるような表示になるような「合成変数」を作成すること」が目的である．合成変数は任意の数（次元）で作成することができる．合成変数を2つ（2次元）にしてしまえば，2次元グラフでも簡単に表示できるようになるはず．

主成分分析の実行

prcomp関数を使う．

各変数の単位が異なっている（ここでは摂氏の温度と，日数の2つの単位が混在している）場合には，「正規化」or「標準化」を行う必要がある．prcomp関数の場合にはscale=TRUEオプションを入れれば良い．

# 各変数の単位（ここでは摂氏何度とか，日数の2つの単位が混在している）が異なっている時は，「正規化」or「標準化」を行わなければならない．
# prcomp関数の場合は，scaleオプションをTRUEにする

# まず目的変数と説明変数に分ける
city.name = df[, 1] # 1列目の全ての行をcity.nameに
city.param = df[, 2:5] # 2列目から5列目までの全ての行をcity.paramに

# 説明変数を正規化をして主成分析にかける
pca.results <- prcomp(city.param, scale=TRUE)
summary(pca.results)

# princomp という関数もあるが，これはサンプルサイズと変数の数に制限があるもの

結果表示は以下の通り．

Importance of components:
                          PC1     PC2     PC3     PC4
Standard deviation     1.9295 0.47388 0.22170 0.05831
Proportion of Variance 0.9307 0.05614 0.01229 0.00085
Cumulative Proportion  0.9307 0.98686 0.99915 1.00000

PCというのは主成分, 合成変数の英語であるPrincipal Componentの略語である．元の変数の数と同じだけ計算されるが，寄与率Proportion of Variance，および累積寄与率Cumulative Proportionに注目する．

寄与率というのは，それぞれの主成分つまり合成変数が「データをどれだけ表現できているか」を表したものである．寄与率の合計（累積寄与率）が1.0に十分に近ければ，合成変数のみでデータを十分表現できていると言っても良い．

今回は2番目の主成分までで累積寄与率が0.987に達しており，計が1に十分に近いので，この2次元だけで元データを十分に表現していると言って良い．つまり「4次元あった変数は2次元に落とし込むことが可能」ということである．

逆に合計寄与率が低かった場合，「この次元数では元のデータを表現できているとは言い難い」となる．この場合，合成変数して用いる次元数を増やす必要がある．1次元ずつ増やしていって，寄与率が1に十分近くなる次元数を探す．次元数を増やしても合計寄与率が上がらなくなった場合はそこで止める．

今回は第3主成分の寄与率が0.012であり，第2主成分までで累積寄与率が0.986つまり98.6%説明できることから，第2主成分までを用いる．（端的にいえば2次元の散布図1枚で描画できるようにするためである）

変換行列の取得と変換，第2主成分までを使った散布図の描画

predict関数に，主成分分析の結果と，新しいデータ（今回は古いデータをそのまま使っている）を食わせることで，主成分を計算することができる．

先にレンダリングに例えて変換行列と書いたが，predict関数はその「特徴をうまくとらえたレンダリング」を行うのと同じである．

レンダリングの結果を第2主成分（2番目の合成変数）までを使って（つまり2次元のみを使って），再び散布図に描画する．つまり「X軸が第1主成分，Y軸が第2主成分」である．

# 結果を使って主成分を計算し
PCs <- predict(pca.results, newdata=city.param)

# 第2主成分までを2次元に散布図で描く
morioka.pcs <- PCs[c(1:42), ]
kofu.pcs <- PCs[c(43:84), ]
yokohama.pcs <- PCs[c(85:126), ]

plot(morioka.pcs[, 1], morioka.pcs[, 2], col=2, xlim=c(-4.0, 3.0), ylim=c(-1.0, 1.0), ann=FALSE)
par(new=TRUE)
plot(kofu.pcs[, 1], kofu.pcs[, 2], col=3, xlim=c(-4.0, 3.0), ylim=c(-1.0, 1.0), ann=FALSE)
par(new=TRUE)
plot(yokohama.pcs[, 1], yokohama.pcs[, 2], col=4, xlim=c(-4.0, 3.0), ylim=c(-1.0, 1.0), xlab='PC1', ylab='PC2')

合成変数の散布図では，盛岡，甲府，横浜が綺麗に分離しており，線形判別分析が容易になっている．この3都市の8月の気温に関しては，主成分分析 - 合成変数が有効であることを示している．

この状態で最初の目的である「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」を達成できたかを考えると，「直線の補助線を引くことで，「この線のこちら側が盛岡で，この線のこちら側が甲府で，この線のこちら側が横浜です」とうように明確に説明できるようになる」ので，目的を達成したと言える．

重回帰分析

重回帰分析とは，説明変数が2変数以上の回帰分析である．多変量解析と呼ばれる分析手法群の中の一つの手法である．重回帰分析は，モデルが線形である場合，つまり簡単に言えば「y = ax1 + bx2 + c」のaとbとcを求めることに相当する．さらに数を増やして「y = ax1 + bx2 + cx3 + dx4 + e」という式のa, b, c, d, eも求められる．

重回帰分析の時，モデルの中の説明変数x1やx2は，実は(x1)^2など二乗の場合でも分析できる．あくまで係数との関係が線形であれば線形回帰として求められる．この場合，説明変数を入力する際にあらかじめ二乗したデータを作っておいて，それを入力する作業が必要になる．

因果関係の推測

相関係数は，あくまで「相関」であり「因果」ではない，と相関係数の項で述べた．

重回帰分析は，線形回帰分析の項で述べた通り，p値を出してくれる．このp値は設定した説明変数ごとに算出されることになるが，このp値は何を言っているか，の帰無仮説と対立仮説は以下の通りである．

• 帰無仮説: この説明変数は，目的変数に影響を与えない．
• 対立仮説: この説明変数は，目的変数に影響を与える変数である．

つまり回帰分析は「因果関係」を推測するために使うことができる．

単純線形回帰分析なら変数が一つなのでそれほど問題はないのだが，重回帰分析の場合，観測できるデータを全て説明変数として回帰分析を行なっても，その変数が「疑似相関」すなわち「因果関係はないのに相関係数が高い」だけの場合もあり，疑似相関関係の変数が説明変数に含まれているときちんとした結果が出ない．

そこで重回帰分析では，「慎重に因果関係を考察，疑似相関と思われる説明変数を除いて，何度も説明変数をとっかえひっかえしながら，結果（主にp値）を見ていく」ことで，因果関係を推測する．

ダミー説明変数

例えば分散分析の群のように，数字で表すことができない説明変数があり，これで重回帰分析を行いたい時は，ダミー説明変数という考え方を使う．

ダミー変数は，グループなどを表す変数の中身を1と0で表すものである．サンプルがA群とB群の2つ（もしくはその属性を持っているか，いないか）だった場合は，説明変数1つで中身は1か0かの2択でよいが，3つの群があった場合「A群である:0or1」「B群である:0or1」「C群である:0or1」の3つの説明変数になることに注意．

データの取得と加工

単純線形回帰分析と同じと同じ気象庁オープンデータの八王子のものを使用するが，取ってくるデータは違う．

「地点を選ぶ」タブで，八王子を選択．

「項目を選ぶ」タブで，「月別値」を選択し，中の「気温」タブで「月平均気温」を選択

「降水」タブで，「降水量の月合計」を選択．

「日照・日射」タブで，「日照時間の月合計」を選択

「風」タブで，「月平均風速」を選択．

「期間」タブで，1977年〜2018年の8月を選択．

「CSVをダウンロード」し，Googleスプレッドシートで開く．

これまで通りに，不要な行と列を削除．今回は八王子だけなので，地点の行が不要な代わりに，平均気温等の項目を列名として用いるために残す．今回は年月もデータとしては食わせないので削除．

半角少かっこはRでは列名として読み込めないようなので，列名の行の少かっこに囲まれた単位を削除する．

この状態で，前回のようにWebにCSV形式で公開し，URLを取得しておく．

url <- 'Web公開したURL'
df <- read.csv(url, header=TRUE)
head(df)

重回帰分析を実行

単純線形回帰分析の時とほぼ同じだが，説明変数として与えるのが「複数の列」であることに注意．

result.lm <- lm(formula='平均気温 ~ .', data=df)
summary(result.lm)

formulaで与えてやる式だが，CSVを読み込んだdata.frame中の「平均気温」列を目的変数とするのは単純線形回帰分析と同じ．~で繋いだあとの「.」は，「data.frameの中の，「平均気温」列を除いた列」という意味になる．

つまり，このdata.frameの中には「平均気温」「降水量の合計」「日照時間」「平均風速」の4つの列があるが，目的変数として「平均気温」を採用して，残りの3列は説明変数として用いる，と言う意味である．

結果を受けた変数をsummary関数に食わせることで，結果を出力することができる．lm結果とsummaryを実行して出力される結果は以下の通り．

Call:
lm(formula = "平均気温 ~ .", data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7948 -0.6126  0.1871  0.5318  1.1036 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.274174   0.642870  34.648  < 2e-16 ***
降水量の合計 -0.000349   0.000866  -0.403  0.68922    
日照時間      0.013820   0.002563   5.393 3.87e-06 ***
平均風速      0.561129   0.174873   3.209  0.00271 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7996 on 38 degrees of freedom
Multiple R-squared:  0.6357,	Adjusted R-squared:  0.6069 
F-statistic:  22.1 on 3 and 38 DF,  p-value: 1.89e-08

結果の考察と再度の実行

単純線形回帰分析と同じく，Interceptが欠片で，Estimateが傾きになる．つまりこの結果は，

「8月の月平均気温 = -0.000349*降水量の月合計 + 0.013820*日照時間の月合計 + 0.0561129*月平均風速 + 22.274174」を表している．

結果を見ると，降水量の傾きが非常に小さいことから，降水量は月平均気温にほとんど影響を与えていないことがわかる．実際にp値（Pr>|t|で表示されている）を見てみると，0.05を超えている．

回帰分析のp値はなんだったかを思い出してみると，帰無仮説：「説明変数は目的変数に影響を与えている」，対立仮説：「説明変数は目的変数に影響を与えているとは言えない」であった．

降水量のp値は0.05を大きく超えているので，「（八王子の場合）8月の降水量は月平均気温には影響を与えているとは言えない」と言うことができる．

では，降水量のデータを外して，再度重回帰分析をしてみよう．

result.lm2 <- lm(formula='平均気温 ~ 日照時間 + 平均風速', data=df)
summary(result.lm2)

今度のformulaは，降水量を外しているので，「目的変数として指定した平均気温の列以外」である「.」は使えないため，列名を直接入れている．

この結果は以下のようになる．

Call:
lm(formula = "平均気温 ~ 日照時間 + 平均風速", data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8292 -0.6002  0.1751  0.5433  1.0780 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 22.127196   0.523674  42.254  < 2e-16 ***
日照時間     0.014232   0.002324   6.124 3.49e-07 ***
平均風速     0.560468   0.172977   3.240  0.00245 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.791 on 39 degrees of freedom
Multiple R-squared:  0.6341,	Adjusted R-squared:  0.6153 
F-statistic: 33.79 on 2 and 39 DF,  p-value: 3.058e-09

降水量を外して計算したので，先ほどとは値が多少違っているが，無事に全てp値が低い値に収まっており，「日照時間と平均風速は月平均気温に影響を与えている」と言うことができる．

面白いのは，傾きが小さい日照時間の方がp値が小さいことである．これは単位というか物理量の違いがあるため．日照時間は「月に何時間日照があったか」という月の延時間，平均風速は気象庁オープンデータの場合「10分ごとの風速の月平均」であることから，傾きが違うのは当然とも言える．

重回帰分析の面白いところは，このように単位が異なっても，また前述のように「ダミー変数」を用いたとしても結果がでることである．ダミー変数は単位すらなく，0か1の値しか取らないが，p値をみることで「それが影響を与えているかどうか」をみることができる．重回帰分析の目的は，この「傾きを出す」ことではなく，p値を見て「影響を与えているか」を探ることである．

「日照時間が気温に影響を与えること」は明白であり，さらにp値も十分に低いことから，日照時間と月の平均気温の間には因果関係があると言っても良いだろう．

さらに「風が強くなると気温に大きく影響を与える」ことが統計的に示されたが，直接の因果関係があるとは考えづらく，間にならんらかの要素を挟んでいる（例えば日照時間に影響を与える雲の量は風速によって影響を受ける）と考えられる．

つまり，p値が直接因果関係を示すわけではなく「因果関係を考えるために統計的なp値を使う」ことができる．

結果のR-squaredという項目を見てみよう．0.6341という値が表示されているが，これはR^2（Rの2乗）やR2という形でよくかかれる「決定係数」というものである．決定係数は「この得られた数式が，観測した値にどれだけフィットしているか」を表す値であり，0〜1の範囲を取り，1であれば非常にフィットしている，0であれば全くフィットしていない，という意味である．

例えば「p値は高いが，R^2が極端に低い時」には，lm関数に与えるformula，つまり式（「モデル」と呼ぶ）を見直してやると良い．今回与えたのは単純な線形多項式であるが，実は説明変数の2乗の項が含まれている二次関数かもしれない．

決定係数が0.5以下だと，数式が大幅に間違っているという判断をされることが多い．

課題

自身の身近な題材を取り上げ，データを取得，相関係数を算出，それを元に線形モデルを立て重回帰分析を行い，因果関係を推測せよ．