データの用意
このページでは、乱数を使って仮想的なデータを作成します。正規分布に従う乱数はrnorm関数で生成できます。rnorm関数は第一引数に生成する個数、第二引数に平均、第三引数に分散を指定します。
x1 <- rnorm(100, 10, 1)
x2 <- rnorm(100, 10, 2)
x3 <- rnorm(100, 10, 3)
y <- x1 + 2 * x2 + 3 * x3 + rnorm(100, 0, 1)
lm関数で重回帰分析をする
説明変数が1つの時もlm関数でしたが、説明変数が増えても説明変数と従属変数の関係が線形結合であればlm関数を使えます。これも、summary関数で回帰の結果の概要を知ることができます。
result <- lm(y ~ x1 + x2 + x3)
summary(result)
以下のように出力されます。
Call:
lm(formula = y ~ x1 + x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-1.96202 -0.62604 -0.01741 0.73180 2.64501
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.80572 1.21700 0.662 0.51
x1 0.85806 0.09992 8.588 1.6e-13 ***
x2 2.03232 0.04711 43.140 < 2e-16 ***
x3 3.04166 0.03769 80.698 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9662 on 96 degrees of freedom
Multiple R-squared: 0.9887, Adjusted R-squared: 0.9883
F-statistic: 2798 on 3 and 96 DF, p-value: < 2.2e-16
y = 0.80572 * x1 + 2.03232 * x2 + 3.04166 * x3 + 0.80572 の式が得られています。乱数を用いたデータのため、同じ結果が得られるとは限りません。今回の例では各係数の有意確率が非常に小さく、係数は0.1%有意水準を満たしているとわかります。Interceptは切片のことで、有意確率が小さくありませんが、これは切片が0と差がないことを否定しにくいだけで通常無視されます。その下に続くもののうち重要なのが、Multiple R-squared: 0.9887 と Adjusted R-squared: 0.9883 のところで、それぞれ決定係数、自由度調節済み決定係数と言い、被説明変数(y)の変化のうち、推定式で説明できた割合のことを指します。今回は説明変数が複数であるため自由度調節済み決定係数を見るべきです。重回帰分析などで説明変数が増えるときには、意味のない説明変数を含めても決定係数は上昇します。この欠点を補うために自由度調整済み決定係数を用いることが重要です。
係数の有意性
意味のない説明変数が含まれている場合には、summary関数の出力で見分けることができます。まず、x4としてyと無関係な乱数を生成し、もう一度lm関数を使ってx1~x4を説明変数とした重回帰分析をします。以下のように出力されます。
x4 <- rnorm(100, 10, 4)
result2 <- lm(y ~ x1 + x2 + x3 + x4)
summary(result2)
以下のように出力されます。
Call:
lm(formula = y ~ x1 + x2 + x3 + x4)
Residuals:
Min 1Q Median 3Q Max
-1.95356 -0.61666 -0.03673 0.74286 2.65140
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.780315 1.227547 0.636 0.527
x1 0.856324 0.100668 8.506 2.55e-13 ***
x2 2.031220 0.047563 42.706 < 2e-16 ***
x3 3.041296 0.037908 80.228 < 2e-16 ***
x4 0.005898 0.024471 0.241 0.810
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.971 on 95 degrees of freedom
Multiple R-squared: 0.9887, Adjusted R-squared: 0.9882
F-statistic: 2078 on 4 and 95 DF, p-value: < 2.2e-16
x4の係数が有意でないことがわかります。すなわち、x4の係数が0でないということを否定できません。
説明変数の選択
説明変数の中で無関係なものが存在する可能性がある場合には、ステップワイズ法により適切と思われる変数を逐一選んでいく方法が有効な場合があり、Rではそれを自動化して行うことができます。以下のように出力されます。最終的にx4が削られているのがわかります。
result3 <- lm(y~1)
step(result3, direction="both", scope=list(upper=~x1+x2+x3+x4))
Start: AIC=417.69
y ~ 1
Df Sum of Sq RSS AIC
+ x3 1 9548.9 1946.6 300.87
+ x2 1 1748.3 9747.2 461.96
11495.5 476.45
+ x4 1 2.8 11492.7 478.43
+ x1 1 0.8 11494.7 478.45
Step: AIC=300.87
y ~ x3
Df Sum of Sq RSS AIC
+ x2 1 1755.6 191.0 70.73
+ x1 1 113.9 1832.7 296.84
1946.6 300.87
+ x4 1 10.2 1936.4 302.34
- x3 1 9548.9 11495.5 476.45
Step: AIC=70.73
y ~ x3 + x2
Df Sum of Sq RSS AIC
+ x1 1 125.6 65.4 -34.43
+ x4 1 3.9 187.1 70.66
191.0 70.73
- x2 1 1755.6 1946.6 300.87
- x3 1 9556.2 9747.2 461.96
Step: AIC=-34.43
y ~ x3 + x2 + x1
Df Sum of Sq RSS AIC
65.4 -34.43
+ x4 1 0.5 65.0 -33.15
- x1 1 125.6 191.0 70.73
- x2 1 1767.2 1832.7 296.84
- x3 1 9679.7 9745.1 463.94
Call:
lm(formula = y ~ x3 + x2 + x1)
Coefficients:
(Intercept) x3 x2 x1
0.8057 3.0416 2.0323 0.8580
多重共線性とは
重回帰分析において、説明変数をやみくもに増やすと多重共線性という問題が生じることがあります。これは、説明変数間に相関がある場合に起こり、多重共線性が生じると回帰式の係数が不安定になります。重回帰分析を行うときには、これに気を配る必要があります。
多重共線性の発生する例
まず、乱数を使って仮想的なデータを生成します。以下の例では、もちろんx1?x3とyの間に線形の関係がありますが、さらにx1とx3の間にも相関があるデータを作成していることに注意してください。
x1 <- rnorm(100, 10, 1)
x2 <- rnorm(100, 10, 2)
x3 <- x1 * 2 + rnorm(100, 10, 0.01)
y <- x1 + 2 * x2 + 3 * x3 + rnorm(100, 0, 1)
ここで散布図行列を書きます。x1とx3の間に相関があるのがわかります。
data <- data.frame(y,x1,x2,x3)
pairs(data)
lm関数で回帰を行い、summary関数で回帰の結果の概要を見ます。以下のように出力されます。
result <- lm(formula = y ~ x1 + x2 + x3)
summary(result)
Call:
lm(formula = y ~ x1 + x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-2.33112 -0.82727 -0.06202 0.88366 2.30364
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -124.64723 119.91432 -1.039 0.301
x1 -23.70632 23.94031 -0.990 0.325
x2 1.96392 0.05202 37.757 <2e-16 ***
x3 15.40270 11.97863 1.286 0.202
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.07 on 96 degrees of freedom
Multiple R-squared: 0.9842, Adjusted R-squared: 0.9837
F-statistic: 1993 on 3 and 96 DF, p-value: < 2.2e-16
仮想的なデータを作った方法によれば、回帰式は y = x1 + 2 * x2 + 3 * x3 に近いものとなるはずですが、全く異なる係数が得られてしまっています。これが多重共線性の影響です。各係数の有意確率も小さくなく、回帰がうまくいっていないことが示唆されます。
VIF値による多重共線性の判別
説明変数の中で多重共線性を生じている可能性がある場合には、VIF(Variance Inflation factor:分散拡大係数)を計算します。10を越えると多重共線性の目安になります。まず、デフォルトのRではVIFは計算できませんので、パッケージをインストールします。
install.packages("fmsb")
library(fmsb)
インストールしたら、以下のようにコマンドを入力すればVIFが計算できます。以下のような出力が得られます。
VIF(lm(x1~x2+x3))
VIF(lm(x2~x1+x3))
VIF(lm(x3~x1+x2))
[1] 44614.76
[1] 1.03132
[1] 44621.33
