重回帰分析
重回帰分析とは，説明変数が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以下だと，数式が大幅に間違っているという判断をされることが多い．
課題
自身の身近な題材を取り上げ，データを取得，相関係数を算出，それを元に線形モデルを立て重回帰分析を行い，因果関係を推測せよ．
