井出草平の研究ノート

多重共線性のシミュレーション

下記エントリーの続き。

ides.hatenablog.com

こちらの教科書から多重共線性について

  • Richard McElreath - Statistical Rethinking_ A Bayesian Course with Examples in R and STAN

6.1. 多重共線性

一般に、回帰モデルに追加する潜在的な予測変数が多くあることは事実である。たとえば、霊長類のミルク・データの場合、我々が結果として選ぶどの列も予測するために利用可能な7つの変数がある。なぜ、7つすべてを含むモデルを構築しないのか? いくつかの危険性がある。

まず、最も心配のない「多重共線性」から始めよう。多重共線性とは、2つ以上の予測変数の間に非常に強い関連があることを意味する。生の相関は重要ではない。むしろ重要なのは、モデル中の他の変数に対する条件付きで、関連性があることである。多重共線性の結果は,たとえすべての変数が実際には結果と強く関連していても,事後分布は,どの変数も結果と確実に関連していないことを示唆しているように見えるということである。 このイライラする現象は、重回帰がどのように機能するかの詳細から発生する。実際、多重共線性には何の問題もない。モデルは予測には問題なく機能するす。ただ、それを理解しようとするとフラストレーションがたまるだけである。多重共線性を理解すれば、一般的な回帰モデルをよりよく理解できるようになることが期待できる。

まずは簡単なシミュレーションから始めよう。その後、再び霊長類のミルクのデータに目を向け、実際のデータセットで多重共線性を見てみよう。

6.1.1. 多重共線性 脚

個人の足の長さを予測変数として、その人の身長を予測しようとすることを想像しよう。確かに身長は脚の長さと正の相関がある、少なくとも我々のシミュレーションではそうだと仮定する。しかし、両足(右足と左足)をモデルに入れると、やっかいなことが起こる。

以下のコードは、100人の身長と脚の長さをシミュレートしている。それぞれについて、まず身長がガウス分布からシミュレートされる。次に、各個体の脚の高さの比率を0.4から0.5までの範囲でシミュレートする。最後に、それぞれの脚は、実際の集団で典型的に見られるように、左右の脚の長さが全く同じにならないように、若干の測定誤差や発育誤差を加味してある。

最後に、このコードは身長と2本の足の長さを共通のデータフレームに入れる。

N <- 100   # 個体数 
set.seed(909)
height <- rnorm(N,10,2)  # 各脚の高さの合計をシミュレートする
leg_prop <- runif(N,0.4,0.5) # 脚長比 
# 左脚を比例+誤差としてシミュレート
leg_left <- leg_prop*height +    rnorm( N , 0 , 0.02 )
# 右足を比率+誤差としてシミュレート
leg_right <- leg_prop*height + rnorm( N , 0 , 0.02 ) 
# 結合してデータフレームにする
d <- data.frame(height,leg_left,leg_right)

さて、これらのデータを分析して、leg_leftとleg_rightの両方の予測変数で結果の高さを予測しようしかし、事後推定を行う前に、我々が何を期待するかを考えてみよう平均して、個人の脚は身長の45%である(このシミュレーション・データでは)。したがって、脚と身長の関連性を測るベータ係数は、平均身長(10)を平均身長の45%(4.5)で割ったあたりになると予想される。これは、10/4.5≒2.2です。では、代わりにどうなるか見てみよう。ここでは、これから起こることの責任が事前分布にないないことを確認するために、非常に曖昧で悪い事前分布を使う。

library(rethinking)
m6.1 <- quap(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + bl*leg_left + br*leg_right , a ~ dnorm( 10 , 100 ) ,
    bl ~ dnorm( 2 , 10 ) , br ~ dnorm( 2 , 10 ) , sigma ~ dexp( 1 )
  ) , data=d )
precis(m6.1)
      mean   sd  5.5% 94.5%
a     0.98 0.28  0.53  1.44
bl    0.21 2.53 -3.83  4.25
br    1.78 2.53 -2.26  5.83
sigma 0.62 0.04  0.55  0.69

この事後平均と標準偏差は狂っているように見える。これはprecis出力のグラフ表示の方が有用なケースである。なぜなら、事後平均と89%区間を一目見て、ここで何かが間違っていることがわかるように表示されるからである。

plot(precis(m6.1))

set.seedの行を省いて、もう何回かシミュレーションをしてほしい。両足の長さがほぼ同じで、身長が足の長さと強く関連しているのなら、この事後分布はなぜこんなに変なのだろうか?事後近似は正しく機能したのだろううか?

