要因間の相関がもたらす問題

　重回帰では多数の説明変数を解析に組み込むことが可能である。こうすることで、同時に多数の説明変数と被説明変数の関係を明らかにできるわけだが、注意をしないとこの解析は大失敗をする。その原因の一つが本項で紹介する極めて強い多重共線性multicollinearityが変数間に存在する時である。多重共線性（マルチコとも呼ばれる）は、説明変数間に相関がある状況を指す。多重共線性があると、どのように解析を失敗してしまうのか、以下にシミュレーションしてみよう。

相関が低いときは問題が起こらない

　以下のような設定でデータを生成し、多重線形回帰を行ってみる。mvrnorm関数を使うことで多変量正規分布に従うデータを生成できる。多変量正規分布とは、各標本が正規分布に従いつつ、標本間に相関があるようなデータである（図1参照）。このとき、標本の分散と標本間の共分散を分散共分散行列で指定する。下記のスクリプトでは、vに格納したものが分散共分散行列である。vの対角要素が各標本の分散、非対角要素が標本間の共分散を指定している。

　注目してほしいのは、非対角要素=標本間の共分散を0としていることで、標本間には全くの相関がない状態にしてデータを生成している。この状態のときは、多重共線性の問題は生じない。

　切片b1は0、x1とx2の係数（b2、b3）は10、交互作用項b4は0として、被説明変数のデータも生成する。その後、多重線形回帰を2000回行う。係数の大きさとしては、十分な検出力が期待される。

------------------------------------------------------

library(plotn)

library(MASS) #mvrnorm関数を使うのに必要

b1 <- 0

b2 <- 10

b3 <- 10

b4 <- 0

s <- 5

n <- 60

ntest <- 2000

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

b1e <- NULL

b2e <- NULL

b3e <- NULL

b4e <- NULL

mu <- c(b2, b3)

v <- matrix(c(s, 0, 0, s), ncol = 2)#分散共分散行列

v

##      [,1] [,2]

## [1,]    5    0

## [2,]    0    5

for (i in 1:ntest){

  data <- mvrnorm(n, mu = mu, Sigma = v) #多変量正規分布に従うデータを生成

 #この関数はnにサンプルサイズ、muに生成する正規分布の平均を格納したベクトル

 #Sigmaに分散共分散行列を引数にとる

  x1 <- data[,1]

  x2 <- data[,2]

  y <- b4 * x1 * x2 + b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = 5)

  d <- data.frame(x1 = x1, x2 = x2, y = y)

  fit1 <- lm(y ~ x1 * x2, data = d)

  fit2 <- summary(fit1)

  b1e <- c(b1e, coef(fit2)[1,1])

  b2e <- c(b2e, coef(fit2)[2,1])

  b3e <- c(b3e, coef(fit2)[3,1])

  b4e <- c(b4e, coef(fit2)[4,1])

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

  p3 <- c(p3, coef(fit2)[3,4])

  p4 <- c(p4, coef(fit2)[4,4])

}

