【R言語】重回帰分析 - ナナシの技術

前回の単回帰分析に続いて今回は重回帰分析です。重回帰分析は一つの目的変数を複数の説明変数で回帰式を表す分析となります。解釈容易性が高い予測モデルを構築することができるため、予測精度よりも解釈がし易いモデルを選ぶ場合には選択肢に入るかもしれませんね。

対象データ
前処理
重回帰分析
おわりに
1. 関連

対象データ

今回の対象データはR言語の標準データセットに含まれる「airquality」を用います。

「airquality」データセットはニューヨークのルーズベルト島で観測したオゾン濃度の時系列データセットです。

データの構成は以下のようになっています。

# ニューヨークのルーズベルト島で観測したオゾン濃度の時系列データセットを使用
libary(datasets)
data(airquality)

#データの構造と項目の一覧を確認する
str(airquality)

'data.frame':	153 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

データセットの一覧を見てわかる通り、観測データと日時データから構成されていることがわかります。

前処理

「airquality」データセットをこのまま重回帰分析するには問題があります。

データセットを一覧を見ると「NA」の表記があります。これは欠損値を表しています。

欠損値のままでは重回帰分析を行うことができないため、どうにか処理をする必要があります。今回は楽な方法として、欠損値がある行を削除していきたいと思います。

また、日付データも今回の重回帰分析では不要になるため削除しておきます。

#欠損値が含まれる行を削除
DFsub <- na.omit(airquality)
str(DFsub)#欠損値行を削除しても111行であるためこのまま分析を進める
#時系列情報を削除
DFsub <- DFsub[, -5:-6]

重回帰分析

相関関係

それでは重回帰分析を行っていきます。今回は目的変数を「Ozone」として分析します。

まずは、相関関係を確認していきましょう。

#ライブラリGGallyのggpairs()で散布図マトリクスを描く
library(ggplot2)
library(GGally)
ggpairs(DFsub,                        #すべての列を含める
        aes(alpha=0.5),                    #透明度
        upper=list(continuous=wrap("cor", size=9))) +
  #相関係数の文字サイズ
  theme(axis.text =element_text(size=9),  #軸の文字サイズ
        strip.text=element_text(size=9))  #項目の文字サイズ

それぞれの変数間の相関係数を観察すると、問題がありそうな説明変数が一つ確認できます。「Temp」は「Solar.R」と「Wind」との間に少なからずも相関関係があるようです。

なぜ相関関係があると問題なのかというと、重回帰分析を行う上では以下の前提があるためです。

独立性（データそれぞれが独立）
等分散性（説明変数にかかわらず分散が一定）
正規性（誤差自体が正規分布している）
線形性（説明変数と目的変数の関係は直線で近似できる）

説明変数間に相関関係があると独立性が保たれません。このように説明変数間に相関関係があり、偏回帰係数の信頼性がなくなることを多重共線性と言います。

多重共線性

次にこの多重共線性の確認をしていきます。VIF統計量を求めることで、各説明変数間の多重共線性について調べることができます。

VIF統計量を求めるには「car」ライブラリをインストールして、vif関数に重回帰モデルを引数に渡すことで算出されます。

library(car)
#3要素用いた重回帰分析
LM3 <- lm(Ozone ~ Solar.R + Wind + Temp, data=DFsub)
vif(LM3)

> vif(LM3)
 Solar.R     Wind     Temp 
1.095253 1.329070 1.431367

一般的にVIF統計量が10以下であれば多重共線性がないとされます。今回の場合は問題はなさそうです。

しかし、今回はあえてTemp変数を減らした重回帰モデルも作成してきたいと思います。

#2要素用いた重回帰分析
LM2 <- lm( Ozone ~ Solar.R + Wind, data=DFsub)

２つの重回帰モデルを比べてどちらが有意な重回帰モデルなのかを比べていきたいと思います。

作成した重回帰モデルのSummaryを出力します。

summary(LM3)
summary(LM2)

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

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
Multiple R-squared:  0.6059,	Adjusted R-squared:  0.5948 
F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

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

Residuals:
    Min      1Q  Median      3Q     Max 
-45.651 -18.164  -5.959  18.514  85.237 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 77.24604    9.06751   8.519 1.05e-13 ***
Solar.R      0.10035    0.02628   3.819 0.000224 ***
Wind        -5.40180    0.67324  -8.024 1.34e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.92 on 108 degrees of freedom
Multiple R-squared:  0.4495,	Adjusted R-squared:  0.4393 
F-statistic: 44.09 on 2 and 108 DF,  p-value: 1.003e-14

２つの重回帰モデルのSumarryを確認すると、P値についてはどちらも各説明変数は0.05以下になっているため有意なようです。

しかし、決定係数については３変数の方が0.6059。2変数の方は0.4495になっています。３変数の方が目的変数を説明できているようです。