正しく動作している。そして、この事後分布は私たちが質問した内容に対する正しい答えである。問題は、その質問である。重回帰は質問に答えるものであることを思い出してほしい他の予測変数がすべてわかっている状態で、各予測変数を知ることの価値は何だろうか? 従って、この場合、質問は次のようになる。各脚の長さを知っていることの価値は何か?

post <- extract.samples(m6.1)
plot( bl ~ br , post , col=col.alpha(rangi2,0.1) , pch=16 )

図6.2. 左:各脚と身長の関連性の事後分布(モデルm6.1より)。両変数はほぼ同じ情報を含むので、事後は負の相関のある値の狭い隆起である。右図。2つのパラメータの和の事後分布は、どちらかの脚と身長の適切な関連性で中央化される。

その結果、図6.2の左側に示すようなプロットが得られた。\ この2つのパラメーターの事後分布は非常に高い相関を持ち、blとbrのもっともらしい値はすべて狭い尾根に沿って横たわっている。blが大きければ、brも小さくなるはずである。ここで起こったことは、両方の脚の変数がほとんど同じ情報を含んでいるので、もしあなたがモデルに両方を含めることにこだわるなら、同じ予測を生み出すblとbrの組み合わせが実質的に無限に存在することになる、ということなのである。

この現象を、このモデルに近似させたと考えるのも一つの方法である

 y_{i} \sim \operatorname{Normal}(\mu i, \sigma)

 \mu_{i} = \alpha + \beta_1 x_i  +  \beta_2 x_i

変数yは例の身長のような結果で、xは例の脚の長さのような1つの予測変数であるここで,x は2回使われており,これは両方の脚の長さを使うことによって起こる問題の完璧な例である.ゴーレムの視点からは,μi のモデルは次のようになる.

 \mu_{i} = \alpha+ \left(\beta_1 + \beta_2 \right) x_i

私が行ったのは、各項からxiを取り除くことだけである。パラメータβ1 とβ2は別々に平均値µに影響を与えることはないので、分離することはできない。

このシミュレーションの例では,事後分布がまさにそうなっている以下は、それらの和の事後分布を計算し、それをプロットする方法である。

sum_blbr <- post$bl + post$br
dens( sum_blbr , col=rangi2 , lwd=2 , xlab="sum of bl and br" )

結果の密度プロットが図6.2の右側に示されている。事後平均は2強と右肩上がりで、標準偏差は和の成分であるblとbrのどちらよりもずっと小さくなっている。脚長変数の片方だけで回帰を当てはめても、ほぼ同じような事後平均が得られる。

m6.2 <- quap(
alist(
height ~ dnorm( mu , sigma ) , mu <- a + bl*leg_left,
a ~ dnorm( 10 , 100 ) , bl ~ dnorm( 2 , 10 ) , sigma ~ dexp( 1 )
) , data=d )
precis(m6.2)
      mean   sd 5.5% 94.5%
a     1.00 0.28 0.54  1.45
bl    1.99 0.06 1.89  2.09
sigma 0.62 0.04 0.55  0.69

その1.99はsum_blbrの平均値とほとんど同じである。
基本的な教訓は、これだけです。2つの予測変数が非常に強く相関しているとき(モデル中の他の変数の条件付きで)、モデルに両方を含めると混乱が生じるかもしれないこのような場合、事後分布は間違っていない。これらのデータでは事後分布は間違っていないのである。そして、「答えられない」と言うことは、モデルにとって素晴らしいことなのである。そして、もしあなたが予測に興味があるだけなら、この足のモデルは素晴らしい予測をすることがわかる。ただ、どちらの脚がより重要であるかは主張しないのである。

この脚の例は明確でキュートである。しかし、これは純粋に統計的なものでもある。私たちはここで深刻な因果関係を問うているのではない。次はもっと因果関係のある面白い例で試してみましょう。

多重共線性 ミルク

脚の長さの例では、両脚をモデルに含めることが少し愚かであることは容易に理解できる。しかし、実際のデータセットで生じる問題は、相関の高い予測変数の間の衝突を予期していないかもしれないということである。そして、それゆえ、我々は、どちらの予測変数も重要でないと言うように事後分布を誤って読んでしまうかもしれない。このセクションでは、実際のデータでこの問題の例を見る。

先ほどの霊長類の乳のデータに戻ろう。

library(rethinking)
data(milk)
d <- milk
d$K <- standardize( d$kcal.per.g )
d$F <- standardize( d$perc.fat )
d$L <- standardize( d$perc.lactose )

