多重共線性の回避策
独立変数 $x_4$ が従属変数 $y$ に与える影響を知りたい。ほかの独立変数 $x_1$ と $x_2$ の間に多重共線性があるが，これはどちらも説明変数に含めたい。$x_1$ と $x_4$，$x_2$ と $x_4$ の相関は低く VIF も小さい。
このような場合にどうしたらよいかという質問に対して，シミュレーションに基づく解答を提示する。
以下では R の場合におけるプログラム例であるが，他の言語，パッケージでも追試可能である。
変数間の相関係数行列は以下のような表のようになっていると仮定しよう。
r <- matrix(c(
1, 0, 0, 0, 0,
0.8, 1, 0, 0, 0,
0.3, 0.3, 1, 0, 0,
0.2, 0.2, 0.3, 1, 0,
0.3, 0.25, 0.3, 0.4, 1), byrow=TRUE, ncol=5)
r <- r+t(r)
diag(r) <- 1
names <- c(paste("$x_", 1:4, "$", sep=""), "$y$")
colnames(r) <- rownames(r) <- names
r
| $x_1$ | $x_2$ | $x_3$ | $x_4$ | $y$ | |
|---|---|---|---|---|---|
| $x_1$ | 1.0 | 0.80 | 0.3 | 0.2 | 0.30 |
| $x_2$ | 0.8 | 1.00 | 0.3 | 0.2 | 0.25 |
| $x_3$ | 0.3 | 0.30 | 1.0 | 0.3 | 0.30 |
| $x_4$ | 0.2 | 0.20 | 0.3 | 1.0 | 0.40 |
| $y$ | 0.3 | 0.25 | 0.3 | 0.4 | 1.00 |
まず，表のような相関係数行列になるように，多変量正規分布にしたがうシミュレーションデータを作成する（$n=500$）。
先頭の 10 行は表\ref{tab7-2} のようになる（小数点以下 4 桁以降は表示していない）。
# install.packages("MASS")
library(MASS)
set.seed(12345)
d <- mvrnorm(500, mu=numeric(5), Sigma=r, empirical=TRUE)
d <- data.frame(t(t(d)*c(10, 9, 6, 8, 2)+c(50, 46, 76, 82, 36)))
colnames(d) <- names
head(d, 10)
| $x_1$ | $x_2$ | $x_3$ | $x_4$ | $y$ | |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 35.87129 | 36.70771 | 73.36168 | 66.53958 | 33.12320 |
| 2 | 37.62556 | 29.15146 | 82.96882 | 76.82022 | 33.95505 |
| 3 | 59.21593 | 55.50001 | 77.28258 | 86.95731 | 36.32317 |
| 4 | 63.24307 | 48.87647 | 84.25619 | 81.25347 | 36.78013 |
| 5 | 48.79484 | 39.84562 | 64.62895 | 98.00760 | 39.80084 |
| 6 | 48.31151 | 51.32432 | 89.10405 | 99.74816 | 39.74805 |
| 7 | 44.33141 | 31.49594 | 74.81446 | 86.10918 | 37.23745 |
| 8 | 68.21829 | 64.98163 | 76.57823 | 86.12623 | 36.14495 |
| 9 | 40.22833 | 44.73737 | 81.74600 | 71.22765 | 33.55838 |
| 10 | 57.57091 | 60.43162 | 80.29413 | 80.86618 | 35.98782 |
それぞれの基礎統計量は以下の通りである。
base <- data.frame(平均値=apply(d, 2, mean), 不偏分散=apply(d, 2, var), 標準偏差=apply(d, 2, sd))
rownames(base) <- names
base
| 平均値 | 不偏分散 | 標準偏差 | |
|---|---|---|---|
| <dbl> | <dbl> | <dbl> | |
| $x_1$ | 50 | 100 | 10 |
| $x_2$ | 46 | 81 | 9 |
| $x_3$ | 76 | 36 | 6 |
| $x_4$ | 82 | 64 | 8 |
| $y$ | 36 | 4 | 2 |
散布図は図のようになり，相関係数行列は正確に表の通りになる。
d3 <- d
colnames(d3) <- c(paste("x", 1:4, sep=""), "y")
plot(d3, cex=0.7, pch=16, col="#00aa0060")
partial.cor <- function(x) # データ行列
{
x <- subset(x, complete.cases(x)) # 欠損値を持つケースを除く
i <- solve(cor(x)) # 相関係数行列の逆行列
d <- diag(i) # 対角成分
i <- -i/sqrt(outer(d, d)) # 偏相関係数行列
diag(i) <- NA # 対角成分は未定義
rownames(i) <- colnames(i) <- paste("Var", 1:ncol(x))
return(i)
}
par.cor <- partial.cor(d)
par.cor
| Var 1 | Var 2 | Var 3 | Var 4 | Var 5 | |
|---|---|---|---|---|---|
| Var 1 | NA | 0.77164731 | 0.06605894 | -0.00737799 | 0.14615638 |
| Var 2 | 0.77164731 | NA | 0.09288327 | 0.04760085 | -0.02200392 |
| Var 3 | 0.06605894 | 0.09288327 | NA | 0.18522770 | 0.15392879 |
| Var 4 | -0.00737799 | 0.04760085 | 0.18522770 | NA | 0.32441680 |
| Var 5 | 0.14615638 | -0.02200392 | 0.15392879 | 0.32441680 | NA |
$x_1$ と $x_2$ をそのまま使用すると，結果は次のようになる。$x_1$ と $x_2$ の VIF（分散拡大要因）は大きなものではないが，$y$ との相関係数が 0.3，0.25 とほぼ同じなのに標準化偏回帰係数の絶対値は 7 倍も違い，符号は逆になっている。この状態が多重共線性か否かについては即断はできないが[^1]，解釈に困るのは確かであろう。
[^1:] $x_2, x_3, x_4$ を制御したときの $x_1$ と $y$ の偏相関係数は round(par.cor[5, 1], 3) = 0.146 であるが，$x_1, x_3, x_4$ を制御したときの $x_2$ と $y$ の偏相関係数は round(par.cor[5, 2], 3) = -0.022 なので，抑制変数の可能性はある。
round(par.cor[5, 1], 3); round(par.cor[5, 2], 3)
0.146
-0.022
x1 と x2 をそのまま使用する場合
ans1 <- lm(y~., d3)
summary(ans1)
Call:
lm(formula = y ~ ., data = d3)
Residuals:
Min 1Q Median 3Q Max
-5.0235 -1.1609 -0.0036 1.2032 4.7367
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.857031 1.136812 20.986 < 2e-16
x1 0.043517 0.013239 3.287 0.001084
x2 -0.007203 0.014710 -0.490 0.624581
x3 0.049655 0.014326 3.466 0.000574
x4 0.079569 0.010428 7.631 1.21e-13
Residual standard error: 1.763 on 495 degrees of freedom
Multiple R-squared: 0.2292, Adjusted R-squared: 0.2229
F-statistic: 36.79 on 4 and 495 DF, p-value: < 2.2e-16
$x_1$ と $x_2$ の相関が高いために困った事態が生じるのであるから，解決策は「どちらかだけを使う」ことであるのは明らかである。表に示すように $x_1$ と $x_2$ の標準化偏回帰係数の値は，それらと $y$ の相関係数に見合ったものであること，また，$x_3$，$x_4$ の標準化偏回帰係数はほぼ同じになることがわかる。
x1 だけを使用する場合
ans2 <- lm(y ~ x1+x3+x4, d3)
summary(ans2)
Call:
lm(formula = y ~ x1 + x3 + x4, data = d3)
Residuals:
Min 1Q Median 3Q Max
-5.082 -1.174 -0.007 1.201 4.726
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.843750 1.135617 20.996 < 2e-16
x1 0.038480 0.008328 4.621 4.89e-06
x3 0.049020 0.014256 3.438 0.000634
x4 0.079350 0.010410 7.622 1.28e-13
Residual standard error: 1.762 on 496 degrees of freedom
Multiple R-squared: 0.2288, Adjusted R-squared: 0.2241
F-statistic: 49.05 on 3 and 496 DF, p-value: < 2.2e-16
x2 だけを使用する場合
ans3 <-lm(y ~ x2+x3+x4, d3)
summary(ans3)
Call:
lm(formula = y ~ x2 + x3 + x4, data = d3)
Residuals:
Min 1Q Median 3Q Max
-5.251 -1.207 -0.027 1.252 4.827
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.860260 1.147992 20.784 < 2e-16
x2 0.030365 0.009352 3.247 0.001245
x3 0.053922 0.014408 3.743 0.000203
x4 0.081036 0.010521 7.703 7.3e-14
Residual standard error: 1.78 on 496 degrees of freedom
Multiple R-squared: 0.2123, Adjusted R-squared: 0.2076
F-statistic: 44.57 on 3 and 496 DF, p-value: < 2.2e-16
もう一つの解決法は，$x_1$ と $x_2$ の相関が高いということはどちらの変数もよく似たものを測定しているのだから両方足してしまう，つまり，$x_1+x_2$ を新たな独立変数として使うということである。この方法は，2 つの変数が平均値も分散もほとんど同じなら単純に和をとることにも違和感がないかもしれないが通常はそのようなことは期待できない。そこで，2 つの変数をそれぞれ標準化すれば平均値と分散は同じになるのだから足し算をすることができる。そして，単純和よりも望ましい方法は重み付けの和をとることであろう。
主成分回帰
さて，解決法が見えてきた。2 つの変数だけで主成分分析を行い，主成分得点を独立変数とすることである。2 つの変数の主成分分析からは 2 つの主成分得点が得られる。もとの 2 つの変数が持つ情報は 2 つの主成分が持つ情報と全く同じである。そして，2 つの主成分得点の相関は 0 である。
主成分分析による回転行列（重み）
pca <- prcomp(d[,1:2], scale.=TRUE)
rotation <- pca$rotation
colnames(rotation) <- paste("第", 1:2, "主成分", sep="")
rownames(rotation) <- paste("$x_", 1:2, "$", sep="")
rotation
| 第1主成分 | 第2主成分 | |
|---|---|---|
| $x_1$ | 0.7071068 | -0.7071068 |
| $x_2$ | 0.7071068 | 0.7071068 |
主成分得点
score <- pca$x
colnames(score) <- paste("第", 1:2, "主成分", sep="")
head(score, 10)
| 第1主成分 | 第2主成分 |
|---|---|
| -1.7291222 | 0.26897984 |
| -2.1987509 | -0.44874140 |
| 1.3980556 | 0.09472678 |
| 1.1624233 | -0.71042984 |
| -0.5687512 | -0.39831639 |
| 0.2989236 | 0.53771256 |
| -1.5403762 | -0.73871630 |
| 2.7795654 | 0.20310977 |
| -0.7901629 | 0.59175950 |
| 1.6691993 | 0.59851146 |
各変数の平均値と不偏分散および標準偏差
base2 <- data.frame(平均値=apply(score, 2, mean), 不偏分散=apply(score, 2, var), 標準偏差=apply(score, 2, sd))
rownames(base2) <- paste("第", 1:2, "主成分", sep="")
base2
| 平均値 | 不偏分散 | 標準偏差 | |
|---|---|---|---|
| <dbl> | <dbl> | <dbl> | |
| 第1主成分 | 2.819522e-15 | 1.8 | 1.3416408 |
| 第2主成分 | 2.793446e-15 | 0.2 | 0.4472136 |
主成分得点ともとの変数の相関
d2 <- d
d2[,1:2] <- score
colnames(d2)[1:2] <- c("$PC_1$", "$PC_2$")
cor(d2)
| $PC_1$ | $PC_2$ | $x_3$ | $x_4$ | $y$ | |
|---|---|---|---|---|---|
| $PC_1$ | 1.000000e+00 | -6.568541e-16 | 3.162278e-01 | 2.108185e-01 | 0.28987545 |
| $PC_2$ | -6.568541e-16 | 1.000000e+00 | -8.953297e-16 | -7.819641e-16 | -0.07905694 |
| $x_3$ | 3.162278e-01 | -8.953297e-16 | 1.000000e+00 | 3.000000e-01 | 0.30000000 |
| $x_4$ | 2.108185e-01 | -7.819641e-16 | 3.000000e-01 | 1.000000e+00 | 0.40000000 |
| $y$ | 2.898755e-01 | -7.905694e-02 | 3.000000e-01 | 4.000000e-01 | 1.00000000 |
では，$x_1$ と $x_2$ の代わりに 2 つの主成分得点を用いて重回帰分析を行ってみよう。結果を表\ref{tab7-10} に示す。
\SKIP{$PC_1$ と $PC_2$ の標準化偏回帰係数（偏回帰係数）の符号がマイナスになっているのは何の問題もない。第 1，第 2 主成分得点の符号を付け替えればよいだけである。見るべき所は，$PC_1$ と $PC_2$ に対する符号が一致しているということである。また，}
$x_3$，$x_4$ に対する標準化偏回帰係数は，前の表の標準化偏回帰係数とほぼ同じになっている。
x1 と x2 の主成分得点を使用する場合
d4 <- d2
colnames(d4) <- c("PC1", "PC2", "x3", "x4", "y")
ans4 <- lm(y ~ ., d4)
summary(ans4)
Call:
lm(formula = y ~ ., data = d4)
Residuals:
Min 1Q Median 3Q Max
-5.0235 -1.1609 -0.0036 1.2032 4.7367
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.70155 1.20597 21.312 < 2e-16
PC1 0.26187 0.06252 4.188 3.33e-05
PC2 -0.35355 0.17648 -2.003 0.045681
x3 0.04966 0.01433 3.466 0.000574
x4 0.07957 0.01043 7.631 1.21e-13
Residual standard error: 1.763 on 495 degrees of freedom
Multiple R-squared: 0.2292, Adjusted R-squared: 0.2229
F-statistic: 36.79 on 4 and 495 DF, p-value: < 2.2e-16
なお，このように独立変数の値をそのまま使うのではなく主成分分析によって得られる主成分得点を使って重回帰分析をするのは主成分回帰（Partial Compnent Regression; PCR）として確立されているものである。しかし，全ての独立変数を主成分分析する必要はなく，必要なもののみを主成分分析し主成分回帰を行うというのが自然なのである。