5 回帰分析

R には統計解析や統計検定などを行う機能が多く実装されている。この資料では、R の使い方を紹介する目的で、lm 関数を使用した回帰分析を中心に紹介する。回帰分析、変数選択や一般化線形モデルなどを詳しく勉強されたい方は、次の資料を参考にしてください。

5.1 単回帰

1 つの説明変数で 1 つの目的変数を説明するような回帰分析は、単回帰とも呼ばれている。R で単回帰分析を行うには lm 関数を使う。ここでは、R の標準パッケージに付属している cars データセットを使用して、回帰分析を行う例を示す。まず、cars データセットの特徴を確認する。

data(cars)
head(cars)
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10

cars データセットは 2 列のデータフレームからなっている。speed 列は車の走行速度（mph）であり、dist 列はその走行速度からブレーキをかけたときに車が止まるまでに必要な距離（フィート）である。両者の関係を図示して確認してみると、走行速度が高ければ、停止に必要な距離も長くなることがわかる。

plot(cars$speed, cars$dist)

このデータセットを使って回帰分析（単回帰）を行い、走行速度と停止距離の関係を定式化する。回帰分析は lm 関数にモデル式を代入して行う。統計分野において、説明変数 speed で目的変数 dist をモデル化するときに、一般的に dist ~ speed のように記述する。R の lm 関数はこの記述法に対応しているので、そのまま lm 関数にこのモデル式を代入すれば回帰分析が行われる。

model <- lm(dist ~ speed, data = cars)

回帰直線の係数やその係数に対する統計検定結果などは、summary 関数の出力内容から確認できる。この例では、dist = 3.9324088 * speed - 17.5790949 という回帰直線の式が求められたことを意味している。この回帰直線のことをモデルとも呼ぶ。

summary(model)
## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

summary 関数の出力情報が多いため、ここでは回帰直線（モデル）の係数だけを出力させてみることにする。

coef(model)
## (Intercept)       speed 
##  -17.579095    3.932409

次に、回帰直線（モデル）を実際のデータ上に重ねて図示してみる。この回帰直線（モデル）が走行速度と停止距離の関係を概ねに説明できていることわかる。

plot(cars$speed, cars$dist)
abline(a = coef(model)[[1]], b = coef(model)[[2]])

回帰直線（モデル）の適合値（予測値）と実際の値の残差を図示した残差プロットを図示する。まず残差（実測値と適合値）を計算する。ここで計算される残差は、summary 関数の出力結果のうち Residuals 項の情報と同じものである。なお、適合値を計算するには、実際の走行速度（speed）の値を predict 関数に代入して求める。

newx <- data.frame(speed = cars$speed)
dist_hat <- predict(model, newdata = newx)

residual_error <- cars$dist - dist_hat
summary(residual_error)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -29.069  -9.525  -2.272   0.000   9.215  43.201

計算された残差で残差プロットを描く。この残差プロットを見ると、走行速度が遅いときは予測がやや正確で、走行速度が速くなると予測がやや不正確になる傾向が確認取れる。

plot(residual_error)
abline(h = 0)

5.2 信頼区間と予測区間

cars データセットを使って求めた回帰直線（モデル）は、cars データセットに基づくものとなっている。cars データセットの測定時と同じ条件でデータを取り直したときに、異なる回帰直線が求められる、と考えられる。このような実験を無数回繰り返していくと、無数の回帰直線が描かれる。走行速度と停止距離の間に真に比例関係が見られるならば、これらの無数の回帰直線はほぼ同じ傾きや切片を持つと考えられる。そこで、無数回の実験により無数本の回帰直線を描いたときに、そのうち 95% の回帰直線が通る領域について見積もってみる。このような領域のことを 95% 信頼区間とよぶ。

信頼区間を求めるには、各係数に対して区間推定を行う必要がある（参考）。しかし、R の predict 関数に interval = 'confidence' オプションを与えると、自分で計算することなく容易に 95% 信頼区間を求めることができる。R で容易に求められるものの、その仕組みを解析者が理解していくことが重要である。信頼区間の求め方を一度確認しておくとよい。なお、predict 関数の出力結果が 3 列になり、1 列目が適合値（予測値）であり、2 列目および 3 列目は 95% 信頼区間の下限および上限となっている。

newx <- data.frame(speed = seq(5, 25, 0.1))
conf.interval <- predict(model, newdata = newx, interval = 'confidence', level = 0.95)
head(conf.interval)
##        fit       lwr      upr
## 1 2.082949 -7.644150 11.81005
## 2 2.476190 -7.176357 12.12874
## 3 2.869431 -6.708713 12.44757
## 4 3.262672 -6.241220 12.76656
## 5 3.655912 -5.773883 13.08571
## 6 4.049153 -5.306705 13.40501
plot(cars$speed, cars$dist)
abline(model)
lines(newx$speed, conf.interval[, 1], col = 'orange')
lines(newx$speed, conf.interval[, 2], col = 'darkgreen')
lines(newx$speed, conf.interval[, 3], col = 'darkgreen')