この例では、perc.fat (脂肪率) と perc.lactose (乳糖率) という変数に注目するこれらを使って、総エネルギー量である kcal.per.g をモデル化する。上のコードでは、これら3つの変数をすでに標準化している。あなたは、多重共線性の自然なケースを探索するために、これらの3つの変数を使用することになる。これらの列には、NAという欠損値がないので、ここでは完全なケースを抽出する必要はないことに注意してほしい。しかし、quapはlmのような無謀な関数とは異なり、黙ってケースを落とすようなことは決してありませんので安心してほしい。
まず、kcal.per.g を perc.fat と perc.lactose の関数として、2つの二変量回帰でモデル化することから始めまよう。これらの事前分布については第5章(147ページ)を参照。

# 6.9 kcal.per.g regressed on perc.fat
m6.3 <-quap(
  alist(
    K ~dnorm(mu,sigma),
    mu <-a+bF*F,
    a ~dnorm(0,0.2),
    bF ~dnorm(0,0.5),
    sigma ~dexp(1)
  ) ,data=d)

# kcal.per.g regressed on perc.lactose
m6.4 <-quap(
  alist(
    K ~dnorm(mu,sigma),
    mu <-a+bL*L,
    a ~dnorm(0,0.2),
    bL ~dnorm(0,0.5),
    sigma ~dexp(1)
  ) ,data=d)

precis( m6.3)
precis( m6.4)
      mean   sd  5.5% 94.5%
a     0.00 0.08 -0.12  0.12
bF    0.86 0.08  0.73  1.00
sigma 0.45 0.06  0.36  0.54


       mean   sd  5.5% 94.5%
a      0.00 0.07 -0.11  0.11
bL    -0.90 0.07 -1.02 -0.79
sigma  0.38 0.05  0.30  0.46

bF と bL の事後分布は、基本的に互いに鏡像である。bFの事後平均は、bLの平均が負であるのと同様に正である。両方とも、ほとんどゼロの片側にある狭い事後分布である。それぞれの予測因子と結果の強い関連性を考えると、両方の変数が、種を問わず、乳中の総エネルギーの信頼できる予測因子であると結論づけられるかもしれない。脂肪が多いほど、牛乳のキロカロリーは多くなる。 乳糖が多いほど、牛乳のキロカロリーは少なくなる。しかし、両方の予測変数を同じ回帰モデルに入れたらどうなるかを見てみよう。

m6.5 <- quap(
  alist(
    K ~ dnorm( mu , sigma ) , mu <- a + bF*F + bL*L ,
    a ~ dnorm( 0 , 0.2 ) ,
    bF ~dnorm(0,0.5),
    bL ~dnorm(0,0.5),
    sigma ~dexp(1)
  ) ,
data=d )
precis( m6.5)
       mean   sd  5.5% 94.5%
a      0.00 0.07 -0.11  0.11
bF     0.24 0.18 -0.05  0.54
bL    -0.68 0.18 -0.97 -0.38
sigma  0.38 0.05  0.30  0.46

このとき、bFとbLの事後平均は0に近くなっている。また、両パラメータの標準偏差は二変量モデル(m6.3、m6.4)のときの2倍になっている。 これは脚の例と同じ統計的現象である。何が起こったかというと、変数perc.fatとperc.lactoseが同じ情報を多く含んでいることである。これらは、ほとんどお互いに代替しているのである。その結果、両方を回帰分析に含めた場合、事後分布はbFとbLの組み合わせの長い尾根を記述することになり、等しく妥当であることがわかる。 脂肪と乳糖の場合、この2つの変数は基本的に1つの変動軸を形成している。これを見る最も簡単な方法は、ペアプロットを使うことである。

pairs( ~kcal.per.g+perc.fat+perc.lactose,data=d,col=rangi2)

このプロットを図 6.3 に示す。対角線に沿って、変数がラベル付けされている。対角線から外れた各散布図では、縦軸の変数は同じ行にラベル付けされた変数で、横軸の変数は同じ列にラベル付けされた変数である。 例えば、図6.3の最初の行の2つの散布図は、kcal.per.g(垂直)対perc.fat(水平)、そしてkcal.per.g(垂直)対perc.lactose(水平)である。 脂肪の割合は結果と正の相関があり、乳糖の割合は結果と負の相関があることに注意。次に、中央の列の一番右の散布図を見てほしい。このプロットは、乳糖(水平)に対する脂肪分(垂直)の散布図です。点がほぼ直線上に並んでいることに注目してください。 この2つの変数には負の相関があり、非常に強い相関があるため、ほとんど冗長になっています。どちらかがkcal.per.gを予測するのに役立ちますが、もう一方が分かっていれば、どちらもそれほど役に立たない。

