重回帰分析2

# 重回帰分析2

なぜ，重回帰分析についての記事が何回も何回も掲載されるのか。それも，重回帰分析すらでなくて単回帰（直線回帰）分析。しかも，不十分な記事が。<br>

重回帰分析の導出なんて，誰が何回やっても同じなんだから，それを個別に何回も何回も記述しても，意味ないじゃん。

その割には，注意すべき点についてちゃんと記述できている記事も少ないし。

重回帰分析は一筋縄ではいかない。数式だけ眺めていても，わからないことがたくさんある。

プラットフォームはいろいろあるが，R でやってみる（Python の人は書き換えてみてね）。

まずは，データフレームを作る。ただし，ある意図を持ったデータなので，注意深く分析することを心がけて欲しい。


```R
library(MASS)
options(digits=16)
set.seed(1234567890)
nc = 20
r = matrix(c(
        1.0, 0.6, 0.5, 0.1,
        0.6, 1.0, 0.4, 0.7,
        0.5, 0.4, 1.0, 0.6,
        0.1, 0.7, 0.6, 1.0), ncol=4)
z <- mvrnorm(nc, mu=rep(0, 4), r, empirical=TRUE)
z = round(t(t(z)*c(3, 4, 5, 6) + c(50, 40, 30, 60)), 1)
df = data.frame(z)
colnames(df) = c("y", "x1", "x2", "x3")
head(df)
```


<table class="dataframe">
<caption>A data.frame: 6 × 4</caption>
<thead>
	<tr><th></th><th scope=col>y</th><th scope=col>x1</th><th scope=col>x2</th><th scope=col>x3</th></tr>
</thead>
<tbody>
	<tr><th scope=row>1</th><td>45.4</td><td>38.3</td><td>27.7</td><td>63.5</td></tr>
	<tr><th scope=row>2</th><td>50.4</td><td>44.0</td><td>37.1</td><td>66.7</td></tr>
	<tr><th scope=row>3</th><td>48.1</td><td>39.2</td><td>31.8</td><td>61.8</td></tr>
	<tr><th scope=row>4</th><td>47.9</td><td>31.3</td><td>24.1</td><td>46.9</td></tr>
	<tr><th scope=row>5</th><td>49.2</td><td>40.4</td><td>32.4</td><td>67.6</td></tr>
	<tr><th scope=row>6</th><td>49.1</td><td>43.1</td><td>29.7</td><td>65.1</td></tr>
</tbody>
</table>




相関係数を求めてみよう。実は，データフレームを作るときに指定した r に他ならない（データを作るときに丸めたのでちょっとだけ違う）

```R
round(cor(df), 4)
```


<table class="dataframe">
<caption>A matrix: 4 × 4 of type dbl</caption>
<thead>
	<tr><th></th><th scope=col>y</th><th scope=col>x1</th><th scope=col>x2</th><th scope=col>x3</th></tr>
</thead>
<tbody>
	<tr><th scope=row>y</th><td>1.0000</td><td>0.5982</td><td>0.4994</td><td>0.0979</td></tr>
	<tr><th scope=row>x1</th><td>0.5982</td><td>1.0000</td><td>0.4017</td><td>0.7019</td></tr>
	<tr><th scope=row>x2</th><td>0.4994</td><td>0.4017</td><td>1.0000</td><td>0.5991</td></tr>
	<tr><th scope=row>x3</th><td>0.0979</td><td>0.7019</td><td>0.5991</td><td>1.0000</td></tr>
</tbody>
</table>




散布図を描いてみても，相関係数が 0.5 とか 0.7 くらいでも，あまりぱっとしない。
こんなので予測できるか？（機械学習の教育用データは，もっと相関の低い変数が一杯ある）

```R
plot(df)
```


    
![output_3_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/68a75ee4-ac00-b73b-39ca-7fc1b3b12207.png)

    


データについての所見

従属変数 y ともっとも相関の高いのは x1。次いで x2。x3 とはほとんど相関がないといえるレベルだ。

普通の人は，「重回帰分析するときに，x3 なんか使わないよね！」と思うかも知れない。


```R
lm.ans1 = lm(y ~ x1 + x2, data=df)
summary(lm.ans1)
```


    
    Call:
    lm(formula = y ~ x1 + x2, data = df)
    
    Residuals:
                 Min               1Q           Median               3Q 
    -3.5656990721298 -1.8009686436998 -0.9731367460045  2.2515107463428 
                 Max 
     3.1338362693141 
    
    Coefficients:
                        Estimate       Std. Error t value   Pr(>|t|)    
    (Intercept) 30.1508994389561  5.6250775499563 5.36009 5.1956e-05 ***
    x1           0.3566805720560  0.1493963265573 2.38748   0.028847 *  
    x2           0.1860625892935  0.1195547419975 1.55630   0.138057    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 2.387951045993 on 17 degrees of freedom
    Multiple R-squared:  0.437902856506,	Adjusted R-squared:  0.3717737808009 
    F-statistic:  6.62194128432 on 2 and 17 DF,  p-value: 0.007471309290127


```R
plot(df$y, predict(lm.ans1))
abline(0, 1, col="red")
```

    
![output_5_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/7b9f453f-b46e-c0af-f7de-60e483563b62.png)

    




どうですか，分析結果の解釈は。

決定係数 0.4379 とな。これじゃあ，重回帰分析は失敗だ。実測値と予測値の図を見てもわかる。

回帰の分散分析の p 値が 0.0074 であってもだ。

y との相関係数が 0.5 くらいもあった x2 すら，偏回帰係数の検定の p 値は 0.1381 で，偏回帰係数が 0 であるという帰無仮説を棄却できていない。

ところで，x3 も使ったらどうなると思う？（そもそも，使おうなんて思っていないだろうが）。



```R
lm.ans2 = lm(y ~ x1 + x2 + x3, data=df)
summary(lm.ans2)
```


    
    Call:
    lm(formula = y ~ x1 + x2 + x3, data = df)
    
    Residuals:
                 Min               1Q           Median               3Q 
    -2.2483635860615 -0.4904258216178 -0.2910102631902  0.4773558776720 
                 Max 
     1.9745355334776 
    
    Coefficients:
                         Estimate        Std. Error  t value   Pr(>|t|)    
    (Intercept) 37.45522309865842  2.67953807205067 13.97824 2.1903e-10 ***
    x1           0.80518187920135  0.08640552755012  9.31864 7.2650e-08 ***
    x2           0.43287512241041  0.06151436784887  7.03698 2.8081e-06 ***
    x3          -0.54414586565041  0.06594139634562 -8.25196 3.7001e-07 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 1.073655938906 on 16 degrees of freedom
    Multiple R-squared:  0.8930546607771,	Adjusted R-squared:  0.8730024096729 
    F-statistic: 44.53637928892 on 3 and 16 DF,  p-value: 5.434308586e-08



なんということだ。x1, x2, x3 ともに，偏回帰係数は有意ということになった。それもすさまじい p 値で。

決定係数は 0.9 近いし，自由度調整済みの決定係数だって 0.87 だ。

当然，回帰の分散分析だって，有意だ。

実測値と予測値の図を見るまでもない。


```R
plot(df$y, predict(lm.ans2))
abline(0, 1, col="red")
```


    
![output_9_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/d4313725-e606-797b-71cd-468345e678fc.png)

    




ところで，重回帰分析にはまだまだ深い闇（なぞ）がある。

上の重回帰分析の結果で，x3 の偏回帰係数は -0.5355 と，負だよね。y と x3 の相関係数は 0.0979 だったのにね。

このなぞをちゃんと説明できないと，先には進めないよ。