データの用意

このページでは、乱数を使って仮想的なデータを作成します。正規分布に従う乱数は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関数を使えます。

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

summary関数で回帰の結果の概要を知ることができます。

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ではそれを自動化して行うことができます。

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

最終的にx4が削られているのがわかります。