科学文献の中には、多重共線性に対処するための様々な方法が紹介されている。その中で因果関係を考慮したものはほとんどない。実際に、モデルを当てはめる前に一対の相関を調べ、相関の高い予測因子を特定して削除するよう学生に教えている分野もある。これは間違いである。一対相関が問題なのではない。相関関係ではなく、条件付きの関連性が問題なのである。 そして、その場合でも、何をするのが正しいかは、何が共線性を引き起こしているかによる。データ内の関連性だけでは、何をすべきかを決めることはできない。
ミルクの例で起こっていることは、哺乳類の母親が従わなければならない、ミルクに関するトレードオフが存在することです。ある種が頻繁に授乳する場合、そのミルクは水っぽく、エネルギーが低い傾向がある。このようなミルクには糖分(ラクトース)が多く含まれている。そのかわり、授乳回数が少ない種は、母乳は高エネルギーである必要がある。そのようなミルクは脂肪分が非常に高い。これは、次のような因果関係を意味する。

中心的なトレードオフは、ミルクの濃さ(D)を決めることである。この変数は観測されていないので、丸で囲んである次に脂肪Fと乳糖Lが決定される。最後に、F と L の組み合わせでキロカロリー K が決まる。もし D を測定できたり、種の他の側面に基づいて Dを予測する進化的・経済的モデルがあれば、回帰分析でつまずくよりましだろう。
多重共線性の問題は、モデルの適合に関する問題の一つで、非同定性(non-identifiability)と呼ばれるものである。あるパラメータが識別不可能な場合、データとモデルの構造上、そのパラメータの値を推定することができないことを意味する。 この問題はモデルのコーディングの間違いに起因することもあるが、多くの重要なタイプのモデルは、たとえ完全に正しくコーディングされていても、識別不可能な、あるいは弱く識別可能なパラメータを有している。自然は、たとえモデルが正しくても、我々に簡単に推論する義務を負ってはいないのである。
一般的に、利用可能なデータが目的のパラメータについて多くの情報を含んでいるという保証はない。 そのような場合、ベイズ計算機は事前分布に非常に近い事後分布を返す。 したがって、事後分布を事前分布と比較することは、モデルがデータからどの程度の情報を抽出したかを確認する良い方法となる。事後分布と事前分布が似ているということは、計算が間違っているということではない。しかし、もっと良い質問をするきっかけになるかもしれない。

Rethinking: 識別は保証する。理解はあなた次第。技術的に言えば、ベイズモデルにおいて同定性は問題ではない。なぜなら、事後分布が適切である限り、つまり、1へ積分する限り、すべてのパラメータが同定されるからである。しかし、この技術的な事実は、事後分布の意味を理解できることを意味するものではない。だから、ベイズ的な文脈では、パラメータの同定が弱いという言い方をした方がいいかもしれない。しかし、この違いは技術的なものに過ぎないかもしれない。実は、DAGで因果関係が特定できるはずだと言われても、統計的に特定できない場合がある。私たちは、デザインを考えるのと同じように、統計にも力を入れなければならないのである。

Overthinking: 共線性のシミュレーション。2つの予測変数の間の関連性によって事後分布の不正確さがどのように増加するかを見るために、シミュレーションを使おう。以下のコードは、相関のある予測変数を生成し、モデルを当てはめ、perc.fatとkcal.per.gを関連付ける傾きの事後分布の標準偏差を返す関数を作成する。

library(rethinking)
data(milk)
d <- milk
sim.coll <- function( r=0.9 ) {
  d$x <- rnorm( nrow(d) , mean=r*d$perc.fat , 
    sd=sqrt( (1-r^2)*var(d$perc.fat) ) )
  m <- lm( kcal.per.g ~ perc.fat + x , data=d )
  sqrt( diag( vcov(m) ) )[2] # stddev of parameter
}
rep.sim.coll <- function( r=0.9 , n=100 ) {
  stddev <- replicate( n , sim.coll(r) ) 
  mean(stddev)
}

r.seq <- seq(from=0,to=0.99,by=0.01)
stddev <- sapply( r.seq , function(z) rep.sim.coll(r=z,n=100) )
plot( stddev ~ r.seq , type="l" , col=rangi2, lwd=2 , xlab="correlation" )

r.seqの各相関値に対して、このコードは100回の回帰を生成し、それらからの平均標準偏差を返す。このコードでは、暗黙の平行事前分布を使用しているが、これは悪い事前分布である。そのため、共線性変数の効果を誇張してしまう。情報量豊富な事前分布を用いると、標準偏差の増大はより緩やかになる。