pythonが流行する以前、Rで統計モデルを求めていた頃を思い出して、今日はRで重回帰分析(線形回帰)のプログラムを書いてみようと思います。実践的な流れとして、
を説明していきます。
線形回帰モデルを求める理論的な話は、この過去の記事を参照してください。
(理論編)線形回帰モデルの計算方法 - 数学やプログラミングの備忘録
また、今日のソースコードとデータは、以下にまとめているので、ご自由にお持ち帰りください。
GitHub - cresselia2012/R-samples
目次
CSVを読み込む
線形回帰の対象となる以下のデータを読むことを考えます。
5人の学生の「五教科合計」「数学」「理科」「英語」の点数のデータです(本来ならもっとたくさんの学生のデータが欲しいところですが、例ですので)。
この記事では、以下の線形回帰モデルを推定します。
では、データを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と引けを取らない気がしましたね。
おわり