Rで単回帰から重回帰-#2【重回帰】(終)

[前回](https://qiita.com/H_Y_J/items/28063931d2f193e41f3a)の単回帰の流れから重回帰もまとめてしまう。

#重回帰

重回帰ではニューヨークの大気データの「airquality」を使います。

ニューヨークの大気データ154日
オゾン量、太陽放射強度、風、気温のデータが入っています。
たとえば体に降り注ぐオゾンの量は簡単に測定できないものとして、
太陽の強さと風と気温という身近な指標から予測するモデルを作りたいと思います。

時系列に並んでいるデータなので、時期によって影響は無いのかも含めてとりあえずplotです。

```R
aq=airquality


times=as.POSIXct(paste0("1973-",airquality$Month,"-",airquality$Day),format="%Y-%m-%d")
aq$time=times

y_limit<-c(0,350)
matplot(aq$time, aq$Temp, 
        type="n", lty=0, xaxt="n", xlab="", 
        las=1, ylim=y_limit, main="plot ts")

matpoints(aq$time, aq$Temp, type="l", lty=1,col="red")
matpoints(aq$time, aq$Wind, type="l", lty=1,col="blue")
matpoints(aq$time, aq$Ozone, type="l", lty=1,col="green")
matpoints(aq$time, aq$Solar.R, type="l", lty=1,col="skyblue")

axis.POSIXct(side=1,
             at=seq(as.POSIXct(aq$time[1]), 
             as.POSIXct(aq$time[153]), 
             "5 day"), las=2, format="%y-%m-%d")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/743e4192-1dc7-0c18-c58e-b4a97afbd3e1.png)


ところどころデータが欠損してます。

```R
summary(aq)
```

```
     Ozone           Solar.R           Wind             Temp      
 Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
 1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
 Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
 Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
 3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
 Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
 NA's   :37       NA's   :7                                       
     Month            Day            time           
 Min.   :5.000   Min.   : 1.0   Min.   :1973-05-01  
 1st Qu.:6.000   1st Qu.: 8.0   1st Qu.:1973-06-08  
 Median :7.000   Median :16.0   Median :1973-07-16  
 Mean   :6.993   Mean   :15.8   Mean   :1973-07-16  
 3rd Qu.:8.000   3rd Qu.:23.0   3rd Qu.:1973-08-23  
 Max.   :9.000   Max.   :31.0   Max.   :1973-09-30  
```
OzoneとSolar.Rが欠損値を含んでいるようです。
NAのある行を単純に削除するなら

```R
aq_del_na=na.omit(aq)
```

ですがplotしてみたら波形になっていて、前後の値と関係していそうなので、
前後の値の合計の平均を使って埋めたいと思います。
より正確に埋めたかったらSVMとか使って欠損値の推測とかもするんでしょうか。

```R
library(zoo)
aq$Ozone=na.approx(aq$Ozone)
aq$Solar.R=na.approx(aq$Solar.R)
```

Ozoneの上から五行目を見てみるとNAが23になっています
Ozoneの4行目(18)と6行目(28)を足した平均が23なのでできていそうです。

欠損値を埋めた状態でplotしてみます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/310bd07c-07e4-9d9c-3345-2320f305a642.png)

一見して埋めた値におかしなところはなさそうです。

データの列に月や日が入っているので、単純に変数だけにします。

```R
aq_relation_plot=data.frame(Ozone=aq$Ozone,Solar=aq$Solar.R,Wind=aq$Wind,Temp=aq$Temp)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/4a54a4b1-1ea5-04fc-7bc5-86b359965118.png)

関係を各変数でplotして比較してみたら
・windとtempで負の相関
・windとozoneに微妙に負の相関
・tempとozoneに正の相関
がありそうでしょうか？

```R
cor(aq_relation_plot)
```

```
           Ozone       Solar        Wind       Temp
Ozone  1.0000000  0.15947312 -0.50753419  0.6016359
Solar  0.1594731  1.00000000 -0.02927963  0.2296967
Wind  -0.5075342 -0.02927963  1.00000000 -0.4579879
Temp   0.6016359  0.22969669 -0.45798788  1.0000000
```

実際に相関係数を確認してみます
「ozoneとsolar」「windとsolar」の間に相関はなさそうです

実際に重回帰してモデルにしてみます。

```R
model=lm(Ozone~Solar+Wind+Temp,data=aq_relation_plot)
summary(model)
```

