重回帰分析
重回帰分析とは，説明変数が2変数以上の回帰分析である．多変量解析と呼ばれる分析手法群の中の一つの手法である．重回帰分析は，モデルが線形である場合，つまり簡単に言えば「y = ax1 + bx2 + c」のaとbとcを求めることに相当する．さらに数を増やして「y = ax1 + bx2 + cx3 + dx4 + e」という式のa, b, c, d, eも求められる．
重回帰分析の時，モデルの中の説明変数x1やx2は，実は(x1)^2など二乗の場合でも分析できる．あくまで係数との関係が線形であれば線形回帰として求められる．この場合，説明変数を入力する際にあらかじめ二乗したデータを作っておいて，それを入力する作業が必要になる．
Scikit-Learnでの重回帰分析
説明変数の単純回帰分析と同じくScikit-Learnライブラリでは，説明変数が多い場合の重回帰分析を全く同じようにで行うことができる．違うのは説明変数の数であり，前回の単純回帰分析では説明変数，目的変数に与えるNumpy配列のshapeは，両方とも(41, 1)であった．
重回帰分析の場合，説明変数のshapeが説明変数の個数の分だけ大きくなる．
つまり，y = ax1 + bx2 + cx3 + d という式の, a, b, c, dを求めたい時の，説明変数のshapeは(41, 3)ということになる．
単純泉系回帰分析と全く同じなので，以下のようにする．
import scipy
import sklearn
from sklearn import linear_model
ols_model = linear_model.LinearRegression()
old_model.fit(説明変数の2次元配列のDataFrame, 目的変数の配列DataFrame)
print('回帰係数(傾き):', ols_model.coef_)
print('切片:', ols_model.intercept_)
print('決定係数', ols_model.score(X, Y))statsmodelsでの重回帰分析
statsmodelsを使う場合でも，単純線形回帰分析と全く同じように行う．
statsmodelsを使った分散分析でも回帰分析が出てくるが，分散分析で使った回帰分析の関数は，式を与えることにより線形以外も可能なものなので，線形のみを前提にしている場合，それとは別の簡単な方法つまり単純線形回帰のプログラムと全く同じプログラムで実行することができる．
import statsmodels.api as sm
ols_model = sm.OLS(目的変数のDataFrame, sm.add_constant(説明変数のDataFrame))
ols_result = ols_model.fit()
print(ols_result.summary())DataFrameで複数の列を持ってくる
単純泉系回帰分析では説明変数は1次元配列であり，プログラム上は無理やり2次元配列にしたものをモデルに食わせるが，重回帰分析では説明変数が最初から2次元配列である．
CSVを読み込んだ後2次元配列テーブルを取り出す方法は，「「列の名前」をカンマで区切った配列」を使えば取得できる．
df = pandas.read_csv('CSVファイル名')
df = df.loc[:, ['平均気温(℃)', '日最高気温の平均(℃)']]
とすれば，平均気温の列と日最高気温の平均の2つの説明変数を格納したDataFrameができ，このDataFrameのshapeは(41, 2)となる．
後は全く同じように回帰分析クラスのインスタンスにこの説明変数と目的変数でfit()してやれば良い．
因果関係の推測
相関係数は，あくまで「相関」であり「因果」ではない，と相関係数の項で述べた．
statsmodelsの重回帰分析は，線形回帰分析の項で述べた通り，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をダウンロード」
データの加工
まずいつものようにCSVをGoogle Colabに読み込み，
from google.colab import files
uploaded = files.upload()ヘッダ行を飛ばしてpandas.DataFrame形式にしておく．
import numpy as np
import pandas as pd
from pandas import Series, DataFrame
df = pd.read_csv('hachioji-temp-data.csv', encoding='SHIFT-JIS', header=2)
# PandasのDataFrameという形式で読み込まれる
# 文字コードはSHIFT-JIS, 官公庁のオープンデータや，歴史が古い企業のデータはエクセルから書き出すものが多いからか，SHIFT-JISが多い．
# このheader=2は頭から2行分は不要な情報なので読み飛ばす，という意味
df.head()
# チェック用として，頭の部分だけ表示
statsmodelsの方が情報量が多いのでそちらを使う．読み込んだpandas.DataFrameから，不要な0行目をまず落とし，その後，目的変数と説明変数になる列を抽出する．
# まず0行目（品質情報や均質番号と書かれている行）を落とす
df = df.drop(0)
# 平均気温の列を目的変数とする
Y = df.loc[:, '平均気温(℃)']
# 品質情報，均質番号の列を落とす
X = df.loc[:, ['降水量の合計(mm)', '日照時間(時間)', '平均風速(m/s)']]
X.head()あとはstatsmodelsのOLSクラスを使って重回帰分析を実行する．
# 重回帰分析をこなう
import statsmodels.api as sm
ols_model = sm.OLS(Y, sm.add_constant(X))
ols_result = ols_model.fit()
print(ols_result.summary())結果は以下のようになる．
OLS Regression Results
==============================================================================
Dep. Variable: 平均気温(℃) R-squared: 0.636
Model: OLS Adj. R-squared: 0.607
Method: Least Squares F-statistic: 22.10
Date: Wed, 18 Mar 2020 Prob (F-statistic): 1.89e-08
Time: 07:23:17 Log-Likelihood: -48.103
No. Observations: 42 AIC: 104.2
Df Residuals: 38 BIC: 111.2
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 22.2742 0.643 34.648 0.000 20.973 23.576
降水量の合計(mm) -0.0003 0.001 -0.403 0.689 -0.002 0.001
日照時間(時間) 0.0138 0.003 5.393 0.000 0.009 0.019
平均風速(m/s) 0.5611 0.175 3.209 0.003 0.207 0.915
==============================================================================
Omnibus: 2.936 Durbin-Watson: 1.391
Prob(Omnibus): 0.230 Jarque-Bera (JB): 2.701
Skew: -0.552 Prob(JB): 0.259
Kurtosis: 2.431 Cond. No. 1.62e+03
==============================================================================考察と再度の実行
単純線形回帰分析と同じく，constが欠片で，coefが傾きになる．つまりこの結果は，
「8月の月平均気温 = -0.0003*降水量の月合計 + 0.0138*日照時間の月合計 + 0.05611*月平均風速 + 22.2742」を表している．
結果を見ると，降水量の傾きが非常に小さいことから，降水量は月平均気温にほとんど影響を与えていないことがわかる．実際にp値（Pr>|t|で表示されている）を見てみると，0.05を超えている．
回帰分析のp値はなんだったかを思い出してみると，帰無仮説：「説明変数は目的変数に影響を与えている」，対立仮説：「説明変数は目的変数に影響を与えているとは言えない」であった．
降水量のp値は0.05を大きく超えているので，「（八王子の場合）8月の降水量は月平均気温には影響を与えているとは言えない」と言うことができる．
では，降水量のデータを外して，再度重回帰分析をしてみよう．
目的変数は月平均気温まので変える必要がなく，元のpandas.DataFrameから抽出する列名指定を変えるだけである．
# 降水量を外して，再度重回帰分析を行う
X2 = df.loc[:, ['日照時間(時間)', '平均風速(m/s)']]
ols_model2 = sm.OLS(Y, sm.add_constant(X2))
ols_result2 = ols_model2.fit()
print(ols_result2.summary())結果は以下のようになる．
OLS Regression Results
==============================================================================
Dep. Variable: 平均気温(℃) R-squared: 0.634
Model: OLS Adj. R-squared: 0.615
Method: Least Squares F-statistic: 33.79
Date: Wed, 18 Mar 2020 Prob (F-statistic): 3.06e-09
Time: 07:40:14 Log-Likelihood: -48.192
No. Observations: 42 AIC: 102.4
Df Residuals: 39 BIC: 107.6
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 22.1272 0.524 42.254 0.000 21.068 23.186
日照時間(時間) 0.0142 0.002 6.124 0.000 0.010 0.019
平均風速(m/s) 0.5605 0.173 3.240 0.002 0.211 0.910
==============================================================================
Omnibus: 2.704 Durbin-Watson: 1.376
Prob(Omnibus): 0.259 Jarque-Bera (JB): 2.495
Skew: -0.527 Prob(JB): 0.287
Kurtosis: 2.439 Cond. No. 788.
==============================================================================降水量を外して計算したので，先ほどとは値が多少違っているが，無事に全てp値が低い値に収まっており，「日照時間と平均風速は月平均気温に影響を与えている」と言うことができる．
面白いのは，傾きが小さい日照時間の方がp値が小さいことである．これは単位というか物理量の違いがあるため．日照時間は「月に何時間日照があったか」という月の延時間，平均風速は気象庁オープンデータの場合「10分ごとの風速の月平均」であることから，傾きが違うのは当然とも言える．
重回帰分析の面白いところは，このように単位が異なっても，また前述のように「ダミー変数」を用いたとしても結果がでることである．ダミー変数は単位すらなく，0か1の値しか取らないが，p値をみることで「それが影響を与えているかどうか」をみることができる．重回帰分析の目的は，この「傾きを出す」ことではなく，p値を見て「影響を与えているか」を探ることである．
結果のR-squaredという項目を見てみよう．0.6341という値が表示されているが，これはR^2（Rの2乗）やR2という形でよくかかれる「決定係数」というものである．決定係数は「この得られた数式が，観測した値にどれだけフィットしているか」を表す値であり，0〜1の範囲を取り，1であれば非常にフィットしている，0であれば全くフィットしていない，という意味である．
例えば「p値は高いが，R^2が極端に低い時」には，lm関数に与えるformula，つまり式（「モデル」と呼ぶ）を見直してやると良い．今回与えたのは単純な線形多項式であるが，実は説明変数の2乗の項が含まれている二次関数かもしれない．
決定係数が0.5以下だと，数式が大幅に間違っているという判断をされることが多い．
課題
自身の身近な題材を取り上げ，データを取得，相関係数を算出，それを元に線形モデルを立て重回帰分析を行い，因果関係を推測せよ．
