独立変数 \(x_4\) が従属変数 \(y\) に与える影響を知りたい。ほかの独立変数 \(x_1\) と \(x_2\) の間に多重共線性があるが，これはどちらも説明変数に含めたい。\(x_1\) と \(x_4\)，\(x_2\) と \(x_4\) の相関は低く VIF も小さい。
このような場合にどうしたらよいかという質問に対して，シミュレーションに基づく解答を提示する。
変数間の相関係数行列は表 1 のようになっていると仮定しよう。
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(y\) | |
|---|---|---|---|---|---|
| \(x_1\) | 1.00 | 0.80 | 0.30 | 0.20 | 0.30 |
| \(x_2\) | 0.80 | 1.00 | 0.30 | 0.20 | 0.25 |
| \(x_3\) | 0.30 | 0.30 | 1.00 | 0.30 | 0.30 |
| \(x_4\) | 0.20 | 0.20 | 0.30 | 1.00 | 0.40 |
| \(y\) | 0.30 | 0.25 | 0.30 | 0.40 | 1.00 |
まず，表 1 のような相関係数行列になるように，多変量正規分布に従うシミュレーションデータを作成する（\( n=500\)）。
先頭の 10 行は表 2 のようになる（小数点以下4桁以降は表示していない）。
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(y\) | |
|---|---|---|---|---|---|
| 1 | 35.871 | 36.708 | 73.362 | 66.540 | 33.123 |
| 2 | 37.626 | 29.151 | 82.969 | 76.820 | 33.955 |
| 3 | 59.216 | 55.500 | 77.283 | 86.957 | 36.323 |
| 4 | 63.243 | 48.876 | 84.256 | 81.253 | 36.780 |
| 5 | 48.795 | 39.846 | 64.629 | 98.008 | 39.801 |
| 6 | 48.312 | 51.324 | 89.104 | 99.748 | 39.748 |
| 7 | 44.331 | 31.496 | 74.814 | 86.109 | 37.237 |
| 8 | 68.218 | 64.982 | 76.578 | 86.126 | 36.145 |
| 9 | 40.228 | 44.737 | 81.746 | 71.228 | 33.558 |
| 10 | 57.571 | 60.432 | 80.294 | 80.866 | 35.988 |
それぞれの基礎統計量は表 3 の通りである。
| 平均値 | 不偏分散 | 標準偏差 | |
|---|---|---|---|
| \(x_1\) | 50.000 | 100.000 | 10.000 |
| \(x_2\) | 46.000 | 81.000 | 9.000 |
| \(x_3\) | 76.000 | 36.000 | 6.000 |
| \(x_4\) | 82.000 | 64.000 | 8.000 |
| \(y\) | 36.000 | 4.000 | 2.000 |
散布図は図 1 のようになり，相関係数行列は正確に表 1 の通りになる。
\(x_1\) と \(x_2\) をそのまま使用すると，結果は表 4 のようになる。\(x_1\) と \(x_2\) の VIF（分散拡大要因）は大きなものではないが，\(y\) との相関係数が 0.3，0.25 とほぼ同じなのに標準化偏回帰係数の絶対値は 7 倍も違い，符号は逆になっている。この状態が多重共線性か否かについては即断はできないが，解釈に困るのは確かであろう。
| 偏回帰係数 | 標準誤差 | \(t\) 値 | \(P\) 値 | 標準化偏回帰係数 | VIF | |
|---|---|---|---|---|---|---|
| \(x_1\) | 0.044 | 0.013 | 3.287 | 0.001 | 0.218 | 2.814 |
| \(x_2\) | -0.007 | 0.015 | -0.490 | 0.625 | -0.032 | 2.814 |
| \(x_3\) | 0.050 | 0.014 | 3.466 | 0.001 | 0.149 | 1.186 |
| \(x_4\) | 0.080 | 0.010 | 7.631 | 0.000 | 0.318 | 1.117 |
| 定数項 | 23.857 | 1.137 | 20.986 | 0.000 |
\(x_1\) と \(x_2\) の相関が高いために困った事態が生じるのであるから，解決策は「どちらかだけを使う」ことであるのは明らかである。表 5，6 に示すように \(x_1\) と \(x_2\) の標準化偏回帰係数の値は，それらと \(y\) の相関係数に見合ったものであること，また，\(x_3\)，\(x_4\) の標準化偏回帰係数はほぼ同じになることがわかる。
| 偏回帰係数 | 標準誤差 | \(t\) 値 | \(P\) 値 | 標準化偏回帰係数 | VIF | |
|---|---|---|---|---|---|---|
| \(x_1\) | 0.038 | 0.008 | 4.621 | 0.000 | 0.192 | 1.115 |
| \(x_3\) | 0.049 | 0.014 | 3.438 | 0.001 | 0.147 | 1.176 |
| \(x_4\) | 0.079 | 0.010 | 7.622 | 0.000 | 0.317 | 1.115 |
| 定数項 | 23.844 | 1.136 | 20.996 | 0.000 |
| 偏回帰係数 | 標準誤差 | \(t\) 値 | \(P\) 値 | 標準化偏回帰係数 | VIF | |
|---|---|---|---|---|---|---|
| \(x_2\) | 0.030 | 0.009 | 3.247 | 0.001 | 0.137 | 1.115 |
| \(x_3\) | 0.054 | 0.014 | 3.743 | 0.000 | 0.162 | 1.176 |
| \(x_4\) | 0.081 | 0.011 | 7.703 | 0.000 | 0.324 | 1.115 |
| 定数項 | 23.860 | 1.148 | 20.784 | 0.000 |
もう一つの解決法は，\(x_1\) と \(x_2\) の相関が高いということはどちらの変数もよく似たものを測定しているのだから両方足してしまう，つまり，\(x_1+x_2\) を新たな独立変数として使うということである。この方法は，2 つの変数が平均値も分散もほとんど同じなら単純に和をとることにも違和感がないかもしれないが通常はそのようなことは期待できない。そこで，2 つの変数をそれぞれ標準化すれば平均値と分散は同じになるのだから足し算をすることができる。そして，単純和よりも望ましい方法は重み付けの和をとることであろう。
さて，解決法が見えてきた。2 つの変数だけで主成分分析を行い，主成分得点を独立変数とすることである。2 つの変数の主成分分析からは 2 つの主成分得点が得られる。もとの 2 つの変数が持つ情報は 2 つの主成分が持つ情報と全く同じである。そして，2 つの主成分得点の相関は 0 である。
| 第1主成分 | 第2主成分 | |
|---|---|---|
| \(x_1\) | 0.70711 | -0.70711 |
| \(x_2\) | 0.70711 | 0.70711 |
| 第1主成分 | 第2主成分 | |
|---|---|---|
| 1 | -1.72912 | 0.26898 |
| 2 | -2.19875 | -0.44874 |
| 3 | 1.39806 | 0.09473 |
| 4 | 1.16242 | -0.71043 |
| 5 | -0.56875 | -0.39832 |
| 6 | 0.29892 | 0.53771 |
| 7 | -1.54038 | -0.73872 |
| 8 | 2.77957 | 0.20311 |
| 9 | -0.79016 | 0.59176 |
| 10 | 1.66920 | 0.59851 |
| 平均値 | 不偏分散 | 標準偏差 | |
|---|---|---|---|
| 第1主成分 | 0.00000 | 1.80000 | 1.34164 |
| 第2主成分 | -0.00000 | 0.20000 | 0.44721 |
| \(PC_1\) | \(PC_2\) | \(x_3\) | \(x_4\) | \(y\) | |
|---|---|---|---|---|---|
| \(PC_1\) | 1.00000 | -0.00000 | 0.31623 | 0.21082 | 0.28988 |
| \(PC_2\) | -0.00000 | 1.00000 | -0.00000 | -0.00000 | -0.07906 |
| \(x_3\) | 0.31623 | -0.00000 | 1.00000 | 0.30000 | 0.30000 |
| \(x_4\) | 0.21082 | -0.00000 | 0.30000 | 1.00000 | 0.40000 |
| \(y\) | 0.28988 | -0.07906 | 0.30000 | 0.40000 | 1.00000 |
では，\(x_1\) と \(x_2\) の代わりに 2 つの主成分得点を用いて重回帰分析を行ってみよう。結果を表 11 に示す。\(PC_1\) と \(PC_2\) の標準化偏回帰係数（偏回帰係数）の符号がマイナスになっているのは何の問題もない。第 1，第 2 主成分得点の符号を付け替えればよいだけである。見るべき所は，\(PC_1\) と \(PC_2\) に対する符号が一致しているということである。また，\(x_3\)，\(x_4\) に対する標準化偏回帰係数は，表 4，5，6 の標準化偏回帰係数とほぼ同じになっている。
| 偏回帰係数 | 標準誤差 | \(t\) 値 | \(P\) 値 | 標準化偏回帰係数 | VIF | |
|---|---|---|---|---|---|---|
| \(PC_1\) | 0.262 | 0.063 | 4.188 | 0.000 | 0.176 | 1.130 |
| \(PC_2\) | -0.354 | 0.176 | -2.003 | 0.046 | -0.079 | 1.000 |
| \(x_3\) | 0.050 | 0.014 | 3.466 | 0.001 | 0.149 | 1.186 |
| \(x_4\) | 0.080 | 0.010 | 7.631 | 0.000 | 0.318 | 1.117 |
| 定数項 | 25.702 | 1.206 | 21.312 | 0.000 |
なお，このように独立変数の値をそのまま使うのではなく主成分分析によって得られる主成分得点を使って重回帰分析をするのは主成分回帰（Partial Compnent Regression; PCR）として確立されているものである。しかし，全ての独立変数を主成分分析する必要はなく，逆に，必要なもののみを主成分分析し主成分回帰を行うというのが自然なのである。
options(digits=7, width=100) 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 print(xtable(r, caption="表 1. 相関係数行列"), caption.placement="top", type="html") library(MASS) library(xtable) source("http://aoki2.si.gunma-u.ac.jp/R/src/mreg.R", encoding="euc-jp") 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 print(xtable(head(d, 10), caption="表 2. データ行列", digits=3), caption.placement="top", type="html") base <- data.frame(平均値=apply(d, 2, mean), 不偏分散=apply(d, 2, var), 標準偏差=apply(d, 2, sd)) rownames(base) <- names print(xtable(base, caption="表 3. 各変数の平均値と不偏分散および標準偏差", digits=3), caption.placement="top", type="html") d3 <- d colnames(d3) <- c(paste("x", 1:4, sep=""), "y") plot(d3, cex=0.7, pch=19, col="#00008820") ans1 <- data.frame(mreg(d)$result) colnames(ans1)[c(3:4,6)] <- c("\\(t\\) 値", "\\(P\\) 値", "VIF") rownames(ans1)[1:4] <- paste("\\(x_", 1:4, "\\)", sep="") ans1[6] <- 1/ans1[6] print(xtable(ans1, caption="表 4. \\(x_1\\) と \\(x_2\\) をそのまま使用する場合", digits=3), caption.placement="top", type="html") ans2 <- data.frame(mreg(d[,c(1,3,4,5)])$result) colnames(ans2)[c(3:4,6)] <- c("\\(t\\) 値", "\\(P\\) 値", "VIF") rownames(ans2)[1:3] <- paste("\\(x_", c(1, 3, 4), "\\)", sep="") ans2[6] <- 1/ans2[6] print(xtable(ans2, caption="表 5. \\(x_1\\) だけを使用する場合", digits=3), caption.placement="top", type="html") ans3 <- data.frame(mreg(d[,2:5])$result) colnames(ans3)[c(3:4,6)] <- c("\\(t\\) 値", "\\(P\\) 値", "VIF") rownames(ans3)[1:3] <- paste("\\(x_", c(2, 3, 4), "\\)", sep="") ans3[6] <- 1/ans3[6] print(xtable(ans3, caption="表 6. \\(x_2\\) だけを使用する場合", digits=3), caption.placement="top", type="html") pca <- prcomp(d[,1:2], scale.=TRUE) rotation <- pca$rotation colnames(rotation) <- paste("第", 1:2, "主成分", sep="") rownames(rotation) <- paste("\\(x_", 1:2, "\\)", sep="") print(xtable(rotation, caption="表 7. 主成分分析による回転行列（重み）", digits=5), caption.placement="top", type="html") score <- pca$x colnames(score) <- paste("第", 1:2, "主成分", sep="") print(xtable(head(score, 10), caption="表 8. 主成分得点", digits=5), caption.placement="top", type="html") base2 <- data.frame(平均値=apply(score, 2, mean), 不偏分散=apply(score, 2, var), 標準偏差=apply(score, 2, sd)) rownames(base2) <- paste("第", 1:2, "主成分", sep="") print(xtable(base2, caption="表 9. 各変数の平均値と不偏分散および標準偏差", digits=5), caption.placement="top", type="html") d2 <- d d2[,1:2] <- score colnames(d2)[1:2] <- c("\\(PC_1\\)", "\\(PC_2\\)") print(xtable(cor(d2), caption="表 10. 主成分得点ともとの変数の相関", digits=5), caption.placement="top", type="html") ans4 <- data.frame(mreg(d2)$result) colnames(ans4)[c(3:4,6)] <- c("\\(t\\) 値", "\\(P\\) 値", "VIF") rownames(ans4)[1:2] <- paste("\\(PC_", 1:2, "\\)", sep="") rownames(ans4)[3:4] <- paste("\\(x_", 3:4, "\\)", sep="") ans4[6] <- 1/ans4[6] print(xtable(ans4, caption="表 11. \\(x_1\\) と \\(x_2\\) の主成分得点を使用する場合", digits=3), caption.placement="top", type="html")