```
Call:
lm(formula = Ozone ~ Solar + Wind + Temp, data = aq_relation_plot)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.651 -15.622  -4.981  12.422 101.411 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -52.16596   21.90933  -2.381   0.0185 *  
Solar         0.01654    0.02272   0.728   0.4678    
Wind         -2.69669    0.63085  -4.275 3.40e-05 ***
Temp          1.53072    0.24115   6.348 2.49e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.26 on 149 degrees of freedom
Multiple R-squared:  0.4321,    Adjusted R-squared:  0.4207 
F-statistic: 37.79 on 3 and 149 DF,  p-value: < 2.2e-16
```

43%くらいデータを説明できていそうなモデルが出来ました。
結果を見てみると、

ozone = solar*0.01654 + wind*(-2.69669) + 1.53072*temp -52.16596

というモデルが出てきました。
検定結果を見るとsolarはほとんど影響していなさそうなので、

ozone = wind*(-2.69669) + 1.53072*temp -52.16596

という式で十分に説明が出来そうです。

#これにて終わり!・・・・・でいいのか?

そこまで重要でない予測ならいいですが、
せっかく多変量なのでモデル選択やモデルの検証も行っていきます。

#分散拡大係数(variance inflation factor : VIF)

説明変数の相関係数行列の対角要素。

パッケージcarのVIFの関数を使うと多重共線性(multicollinearity)の程度が計算できるようです。
10よりも大きいと変数間の相関が大きいという事になり、そのような変数はモデルから取り除くのが良い。

多変量の重回帰ではデータを説明する変数が多いほどモデルの当てはまりがよくなる。
説明変数間に強い相関のあるデータが混ざっていると検定結果が有意になったり、実はモデルがデータを全く説明できていないのに、決定係数が見かけ上だけ大きな値になってしまう
などが起こるようです。

もし大きければ変数を一部外したり、取り除くようです。

```R
library(car)
print(vif(model))
```

```
   Solar     Wind     Temp 
1.063891 1.275246 1.345122 
```

モデルの説明変数に強い相関を持つデータはありませんでした。

#モデルの最適化 AIC(Akaike's Information Criterion)

モデルの最適さを表すと言えばこのAICが出てきます。
AICはモデルを比較して理想を探してくれる関数です。
Rでは比較して一番いい物を見つけてくれるstep関数があるので使いましょう。


```
best_model=step(model)
```

```
Start:  AIC=979.78
Ozone ~ Solar + Wind + Temp

        Df Sum of Sq    RSS     AIC
- Solar  1     311.9  88035  978.32
<none>                87723  979.78
- Wind   1   10758.1  98481  995.48
- Temp   1   23721.2 111444 1014.40

Step:  AIC=978.32
Ozone ~ Wind + Temp

       Df Sum of Sq    RSS     AIC
<none>               88035  978.32
- Wind  1     10520  98555  993.59
- Temp  1     26642 114677 1016.77
```

モデルとしてはsolarを外したものが最もAICを下げてくれていいモデルであるようですね。
さっきの上で述べた検定とも一致します。


```R
summary(best_model)
```

```
Call:
lm(formula = Ozone ~ Wind + Temp, data = aq_relation_plot)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.563 -15.252  -4.606  12.607 102.384 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -52.8013    21.8576  -2.416   0.0169 *  
Wind         -2.6564     0.6274  -4.234 3.99e-05 ***
Temp          1.5734     0.2335   6.738 3.23e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.23 on 150 degrees of freedom
Multiple R-squared:  0.4301,    Adjusted R-squared:  0.4225 
F-statistic:  56.6 on 2 and 150 DF,  p-value: < 2.2e-16
```

モデルの説明具合は43%くらいのままですね。

今回得られた予測式を使って得られる予測値と、実際の値を並べてみます

```R
plot(x=best_model$fitted.value, y=aq$Ozone,
     main="Ozone prediction",
     xlab="predict --- lm(formula = Ozone ~ Wind + Temp)",
     ylab="real data",
     xlim=c(-30, 120), ylim=c(0,150))
abline(0, 1, col="red")
abline(h=0, v=0, col="gray")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/453f55e5-3e43-250c-69b5-cbce72196bc7.png)


予測値は風の負の符号の影響を受けてオゾンがマイナスになってしまっている箇所がありますが、概ね予測式が出来ていそうです。

#以上

今後は
リッジやPLS回帰などもまとめていきたいと思います。