数学やプログラミングの備忘録

数理最適化, Python, C++をメインに紹介するブログ。

MENU

実践的!Rによる重回帰分析

pythonが流行する以前、Rで統計モデルを求めていた頃を思い出して、今日はRで重回帰分析(線形回帰)のプログラムを書いてみようと思います。実践的な流れとして、

  1. CSV形式のデータを読み込む
  2. 読み込んだデータから必要な部分を抽出する
  3. 重回帰モデルを求める&結果を見る
  4. AICの計算やステップワイズ法の適用

を説明していきます。

線形回帰モデルを求める理論的な話は、この過去の記事を参照してください。

(理論編)線形回帰モデルの計算方法 - 数学やプログラミングの備忘録

また、今日のソースコードとデータは、以下にまとめているので、ご自由にお持ち帰りください。

GitHub - cresselia2012/R-samples

目次

CSVを読み込む

線形回帰の対象となる以下のデータを読むことを考えます。

5人の学生の「五教科合計」「数学」「理科」「英語」の点数のデータです(本来ならもっとたくさんの学生のデータが欲しいところですが、例ですので)。

この記事では、以下の線形回帰モデルを推定します。


(五教科合計) = b + a_1 \times (数学) + a_2 \times (理科) + a_3 \times (英語)

では、データをRで読み込んでみましょう。

> data <- read.csv("/path/data.csv")
  • > より後ろをRコンソールに打ち込みます
  • /path の部分はCSVファイルが格納されている場所に応じて、書き直してください。

ちゃんと読み込まれたか確認してみましょう。Rコンソールに変数 data を打ち込めば、中身が出力されます。

> data
五教科合計 数学 理科 英語
1        450  100   95   80
2        400   80   90   70
3        350   80   85   60
4        300   50   55   60
5        250   30   50   40

データを抽出

(今回は読み込んだデータ全てを使いますが)次に、読み込んだデータの必要な部分を抽出します。まずは、目的変数となる 五教科合計 の値をベクトル y に格納します。

> y <- data[,1]

Rは配列を1から数えるんですねえ、、笑

次に、残りの 数学 理科 英語 を5×3行列に格納します。

# 空の5×3行列を生成
> X <- matrix(NA, 5, 3)
# 2列目, 3列目, 4列目を抽出
> for (i in 1:length(data)-1) {
+     X[,i] <- data[,i+1]
+ }

ここで、一度 y, X の中身を確認してみます。

> y
[1] 450 400 350 300 250
> X
     [,1] [,2] [,3]
[1,]  100   95   80
[2,]   80   90   70
[3,]   80   85   60
[4,]   50   55   60
[5,]   30   50   40

思い通り格納できてますね。

重回帰モデルを求める

$$ (五教科合計) = b + a_1 \times (数学) + a_2 \times (理科) + a_3 \times (英語) $$

いざ、上記の重回帰モデルを求めてみましょう!

> result <- lm(y~X)

たったこれだけです。結果のサマリーを見るときは、summary() 関数を使います。

> summary(result)

Call:
lm(formula = y ~ X)

Residuals:
     1      2      3      4      5
11.981 -4.792 -8.786 -7.987  9.585

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    5.591    103.689   0.054    0.966
X1            -0.639      2.363  -0.270    0.832
X2             2.332      2.259   1.032    0.490
X3             3.434      2.044   1.680    0.342

Residual standard error: 19.98 on 1 degrees of freedom
Multiple R-squared:  0.984,  Adjusted R-squared:  0.9361
F-statistic: 20.53 on 3 and 1 DF,  p-value: 0.1605

結果の見方

  • Residuals: データとの誤差
  • Coefficients: モデルの係数
    • Intercept: モデルの切片 $b$
    • Xk: モデルの係数 $a_k$

AICの計算やステップワイズ法の適用

他にも、求めたモデルのAIC(赤池情報量規準)を計算することができます。

> AIC(result)
[1] 46.09152

より予測精度の高いモデルを見つけたいときは、ステップワイズ法を適用します。

> step(result)
Start:  AIC=29.9
y ~ X

       Df Sum of Sq     RSS    AIC
<none>                399.4 29.902
- X     3     24601 25000.0 44.586

Call:
lm(formula = y ~ X)

Coefficients:
(Intercept)           X1           X2           X3  
      5.591       -0.639        2.332        3.435  

まとめ(ソースコード)

# read
data <- read.csv("/path/data.csv")

# generate y and X
y <- data[,1]
X <- matrix(NA, 5, 3)
for (i in 1:length(data)-1) {
    X[,i] <- data[,i+1]
}

# linear regression
# y = a1 x1 + a2 x2 + a3 x3 + b
result <- lm(y~X)

# summary of result
summary(result)

# AIC
AIC(result)

# stepwise method
step(result)

久しぶりにRで書いてみましたが、重回帰分析に限っていえば、使い勝手に関してpythonと引けを取らない気がしましたね。

おわり