R による重回帰分析

重回帰分析は多変量解析法の一つである。一つの従属変数 y を、複数の説明変数 xi でモデル化し、それらの関係や傾向を分析する方法である。モデルを数式で表すと以下のようになる。

\[ y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + \cdots + \beta_{n}x_{n} \]

ここでは R を利用て重回帰分析を行う例を示す。分析結果がわかりやすいように、シミュレーションデータを利用する。シミュレーションデータは以下のモデルとなるように生成する。説明変数を 3 個とする。

\[ y = 10 + 20x_{1} + 30x_{2} + 40x_{3} \]

シミュレーションデータを生成して、従属変数を y、説明変数を x1x2x3 の 3 つのベクトルにそれぞれ保存する。

x1 <- rnorm(100, 5, 1)
x2 <- rnorm(100, 8, 2)
x3 <- rnorm(100, 13, 3)
y <- 10 + 20 * x1 + 30 * x2 + 40 * x3 + rnorm(100, 0, 1)

このデータを利用して回帰分析を行う。回帰分析を行う関数は lmglm とがある。ここでは lm を用いる。まず、ここではモデルを「y = β0 + β1x1 + β2x2 + β3x3」として解析する。

res <- lm(y ~ x1 + x2 + x3)

解析結果を見る際に summary 関数を利用する。

summary(res)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max
## -2.5127 -0.6486 -0.0172  0.8330  2.3367
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.12927    0.90731   11.16   <2e-16 ***
## x1          19.95131    0.10466  190.63   <2e-16 ***
## x2          29.98828    0.05163  580.88   <2e-16 ***
## x3          40.00809    0.03555 1125.51   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 1.023 on 96 degrees of freedom
## Multiple R-squared:  0.9999,	Adjusted R-squared:  0.9999
## F-statistic: 5.177e+05 on 3 and 96 DF,  p-value: < 2.2e-16

Coefficients の項目の Estimate の列が計算された β の値である。(Intercept) が β0 であり、x1、x2、x3 の Estimate の値はそれぞれ β1、β2、β3 の値である。この解析結果から求められる回帰モデルの式は以下のようになり、シミュレーションデータの生成式と非常に似ていることがわかる。

\[ y = 10.12927 + 19.95131 x_{1} + 29.98828 x_{2} + 40.00809 x_{3} \]

また、Pr(>|t|) の列は、「β = 0」を検定したときに得られる p-value を示している。この p-value がゼロに近づけば、「β ≠ 0」として、従属変数 y に影響を与えていると判断することができる。

重要でない説明変数がある場合

仮に、従属変数が普遍で、説明変数が 4 つになった場合はどうなるのかを検証してみる。3 つまでの独立変数はこれまで同様なものを利用し、4 つ目の独立変数を乱数で生成する。4 つ目の独立変数が他の従属変数と独立変数とはまったく関係を持たないことに注意してください。

x4 <- rnorm(100, 21, 4)

では、lm 関数で回帰分析を行ってみる。

res2 <- lm(y ~ x1 + x2 + x3 + x4)
summary(res2)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max
## -2.58947 -0.64411 -0.02052  0.80106  2.35561
## 
## Coefficients:
##             Estimate Std. Error  t value Pr(>|t|)
## (Intercept)  9.92156    1.06926    9.279 5.75e-15 ***
## x1          19.95096    0.10514  189.762  < 2e-16 ***
## x2          29.98797    0.05187  578.189  < 2e-16 ***
## x3          40.00853    0.03573 1119.833  < 2e-16 ***
## x4           0.01021    0.02748    0.371    0.711
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 1.028 on 95 degrees of freedom
## Multiple R-squared:  0.9999,	Adjusted R-squared:  0.9999
## F-statistic: 3.848e+05 on 4 and 95 DF,  p-value: < 2.2e-16

Pr の列を見ると、x4 の値が 0.711 となっている。これは β4 = 0 の検定の結果である。この結果から、β4 はゼロであるということを示している。すなわち、x4 がほとんど y に寄与していないことがわかる。

重回帰分析と標準化

独立変数と従属変数の間に単位の違い等が見られる場合は、両者それぞれ平均 0、分散 1 のデータに標準化する必要がある。単位尺度の変更によって、回帰係数に影響を与えないようにするためには標準化が欠かせない。

上述の yx1x2x3 のデータを利用して説明する。例えば、これまで通りの解析を行うと、回帰係数は以下のように求まる。

res3 <- lm(y ~ x1 + x2 + x3)
res3$coefficients
## (Intercept)          x1          x2          x3
##   10.12927    19.95131    29.98828    40.00809

ここで、例えば x1 はセンチメートルで測っていたことが判定し、これをメートル単位に修正する必要がでてきたとする。すると、x1 全体には 0.01 をかける必要がある。メートルに変換した後に回帰係数を求めるてみると以下のようになる。

x2m <- x2 * 0.01
res4 <- lm(y ~ x1 + x2m + x3)
res4$coefficients
res4$coefficients
## (Intercept)          x1         x2m          x3
##    10.12927    19.95131  2998.82783    40.00809

このように単位の尺度を変更しただけで、x2 の係数が 100 倍されていることがわかる。これでは、単位を変更する度に、回帰係数が変化するためモデル化には不都合である。そこで、各変数をすべて平均 0 分散 1 のデータに転換(標準化)する作業が必要になってくる。

R の scale 関数を利用することで標準化を行うことができる。

y_std <- scale(y)
x1_std <- scale(x1)
x2_std <- scale(x2)
x3_std <- scale(x3)

正規化後のデータで再解析する。

res5 <- lm(y_std ~ x1_std + x2_std + x3_std)
res5$coefficients
##  (Intercept)       x1_std       x2_std       x3_std
## 7.667968e-16 1.595458e-01 4.807529e-01 9.146514e-01

x2m_std <- x2_std * 0.01
res6 <- lm(y_std ~ x1_std + x2m_std + x3_std)
res6$coefficients
##  (Intercept)       x1_std      x2m_std       x3_std
## 7.501174e-16 1.595458e-01 4.807529e+01 9.146514e-01

このように標準化後のデータを回帰分析に用いることで、測定単位を変更しても回帰係数に影響を与えなくなる。