R で多重回帰分析


[昨日](http://qiita.com/ynakayama/items/36b7c1640e6a02ce2e00)は多重回帰分析について書きました。数学的説明については昨日の説明が詳しいのでそちらを見ていただくとして、今日はその補足的な内容として不自由ソフトウェアである Excel の代わりに R を使ってみたいと思います。

R による回帰分析の方法は[ここ](http://cse.naro.affrc.go.jp/takezawa/r-tips/r/71.html)が詳しいです。やりかたについてはリンク先を見ていただくとして、昨日使ったデータを実際に利用してみましょう。データは[こちら](https://github.com/ynakayama/sandbox/blob/master/r/stats/macbook_price.tsv)に用意しました。

## R による回帰分析

```r
# データの読み込み
df <- read.csv("macbook_price.csv", sep="\t")

# 回帰
result <- lm(df$PRICE ~ df$DISPLAY + df$RETINA + df$CPU + df$CORE + df$MEM + df$SSD)

# 回帰分析の要約を表示
summary(result)
#=>
# Call:
# lm(formula = df$PRICE ~ df$DISPLAY + df$RETINA + df$CPU + df$CORE +
#     df$MEM + df$SSD)

# Residuals:
#     Min      1Q  Median      3Q     Max
# -8525.4 -6416.6  -413.8  4082.0 13507.4

# Coefficients: 
#               Estimate Std. Error t value Pr(>|t|)
# (Intercept) -1.603e+05  6.784e+04  -2.362  0.04578 *
# df$DISPLAY   4.783e+03  5.187e+03   0.922  0.38349
# df$RETINA   -6.295e+04  2.883e+04  -2.184  0.06052 .
# df$CPU       7.973e+01  2.708e+01   2.944  0.01859 *
# df$CORE      3.812e+04  1.136e+04   3.357  0.00998 **
# df$MEM      -5.335e-01  1.676e+00  -0.318  0.75839
# df$SSD       5.744e+01  2.229e+01   2.577  0.03276 *
# ---  
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Residual standard error: 9582 on 8 degrees of freedom
# Multiple R-squared:  0.9712,    Adjusted R-squared:  0.9496
# F-statistic: 44.99 on 6 and 8 DF,  p-value: 9.829e-06
```

R2 が 0.9712 、補正 R2 が 0.94 となりました。
切片を強制的に 0 にしていませんので、昨日の Excel の結果とは異なります。やはり昨日と同様、変数の選択があまり良くないみたいですね。

### 変数の選択

昨日やったように変数を減らしてみると次のようになります。

```r
result = lm(df$PRICE ~ df$CPU + df$CORE + df$MEM + df$SSD)
summary(result)
#=>
# Call:
# lm(formula = df$PRICE ~ df$CPU + df$CORE + df$MEM + df$SSD)

# Residuals:
#    Min     1Q Median     3Q    Max 
# -13059  -7035  -1122   1931  23699 

# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)   
# (Intercept) -19840.669  19092.861  -1.039  0.32320   
# df$CPU          27.307      7.697   3.548  0.00529 **
# df$CORE      28389.394   6898.341   4.115  0.00209 **
# df$MEM           1.206      1.656   0.728  0.48315   
# df$SSD          66.256     25.037   2.646  0.02447 * 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Residual standard error: 11070 on 10 degrees of freedom
# Multiple R-squared:  0.9519,    Adjusted R-squared:  0.9327 
# F-statistic: 49.52 on 4 and 10 DF,  p-value: 1.477e-06
```

良くなりました。

### 適切なモデルの探索

R の便利なところは、階層を保ちながら項を増やしたり減らしたりして適切なモデルを探索する関数が簡単に使えるところです。 step() 関数を使うと AIC (赤池情報量基準) を参照してふさわしいモデルが返却されます。

```r
step(result)
#=>
# Start:  AIC=283.29
# df$PRICE ~ df$CPU + df$CORE + df$MEM + df$SSD

#           Df  Sum of Sq        RSS    AIC
# - df$MEM   1   65042361 1291341482 282.06
# <none>                  1226299122 283.29
# - df$SSD   1  858753814 2085052935 289.25
# - df$CPU   1 1543683270 2769982392 293.51
# - df$CORE  1 2076917909 3303217031 296.15

# Step:  AIC=282.06
# df$PRICE ~ df$CPU + df$CORE + df$SSD

#           Df  Sum of Sq        RSS    AIC
# <none>                  1.2913e+09 282.06
# - df$SSD   1 1.0326e+09 2.3239e+09 288.88
# - df$CPU   1 3.2852e+09 4.5765e+09 299.04
# - df$CORE  1 1.1890e+10 1.3181e+10 314.91

# Call:
# lm(formula = df$PRICE ~ df$CPU + df$CORE + df$SSD)

# Coefficients:
# (Intercept)       df$CPU      df$CORE       df$SSD  
#   -29919.53        30.85     32789.22        70.58  
```

これは便利です。

### モデルの比較

それではそれぞれの回帰式で元の価格に対しどれだけ当てはまりがあるのか可視化してみましょう。

summary() 、 step() 、また昨日の Excel での回帰式から算出した値段と実際の値段を表形式に[こちら](https://github.com/ynakayama/sandbox/blob/master/r/stats/macbook_price_check.tsv)のデータにまとめました。

```py3
df = pd.read_csv('macbook_price_check.tsv', index_col=0, sep="\t")
df.index = np.arange(0, len(df.index))

plt.figure()
df.plot()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/38501/994a8a2d-a2e9-91d1-bc41-100c38ed2317.png)

補正 R2 が高いのでどれも当てはまっている感じがしますね。

## まとめ

多重回帰分析の理論と実践について触れました。