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')

信頼区間と似た概念として予測区間というものもある。信頼区間は、同じ条件下で再実験を行い、その実験結果を利用して回帰直線を求めたとき、その回帰直線が通る区間のことである。これに対して、予測区間は同じ条件下で再実験を行い、その実験結果（各点）そのものが存在する区間のことである。予測区間も predict 関数で求めることができ、その際 interval = 'prediction' オプションを与える必要がある。

newx <- data.frame(speed = seq(5, 25, 0.1))
conf.interval <- predict(model, newdata = newx, interval = 'prediction', level = 0.95)

plot(cars$speed, cars$dist)
abline(model)
lines(newx$speed, conf.interval[, 1], col = 'orange')
lines(newx$speed, conf.interval[, 2], col = 'blue')
lines(newx$speed, conf.interval[, 3], col = 'blue')

5.3 重回帰分析

説明変数が複数個あるときに、単回帰分析と区別するために重回帰分析とよばれることがある。生物学分野で観測されたデータの中に説明変数が 1 つだけのものが少ないことから、一般に回帰分析と書かれた場合は重回帰分析を指すことがほとんどである。

ここでは trees データセットを利用して、R で重回帰分析を行う方法を示す。trees データセットは 3 列からなり、それぞれ周長（Girth）、高さ（Height）、体積（Volume）となっている。ここでは、周長と高さを利用して、体積を説明するモデルの定式化を試みる。

data(trees)
head(trees)
##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7

データの分布などを確認する。散布図からは、周長（Girth）および高さ（Height）のどれを使用しても体積（Volume）を説明できることが確認取れる。とりわけ、体積を説明するのに周長の方が重要と思われる。

pairs(trees)

周長と高さの単位は同じであるものの、値のスケールが異なるので、まず両者を標準化する。

trees_std <- trees
trees_std$Girth <- scale(trees_std$Girth)
trees_std$Height <- scale(trees_std$Height)
head(trees_std)
##        Girth     Height Volume
## 1 -1.5768542 -0.9416472   10.3
## 2 -1.4812561 -1.7263533   10.3
## 3 -1.4175241 -2.0402357   10.2
## 4 -0.8758017 -0.6277648   16.4
## 5 -0.8120696  0.7847060   18.8
## 6 -0.7802036  1.0985884   19.7

次に、標準化された周長と高さを使って体積を説明する回帰モデルを構築してみる。体積（Volume）を目的変数として、周長（Girth）と高さ（Height）を説明変数として回帰分析を行いたいので、統計的な記述方法では Volume ~ Girth + Height と書ける。この式をそのまま lm 関数に代入して実行する。

model <- lm(Volume ~ Girth + Height, data = trees_std)
summary(model)
## 
## Call:
## lm(formula = Volume ~ Girth + Height, data = trees_std)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4065 -2.6493 -0.2876  2.2003  8.4847 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.1710     0.6972  43.275   <2e-16 ***
## Girth        14.7749     0.8293  17.816   <2e-16 ***
## Height        2.1616     0.8293   2.607   0.0145 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
## F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

回帰直線の係数を確認すると、Girth の係数が 14.7748603 であり、Height の係数が 2.1616454 である。Girth および Height 変数のデータはすでに標準化され、同じスケールである。そのため、これらの係数から、Volume を説明するには、Height に比べ Girth が重要であると考えられる。

5.4 変数選択

回帰分析において、モデルの中に複数の説明変数を入れることが可能だが、変数が多くなればなるほどよいモデルができると限らない。 そこで、回帰分析を行うとき、複数の変数を様々な組み合わせでモデル化し、そのうち最適なモデルを選択するという作業を行う必要がある。統計の分野において、モデルの良さを測る指標として AIC が一般的に使われている。回帰モデルの AIC を計算するには AIC 関数を使う。ここで、周長と高さを使って体積を説明する回帰モデルの AIC を計算してみる。

model <- lm(Volume ~ Girth + Height, data = trees_std)
AIC(model)
## [1] 176.91

次に、周長だけを使って体積を説明するモデルと高さだけを使って体積を説明するモデルを作成して、それぞれの AIC を計算してみる。

model_r1 <- lm(Volume ~ Girth, data = trees_std)
AIC(model_r1)
## [1] 181.6447
model_r2 <- lm(Volume ~ Height, data = trees_std)
AIC(model_r2)
## [1] 252.7986

これまでに構築した 3 つのモデルの AIC を比べると、周長と高さの両方の説明変数を使ってモデルを構築したとき AIC が最も低いことがわかった。つまり、このデータセットに関しては、周長と高さの両方を使った方が最適なモデルが得られるということになる。

モデルの選択は変数選択とも呼ばれている。trees データセットを使用した例では、考えられる説明変数の組み合わせのすべてで回帰モデルを作成して、その中から最適なモデルを決めている。これを総当たり法という。説明変数が 3 つのときは 3 通りの組み合わせ（Girth + Height; Girth; Height）を試せばよいが、 説明変数が多いときその組み合わせ数が爆発的に増える。このようなデータセットに対して、変数増加法やスパース回帰などの方法で変数選択が行われる（参考）。統計分野では、適切な変数選択法を利用して解釈可能な最適なモデルを構築していくことが求められる。

CC BY 4.0