多重共線性の回避策

## 多重共線性の回避策

独立変数 $x_4$ が従属変数 $y$ に与える影響を知りたい。ほかの独立変数 $x_1$ と $x_2$ の間に多重共線性があるが，これはどちらも説明変数に含めたい。$x_1$ と $x_4$，$x_2$ と $x_4$ の相関は低く VIF も小さい。

このような場合にどうしたらよいかという質問に対して，シミュレーションに基づく解答を提示する。
以下では R の場合におけるプログラム例であるが，他の言語，パッケージでも追試可能である。

変数間の相関係数行列は以下のような表のようになっていると仮定しよう。


```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
```


<table class="dataframe">
<caption>A matrix: 5 × 5 of type dbl</caption>
<thead>
	<tr><th></th><th scope=col>$x_1$</th><th scope=col>$x_2$</th><th scope=col>$x_3$</th><th scope=col>$x_4$</th><th scope=col>$y$</th></tr>
</thead>
<tbody>
	<tr><th scope=row>$x_1$</th><td>1.0</td><td>0.80</td><td>0.3</td><td>0.2</td><td>0.30</td></tr>
	<tr><th scope=row>$x_2$</th><td>0.8</td><td>1.00</td><td>0.3</td><td>0.2</td><td>0.25</td></tr>
	<tr><th scope=row>$x_3$</th><td>0.3</td><td>0.30</td><td>1.0</td><td>0.3</td><td>0.30</td></tr>
	<tr><th scope=row>$x_4$</th><td>0.2</td><td>0.20</td><td>0.3</td><td>1.0</td><td>0.40</td></tr>
	<tr><th scope=row>$y$</th><td>0.3</td><td>0.25</td><td>0.3</td><td>0.4</td><td>1.00</td></tr>
</tbody>
</table>



まず，表のような相関係数行列になるように，多変量正規分布にしたがうシミュレーションデータを作成する（$n=500$）。

先頭の 10 行は表\ref{tab7-2} のようになる（小数点以下 4 桁以降は表示していない）。


```R
# 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)
```


<table class="dataframe">
<caption>A data.frame: 10 × 5</caption>
<thead>
	<tr><th></th><th scope=col>$x_1$</th><th scope=col>$x_2$</th><th scope=col>$x_3$</th><th scope=col>$x_4$</th><th scope=col>$y$</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>1</th><td>35.87129</td><td>36.70771</td><td>73.36168</td><td>66.53958</td><td>33.12320</td></tr>
	<tr><th scope=row>2</th><td>37.62556</td><td>29.15146</td><td>82.96882</td><td>76.82022</td><td>33.95505</td></tr>
	<tr><th scope=row>3</th><td>59.21593</td><td>55.50001</td><td>77.28258</td><td>86.95731</td><td>36.32317</td></tr>
	<tr><th scope=row>4</th><td>63.24307</td><td>48.87647</td><td>84.25619</td><td>81.25347</td><td>36.78013</td></tr>
	<tr><th scope=row>5</th><td>48.79484</td><td>39.84562</td><td>64.62895</td><td>98.00760</td><td>39.80084</td></tr>
	<tr><th scope=row>6</th><td>48.31151</td><td>51.32432</td><td>89.10405</td><td>99.74816</td><td>39.74805</td></tr>
	<tr><th scope=row>7</th><td>44.33141</td><td>31.49594</td><td>74.81446</td><td>86.10918</td><td>37.23745</td></tr>
	<tr><th scope=row>8</th><td>68.21829</td><td>64.98163</td><td>76.57823</td><td>86.12623</td><td>36.14495</td></tr>
	<tr><th scope=row>9</th><td>40.22833</td><td>44.73737</td><td>81.74600</td><td>71.22765</td><td>33.55838</td></tr>
	<tr><th scope=row>10</th><td>57.57091</td><td>60.43162</td><td>80.29413</td><td>80.86618</td><td>35.98782</td></tr>
</tbody>
</table>



それぞれの基礎統計量は以下の通りである。


```R
base <- data.frame(平均値=apply(d, 2, mean), 不偏分散=apply(d, 2, var), 標準偏差=apply(d, 2, sd))
rownames(base) <- names
base
```


<table class="dataframe">
<caption>A data.frame: 5 × 3</caption>
<thead>
	<tr><th></th><th scope=col>平均値</th><th scope=col>不偏分散</th><th scope=col>標準偏差</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>$x_1$</th><td>50</td><td>100</td><td>10</td></tr>
	<tr><th scope=row>$x_2$</th><td>46</td><td> 81</td><td> 9</td></tr>
	<tr><th scope=row>$x_3$</th><td>76</td><td> 36</td><td> 6</td></tr>
	<tr><th scope=row>$x_4$</th><td>82</td><td> 64</td><td> 8</td></tr>
	<tr><th scope=row>$y$</th><td>36</td><td>  4</td><td> 2</td></tr>
</tbody>
</table>



散布図は図のようになり，相関係数行列は正確に表の通りになる。


```R
d3 <- d
colnames(d3) <- c(paste("x", 1:4, sep=""), "y")
plot(d3, cex=0.7, pch=16, col="#00aa0060")
```


    
![output_7_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/3e0508d7-032e-0f35-6be7-96642d95d323.png)

    



```R
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
```


<table class="dataframe">
<caption>A matrix: 5 × 5 of type dbl</caption>
<thead>
	<tr><th></th><th scope=col>Var 1</th><th scope=col>Var 2</th><th scope=col>Var 3</th><th scope=col>Var 4</th><th scope=col>Var 5</th></tr>