plotn(x2 ~ x1, main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図1の描画

------------------------------------------------------

データは大きくばらつき、一定の関係があるようには見えない。このとき、P値の分布は以下の通り。

------------------------------------------------------

histn(p1, xlab = "b1 P value", ylab = "頻度", 

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図2の描画

histn(p2, xlab = "b2 P value", ylab = "頻度",

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図3の描画

histn(p3, xlab = "b3 P value", ylab = "頻度",

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図4の描画

histn(p4, xlab = "b4 P value", ylab = "頻度",

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図5の描画

------------------------------------------------------

b1とb4は0を指定しているので、図2と5のようなP値の分布となる。一方、b2とb3は十分大きな値なのでP値が低くなる傾向が強く出ている。実際、検出力は100%に近い。

　また、各係数の推定値の分布は以下の通り。

------------------------------------------------------

histn(b1e, xlab = "estimate b1", ylab = "頻度",

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図6の描画

histn(b2e, xlab = "estimate b2", ylab = "頻度",

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図7の描画

histn(b3e, xlab = "estimate b3", ylab = "頻度",

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図8の描画

histn(b4e, xlab = "estimate b4", ylab = "頻度",

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図9の描画

plotn(b3e ~ b2e, xlab = "estimate b2", ylab = "estimate b3", 

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1))

------------------------------------------------------

特にb2とb3の推定値に注目すると、ともに設定どおり10付近を平均にとり、5～15くらいの範囲にほぼすべての推定値が存在していることが見て取れるだろう。実際、推定値の標準偏差は1.5くらいであり、正規分布の性質から10 - 1.5×3 = 5.5 ~ 14.5 = 10 + 1.5×3の範囲に99%のデータが存在することになる。

　合わせて、b2とb3の推定値の相関を確認すると、

------------------------------------------------------

plotn(b3e ~ b2e, xlab = "estimate b2", ylab = "estimate b3", 

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図10の描画

------------------------------------------------------

おおむね正の相関となっていることがわかる。以上のような状態が、正常に推定できているときである。

相関が極めて高いときに問題が起こる

　では、共分散を高めることで相関を高めて、問題を誘発してみよう。今度は、分散共分散行列の標本間の共分散を4.9995にした。このときに期待される標本間の相関係数は(標本1と2の共分散)/{(√標本1の分散)(√標本2の分散)} = 4.9995/(√5)(√5) = 4.9995/5 = 0.9999であり、極めて高い相関を示すデータを生成できる。

------------------------------------------------------

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

b1e <- NULL

b2e <- NULL

b3e <- NULL

b4e <- NULL

mu <- c(b2, b3)

v <- matrix(c(s, 0, 0, s), ncol = 2)#分散共分散行列

v

##        [,1]   [,2]

## [1,] 5.0000 4.9995

## [2,] 4.9995 5.0000

for (i in 1:ntest){

  data <- mvrnorm(n, mu = mu, Sigma = v) #多変量正規分布に従うデータを生成

 #この関数はnにサンプルサイズ、muに生成する正規分布の平均を格納したベクトル

 #Sigmaに分散共分散行列を引数にとる

  x1 <- data[,1]

  x2 <- data[,2]

  y <- b4 * x1 * x2 + b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = 5)

  d <- data.frame(x1 = x1, x2 = x2, y = y)

  fit1 <- lm(y ~ x1 * x2, data = d)

  fit2 <- summary(fit1)

  b1e <- c(b1e, coef(fit2)[1,1])

  b2e <- c(b2e, coef(fit2)[2,1])

  b3e <- c(b3e, coef(fit2)[3,1])

  b4e <- c(b4e, coef(fit2)[4,1])

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

  p3 <- c(p3, coef(fit2)[3,4])

  p4 <- c(p4, coef(fit2)[4,4])

}

plotn(x2 ~ x1, main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1)) #図11の描画

------------------------------------------------------

このように相関が高いとデータは一直線状にならび、ほとんどその直線からばらつかないことがわかる。今回に関しては標本の平均もそろえているので、x1のデータがわかればほぼそれに等しい値がx2で生成されているということだ。もし共分散を5、つまり相関係数を1にすると、データは全くばらつかず、x1とx2では全く同じデータが生成されることになる。さてこのとき、P値の分布は以下のようになる。

------------------------------------------------------

histn(p1, xlab = "b1 P value", ylab = "頻度", 

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図12の描画

histn(p2, xlab = "b2 P value", ylab = "頻度",

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図13の描画

histn(p3, xlab = "b3 P value", ylab = "頻度",

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図14の描画

histn(p4, xlab = "b4 P value", ylab = "頻度",

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図15の描画

------------------------------------------------------

すると、共分散以外の設定は変えていないにもかかわらず、b2とb3の検出力は極めて悪化し、8%になる。このように強い多重共線性は検出力を低下させてしまう。その原因は何だろうか。係数の推定値の分布を見てみよう。

------------------------------------------------------

histn(b1e, xlab = "estimate b1", ylab = "頻度",

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図16の描画

histn(b2e, xlab = "estimate b2", ylab = "頻度",

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図17の描画

histn(b3e, xlab = "estimate b3", ylab = "頻度",

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図18の描画

histn(b4e, xlab = "estimate b4", ylab = "頻度",

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図19の描画

plotn(b3e ~ b2e, xlab = "estimate b2", ylab = "estimate b3", 

      main = "cov(x1,x2) = 0", mar = c(3.8, 3.8, 2.5, 1))

------------------------------------------------------

すると、共分散が小さかった時と比べて、b2とb3は平均は10付近であるものの、そのばらつきは極めて大きくなり、-50～70となってしまった。実際の標準偏差は21.2程度であり、このことから推定値は10 - 21.2×3 = -53.6 ~ 73.6 = 10 + 21.2×3にデータのほとんどが存在する、という状況だ。合わせて、b2とb3の推定値の相関を確認すると、

------------------------------------------------------

plotn(b3e ~ b2e, xlab = "estimate b2", ylab = "estimate b3", 

      main = "cov(x1,x2) = 4.9995", mar = c(3.8, 3.8, 2.5, 1)) #図20の描画

------------------------------------------------------

すると、係数間に極めて強い負の相関が確認できる。実は推定された係数のペアの和は大体20程度になるのだが、一方の推定値が大きくぶれると、それを調整するようにもう一方の推定値も大きくぶれるのである。つまり、係数の推定が安定しないということである。

　説明変数間に強い相関があることで、係数の推定が安定しなくなる原因を2通りで説明してみよう。一つ目は直感的な説明である。今、相関の強い2変数をx1とx2として、これらを説明変数に組み込むことを考える。2変数は相関するので、誤差をuとして

　x2 = β x1 + α + u

とあらわせる。極めて強い相関を示すので、誤差は極めて小さく、実質、

　x2 = β x1 + α……(1)

である。ここで、交互作用を考えない多重線形回帰モデルは、

　y ~ b1 + b2 x1 + b3 x2……(2)

である。ここに(1)を代入すると

　y ~ b1 + b2 x1 + b3 (β x1 + α)

整理して、

　y ~ (b1 + α b3) + (b2 + β b3) x1……(3)

となる。b1 + α b3 = b1'、b2 + β b3 = b2'と置きなおし、回帰してある係数が決まったとしよう。

　^b1' = ^b1 + α (^b3)……(4)

　^b2' = ^b2 + β (^b3)……(5)

ところでこの問題は、推定値^b1～^b3を得る問題だったはずだ。しかし、3つの未知数に対し、式は2つしかなく、一意に未知数を決定できない。もっと言えば、(4)(5)から^b3を消去して

　^b2' = ^b2 + (β/α) (^b1' － ^b1)

整理して、

　^b2 － (β/α) (^b1) = ^b2' －(β/α) (^b1')  = (定数)……(6)

とでき、(6)は^b3を含まないので、どんな^b3をとることも可能である。^b3を決定できたとき、はじめて残りの推定値を決定できる。実際はわずかに誤差があるので何かしらの答えは出るが、わずかな誤差の違いで大きく答えが変わってしまうのだ。

　もうひとつ、もう少し数学的に解説しよう。多重線形回帰の項でパラメータは次のように計算できるのだった。

　^b = {(X^T)X}^(-1) (X^T) y

パラメータの推定のために(X^T)Xの逆行列の計算が必要になる。ここで(X^T)Xの実態は何だったかというと、以下のような行列だった。

とりあえず、3パラメータを推定する状況を考える。直感的説明と同様に、x2はx1で以下のようにあらわせるとする。

　　x2,i = β x1,i + α……(1)

すると、上記の行列は、

行列式determinantを|A|と表せば、ある行列Aについて逆行列が存在するためには|A| ≠ 0である必要がある。行列式の計算を詳しく紹介しないが、行列式の値はもとの行列の行間や列間で足し引きしても変わらない性質がある。また、その計算の結果、ある行や列を0にすることができれば、その行列式は0になる。したがって、上記の行列の行列式は下記のように計算できる。

このように変数間に極めて強い相関があると、行列式が0になり、逆行列が定義できない。つまり、パラメータの推定に失敗するわけである。実際は、わずかな誤差があるため0にはならないが、行列式は0に極めて近くなる。 逆行列はもとの行列に関して行列式で除した値で構成されているので、行列式が0に近い値をとるということは、行列式の各要素は大きな値を示すようになる。{(X^T)X}^(-1) にσ^2を掛けたものがパラメータの分散だったから、多重共線性があるとパラメータの分散が大きくなるのは必然なのだ。

様々な相関でシミュレートしてみる

　では、以下のように様々な相関で推定値のばらつきや検出力がどう変化するか確認しよう。

------------------------------------------------------

b1 <- 0

b2 <- 10

b3 <- 10

b4 <- 0

s <- 5

n <- 60

ntest <- 2000

r <- c(seq(0, 0.9, by = 0.1), 0.99, 0.999, 0.9999)#シミュレートする相関係数

m2 <- NULL

s2 <- NULL

m3 <- NULL

s3 <- NULL

p1t <- NULL

p2t <- NULL

p3t <- NULL

p4t <- NULL

for(j in r){

  p1 <- NULL

  p2 <- NULL

  p3 <- NULL

  p4 <- NULL

  b1e <- NULL

  b2e <- NULL

  b3e <- NULL

  b4e <- NULL

  mu <- c(b2, b3)

  v <- matrix(c(s, j*s, j*s, s), ncol = 2)

  #相関係数に分散を掛けたものが共分散になる。

  for (i in 1:ntest){

    data <- mvrnorm(n, mu = mu, Sigma = v)

    x1 <- data[,1]

    x2 <- data[,2]

    y <- b4 * x1 * x2 + b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = 5)

    d <- data.frame(x1 = x1, x2 = x2, y = y)

    fit1 <- lm(y ~ x1 * x2, data = d)

    fit2 <- summary(fit1)

    b1e <- c(b1e, coef(fit2)[1,1])

    b2e <- c(b2e, coef(fit2)[2,1])

    b3e <- c(b3e, coef(fit2)[3,1])

    b4e <- c(b4e, coef(fit2)[4,1])

    p1 <- c(p1, coef(fit2)[1,4])

    p2 <- c(p2, coef(fit2)[2,4])

    p3 <- c(p3, coef(fit2)[3,4])

    p4 <- c(p4, coef(fit2)[4,4])

  }

  p1t <- c(p1t, sum(p1 < 0.05))#b1の有意なP値の個数の計算

  p2t <- c(p2t, sum(p2 < 0.05))#b2の有意なP値の個数の計算

  p3t <- c(p3t, sum(p3 < 0.05))#b3の有意なP値の個数の計算

  p4t <- c(p4t, sum(p4 < 0.05))#b4の有意なP値の個数の計算

  m2 <- c(m2, mean(b2e))#b2の推定値の平均

  m3 <- c(m3, mean(b3e))#b3の推定値の標準偏差

  s2 <- c(s2, sd(b2e))#b2の推定値の平均

  s3 <- c(s3, sd(b3e))#b3の推定値の標準偏差

}

plotn(0, 0, xlab = "cov(x1, x2)", ylab = "Average of estimated b2", col.dot = NA,

      xlim = range(r*s), ylim = c(mean(m2) - max(s2), mean(m2) + max(s2)))#図21

overdraw(Mean_pt(cbind(m2,s2), at = r*s, SD = T, pch = 16, cex = 1.2, length = 0.1))

plotn(0, 0, xlab = "cov(x1, x2)", ylab = "Average of estimated b3", col.dot = NA,

      xlim = range(r*s), ylim = c(mean(m3) - max(s3), mean(m3) + max(s3)))#図22

overdraw(Mean_pt(cbind(m3,s3), at = r*s, SD = T, pch = 16, cex = 1.2, length = 0.1))

s2

##  [1]  1.504063  1.471013  1.448153  1.464267  1.413339  1.385706  1.327566

##  [8]  1.312858  1.282832  1.304819  2.324880  6.900988 21.195559

s3

##  [1]  1.522712  1.471367  1.431465  1.466504  1.412916  1.403202  1.322107

##  [8]  1.301906  1.281413  1.294064  2.351096  6.761350 21.122618

------------------------------------------------------

確認すると共分散に対して推定値の平均は10と一定のままだが、標準偏差は5付近に近づくと劇的に大きくなる。

------------------------------------------------------

cv <- r*s

plotn(p2t/ntest ~ cv, xlab = "cov(x1, x2)", ylab = "Power in b2")#図23の作図

plotn(p3t/ntest ~ cv, xlab = "cov(x1, x2)", ylab = "Power in b3")#図24の作図

p2t/ntest

##  [1] 0.9995 0.9990 1.0000 1.0000 1.0000 1.0000 1.0000 0.9995 1.0000 1.0000

## [11] 0.9840 0.3200 0.0800

p3t/ntest

##  [1] 0.9990 0.9990 1.0000 1.0000 1.0000 1.0000 0.9995 1.0000 1.0000 1.0000

## [11] 0.9805 0.3025 0.0740

------------------------------------------------------

検出力は共分散が4.5くらいまでは100%近いが、5に極めて近くなると急激に低下する。このように説明変数間の強い相関は解析上危険になりうる。すべての変数を解析に入れるのではなく、解析の前に変数間の相関を確認しておく方が良いだろう。

実際の解析

　では、以下のようなデータを解析をしてみよう。まずはデータの図示から。

------------------------------------------------------

x1 <- c(19, 0, 15, 28, 5, 26, 21, 10, 13, 20, 10, 

        19, 19, 5, 5, -2, 11, 22, -9, 21, -6, 10, 

        26, 27, 28, 14, 4, 29, 18, 11, -6, -5, 9, 

        27, 13, -1, 15, 11, 8, 26, 27, 24, 2, -10, 

        27, 13, 4, 9, 11, -6, 9, -3, 3, 8, 5, 2, 20, -4, 7, 22)

x2 <- c(27, 8, 30, 47, 18, 35, 28, 23, 22, 33, 

        20, 23, 24, 10, 7, 17, 26, 36, 6, 32,

        0, 22, 42, 45, 47, 27, 14, 43, 26, 25,

        12, 6, 26, 39, 27, 9, 35, 21, 11, 33,

        28, 30, 17, 2, 42, 30, 10, 15, 24, 4,

        17, 7, 19, 18, 14, 17, 28, -6, 9, 31)

x3 <- c(7, 1, -18, -18, -7, 6, -18, -14, -6, 9, 

        -19, 3, -2, -7, -5, -18, -1, 3, -5, -8,

        -3, 9, -20, -14, -9, 6, 0, -14, -2, -7, 

        -4, -3, -2, -17, -13, 6, -17, 5, 2, -14,

        -17, 5, 2, -9, -12, 7, -3, -9, -5, 8, 5, 

        5, -14, -20, 6, -12, -17, -10, -8, -4)

x4 <- c(2, -4, -23, -23, -12, 1, -23, -19, -11, 4, 

        -23, -2, -7, -12, -10, -23, -6, -2, -10, 

        -13, -8, 4, -25, -19, -14, 2, -5, -19, -8,

        -12, -9, -8, -7, -22, -18, 1, -22, 0, -3, 

        -19, -22, 0, -3, -14, -17, 3, -8, -14, -10, 

        3, 0, 1, -19, -25, 1, -17, -22, -15, -13, -9)

y <- c(568, 64, 48, 373, 54, 706, 99, 20, 206, 663, 

       -90, 437, 343, -34, -2, -244, 326, 602, -157, 

       342, -141, 473, 244, 408, 545, 520, 143, 410, 

       355, 195, -39, -80, 280, 295, 103, 157, 128, 

       390, 194, 295, 169, 616, 191, -293, 415, 558,

       62, 33, 215, 108, 340, 129, -91, -165, 270, 

       -71, 114, -332, -38, 407)

d <- data.frame(x1 = x1, x2 = x2, x3 = x3, x4 = x4, y = y)

plotn(d)#図25の作図

------------------------------------------------------

するとx1とx2、x3とx4の間に相関がみられる。特にx3とx4は極めて相関が強そうだ。実際、x3とx4を同時に解析に入れると、x4側は有意になるがパラメータの標準誤差が大きくなる。

------------------------------------------------------

fit1 <- lm(y ~ x1 + x2 + x3 + x4, data = d)

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 + x2 + x3 + x4, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -14.4046  -5.7981   0.3572   5.1785  16.1578 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  62.2521    18.1404   3.432  0.00115 ** 

## x1            9.8142     0.2350  41.756  < 2e-16 ***

## x2           10.1448     0.2092  48.495  < 2e-16 ***

## x3            0.8741     3.6134   0.242  0.80976    

## x4           19.0064     3.5969   5.284 2.23e-06 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 7.767 on 55 degrees of freedom

## Multiple R-squared:  0.9991, Adjusted R-squared:  0.999 

## F-statistic: 1.498e+04 on 4 and 55 DF,  p-value: < 2.2e-16

------------------------------------------------------

こんな時はx3とx4のどちらかを解析から除外してもよいだろう。

------------------------------------------------------

fit3 <- lm(y ~ x1 + x2 + x4, data = d)

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = y ~ x1 + x2 + x4, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -14.4021  -5.7287   0.0615   5.2264  16.1096 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  66.5933     2.6233   25.39   <2e-16 ***

## x1            9.8233     0.2300   42.71   <2e-16 ***

## x2           10.1374     0.2052   49.41   <2e-16 ***

## x4           19.8760     0.1178  168.74   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 7.701 on 56 degrees of freedom

## Multiple R-squared:  0.9991, Adjusted R-squared:  0.999 

## F-statistic: 2.032e+04 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

するとx4の標準誤差も改善される。