【Rによるデータサイエンス】線形回帰分析　〜　重回帰分析

前回は線形回帰分析のうち、線形単回帰分析について勉強した。今回は、重回帰分析について勉強する。

重回帰分析とは、説明変数が複数である回帰分析のこと。

次のように表現できることを想定する。

$y=a_{0}+a_{1}x_{1}+a_{2}x_{2}+\dots+a_{p}x_{p}+\epsilon$

回帰分析で求める回帰式は次。

$\hat{y}=a_{0}+a_{1}x_{1}+a_{2}x_{2}+\dots+a_{p}x_{p}$

説明変数のデータをX、目的変数をY、誤差をE、係数をAとすると、回帰モデルは次のようになる。

$Y=XA+E$

係数は最小2乗法などで求めることができて次のようになる。

$\hat{A}=(X^{t}X)^{-1}X^{t}Y$

導出はここに分かりやすいPDFがあるので参照。

◎ケーススタディ
Rに標準で用意されているデータairqualityを使用する。

> airquality
    Ozone Solar.R Wind Temp Month Day
1      41     190  7.4   67     5   1
2      36     118  8.0   72     5   2
3      12     149 12.6   74     5   3
4      18     313 11.5   62     5   4
5      NA      NA 14.3   56     5   5
6      28      NA 14.9   66     5   6
7      23     299  8.6   65     5   7

> colnames(airquality)
[1] "Ozone"   "Solar.R" "Wind"    "Temp"    "Month"   "Day"

それぞれ、オゾンの量（Ozone)、太陽放射量(Solar.R)、風力(Wind)、華氏温度（Temp)、月（Month)、日(Day)である。
ここでは、太陽放射量、風力、温度の値でオゾンの量を説明する重回帰モデル考える。

関数 pairs(行列) で各列同士の組合せ全てについて散布図を描いてみる。

> pairs(airquality[,1:4])
>

オゾンの量を目的変数、太陽放射量、風力、気温を説明変数とした重回帰分析を行う。

> air.lm<-lm(Ozone~Solar.R+Wind+Temp,data=airquality)
> summary(air.lm)

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp, data = airquality)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.485 -14.219  -3.551  10.097  95.619 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
Solar.R       0.05982    0.02319   2.580  0.01124 *  
Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
Temp          1.65209    0.25353   6.516 2.42e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.18 on 107 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared:  0.6059,    Adjusted R-squared:  0.5948 
F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

回帰係数(Coefficients)も表示されているが、小数点第2位までに丸めたものを取り出すと次のようになる。

> round(coefficients(air.lm),2)
(Intercept)     Solar.R        Wind        Temp 
     -64.34        0.06       -3.33        1.65

この結果から、オゾンの量を目的変数、太陽放射量、風力、気温を説明変数とした重回帰分析を行うと次の回帰式が求まる。

$Ozone=-64.34+0.06*{Solar.R}-3.33*Wind+1.65*Temp$

残差を視覚的に考察するために回帰診断図を作成する。

> par(mfrow=c(2,2))
> plot(air.lm)

Normal Q-Q（残差の正規Q-Qプロット）を見ると、点がほぼ直線上にのっているので、残差は正規分布に従っているとみなせる（たぶん・・・）。
また、残差の平方根プロットを見ると、0.5と1.5の間に収まっている。これで分散が説明変数や目的変数の値によらずに一定になっている（と見なせる？？？）。