</thead>
<tbody>
	<tr><th scope=row>Var 1</th><td>         NA</td><td> 0.77164731</td><td>0.06605894</td><td>-0.00737799</td><td> 0.14615638</td></tr>
	<tr><th scope=row>Var 2</th><td> 0.77164731</td><td>         NA</td><td>0.09288327</td><td> 0.04760085</td><td>-0.02200392</td></tr>
	<tr><th scope=row>Var 3</th><td> 0.06605894</td><td> 0.09288327</td><td>        NA</td><td> 0.18522770</td><td> 0.15392879</td></tr>
	<tr><th scope=row>Var 4</th><td>-0.00737799</td><td> 0.04760085</td><td>0.18522770</td><td>         NA</td><td> 0.32441680</td></tr>
	<tr><th scope=row>Var 5</th><td> 0.14615638</td><td>-0.02200392</td><td>0.15392879</td><td> 0.32441680</td><td>         NA</td></tr>
</tbody>
</table>



$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 なので，抑制変数の可能性はある。


```R
round(par.cor[5, 1], 3); round(par.cor[5, 2], 3)
```


0.146



-0.022


### x1 と x2 をそのまま使用する場合


```R
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 だけを使用する場合


```R
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 だけを使用する場合


```R
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 である。

### 主成分分析による回転行列（重み）


```R
pca <-  prcomp(d[,1:2], scale.=TRUE)
rotation <- pca$rotation
colnames(rotation) <- paste("第", 1:2, "主成分", sep="")
rownames(rotation) <-  paste("$x_", 1:2, "$", sep="")
rotation
```


<table class="dataframe">
<caption>A matrix: 2 × 2 of type dbl</caption>
<thead>
	<tr><th></th><th scope=col>第1主成分</th><th scope=col>第2主成分</th></tr>
</thead>
<tbody>
	<tr><th scope=row>$x_1$</th><td>0.7071068</td><td>-0.7071068</td></tr>
	<tr><th scope=row>$x_2$</th><td>0.7071068</td><td> 0.7071068</td></tr>
</tbody>
</table>



### 主成分得点


```R
score <- pca$x
colnames(score) <- paste("第", 1:2, "主成分", sep="")
head(score, 10)
```


<table class="dataframe">
<caption>A matrix: 10 × 2 of type dbl</caption>
<thead>
	<tr><th scope=col>第1主成分</th><th scope=col>第2主成分</th></tr>
</thead>
<tbody>
	<tr><td>-1.7291222</td><td> 0.26897984</td></tr>
	<tr><td>-2.1987509</td><td>-0.44874140</td></tr>
	<tr><td> 1.3980556</td><td> 0.09472678</td></tr>
	<tr><td> 1.1624233</td><td>-0.71042984</td></tr>
	<tr><td>-0.5687512</td><td>-0.39831639</td></tr>
	<tr><td> 0.2989236</td><td> 0.53771256</td></tr>
	<tr><td>-1.5403762</td><td>-0.73871630</td></tr>
	<tr><td> 2.7795654</td><td> 0.20310977</td></tr>
	<tr><td>-0.7901629</td><td> 0.59175950</td></tr>
	<tr><td> 1.6691993</td><td> 0.59851146</td></tr>
</tbody>
</table>



### 各変数の平均値と不偏分散および標準偏差


```R
base2 <- data.frame(平均値=apply(score, 2, mean), 不偏分散=apply(score, 2, var), 標準偏差=apply(score, 2, sd))
rownames(base2) <- paste("第", 1:2, "主成分", sep="")
base2
```


<table class="dataframe">
<caption>A data.frame: 2 × 3</caption>
<thead>
	<tr><th></th><th scope=col>平均値</th><th scope=col>不偏分散</th><th scope=col>標準偏差</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>第1主成分</th><td>2.819522e-15</td><td>1.8</td><td>1.3416408</td></tr>
	<tr><th scope=row>第2主成分</th><td>2.793446e-15</td><td>0.2</td><td>0.4472136</td></tr>
</tbody>
</table>



### 主成分得点ともとの変数の相関


```R
d2 <- d
d2[,1:2] <- score
colnames(d2)[1:2] <- c("$PC_1$", "$PC_2$")
cor(d2)
```


<table class="dataframe">
<caption>A matrix: 5 × 5 of type dbl</caption>
<thead>
	<tr><th></th><th scope=col>$PC_1$</th><th scope=col>$PC_2$</th><th scope=col>$x_3$</th><th scope=col>$x_4$</th><th scope=col>$y$</th></tr>
</thead>
<tbody>
	<tr><th scope=row>$PC_1$</th><td> 1.000000e+00</td><td>-6.568541e-16</td><td> 3.162278e-01</td><td> 2.108185e-01</td><td> 0.28987545</td></tr>
	<tr><th scope=row>$PC_2$</th><td>-6.568541e-16</td><td> 1.000000e+00</td><td>-8.953297e-16</td><td>-7.819641e-16</td><td>-0.07905694</td></tr>
	<tr><th scope=row>$x_3$</th><td> 3.162278e-01</td><td>-8.953297e-16</td><td> 1.000000e+00</td><td> 3.000000e-01</td><td> 0.30000000</td></tr>
	<tr><th scope=row>$x_4$</th><td> 2.108185e-01</td><td>-7.819641e-16</td><td> 3.000000e-01</td><td> 1.000000e+00</td><td> 0.40000000</td></tr>
	<tr><th scope=row>$y$</th><td> 2.898755e-01</td><td>-7.905694e-02</td><td> 3.000000e-01</td><td> 4.000000e-01</td><td> 1.00000000</td></tr>
</tbody>
</table>



では，$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 の主成分得点を使用する場合


```R
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）として確立されているものである。しかし，全ての独立変数を主成分分析する必要はなく，必要なもののみを主成分分析し主成分回帰を行うというのが自然なのである。