次にAIC（赤池情報量基準）を用いて比べてみます。AICはモデルの当てはまり度を表す統計量です。小さい値の方が良いとされていますが、相対的な評価として用いられるため、ある一定値以下が望ましいといったことはありません。

#AIC(赤池情報量基準)を使ってモデルを比較
#モデルの複雑さとデータへの適合のバランスを見る指標
#小さい方がよい
AIC(LM3)
AIC(LM2)

> AIC(LM3)
[1] 998.7171
> AIC(LM2)
[1] 1033.816

３変数の重回帰モデルの方が値が小さくなりました。AICを用いた比較でも３変数重回帰モデルの方が良いモデルと言えそうです。

補足として、BIC（ベイズ情報基準量）でも比較してみます。AICと同様に値が小さいほうがモデルに適合していると言えます。

#参考：BIC(ベイズ情報量基準)を使った比較
#指標値の見方はAICと同じ
BIC(LM3)
BIC(LM2)

> BIC(LM3)
[1] 1012.265
> BIC(LM2)
[1] 1044.654

こちらの結果もAICと同様に３変数重回帰モデルの方が良いと言えそうです。

モデル評価

次に予測値と残差を用いて図からモデルの比較を行っていきたいと思います。予測値は重回帰モデルのを格納した変数のfitted.valuesに格納されています。また、残差についてはresidualsに格納されているのでこれらを用いていきます。

#予測値はモデル作成の際にfitted.valuesとして格納されている
#実測値との差（残差）は同様にresidualsとして格納されている
head(DFsub$Ozone)       #実測値
head(LM3$fitted.values) #予測値
head(LM3$residuals)     #残差（予測値－実測値）

> head(DFsub$Ozone)       #実測値
[1] 41 36 12 18 23 19
> head(LM3$fitted.values) #予測値
        1         2         3         4         7         8 
33.045483 34.998710 24.822814 18.475226 32.261431 -6.949919 
> head(LM3$residuals)     #残差（予測値－実測値）
          1           2           3           4           7           8 
  7.9545175   1.0012902 -12.8228139  -0.4752262  -9.2614315  25.9499188

まずは、横軸を実測値、縦軸をモデル上の予測値にした散布図を作成していきます。散布図を作成する目的は、実測値と予測値がどれだけ乖離しているのかを視覚的に確認するためです。実測値と予測値が完全に一致していればデータは対角線の上に一直線に並ぶはずです。

#実測値を横軸に、予測値を縦軸にプロットする
#　モデルが完璧にデータに適合していれば、点は対角線（y=x）
#　の上に並ぶはず（対角線からの縦方向のずれが残差を表す）
ggplot()+
  geom_point(aes(x=DFsub$Ozone, y=LM3$fitted.values),
             colour="orange", size=4, #LM3をオレンジで表示
             shape=16, alpha=.6 ) +　 #shape=16：●
  geom_point(aes(x=DFsub$Ozone, y=LM2$fitted.values),
             colour="brown",  size=3, #LM2をブラウンで表示
             shape=17, alpha=.6 ) +   #shape=17：▲
  xlab("実測値") +                    #軸ラベル
  ylab("モデル上の予測値(元のデータに基づく予測値)") +
  stat_function(colour="black", fun=function(x)x)   #y=Xに沿って線を引く

作成した散布図を見ると、実測値の値が大きくなるにつれて両モデル共に対角線より下に乖離しています。つまり、予測された値が実際より小さくなるようですね。

次に２つのモデルの残差を密度プロットで表してみましょう。残差の値は重回帰モデルのresidualsに格納されています。

#残差の分布を密度プロットで比較
ggplot()+
  geom_density( aes(x=LM3$residuals), #LM3の残差
                color="orange",       #枠線をオレンジ
                fill ="orange",       #塗りもオレンジ
                alpha=0.2  ) +        #透明度を指定
  geom_density( aes(x=LM2$residuals), #LM2の残差
                color="brown",        #枠線をブラウン
                fill ="brown",        #塗りもブラウン
                alpha=0.2  ) +        #透明度を指定
  xlab("残差(モデル上の予測値－実測値)")  #軸ラベル

残差は０に近いほうが良く、分布が左右に広がらずに中央に寄っている方が良いと言えます。逆に山が低く裾が広がっているとモデルとデータの乖離が大きいということになります。以上のことにより、３変数重回帰モデルの優れていると言えそうです。

折角なので残差を二乗した密度プロットも出力してみましょう。

#残差の2乗の分布を密度プロットで比較
ggplot()+
  geom_density( aes(x=LM3$residuals^2), #LM3の残差(2乗)
                color="orange",         #枠線をオレンジ
                fill ="orange",         #塗りもオレンジ
                alpha=0.2  ) +          #透明度を指定
  geom_density( aes(x=LM2$residuals^2), #LM2の残差(2乗)
                color="brown",          #枠線をブラウン
                fill ="brown",          #塗りもブラウン
                alpha=0.2  ) +          #透明度を指定
  xlab("残差の2乗")  #軸ラベル