こんにちは。林岳彦です。統計的因果推論の本の初稿を書き上げるまで髪の毛を切らないぞ、と願掛けしましたが、書けないままどんどん髪の毛だけが伸びてきています。いつか塔に籠もってラプンツェルになるかもしれません。あるいは突然全てが嫌になって前田大然になるかです。今日は久々のリアル外勤のスキマ時間でエイヤッとこの記事を書いています。
さて。
今日は、ややマニアックな話として、重回帰分析が突然バグる状況について書きたいと思います。結論から言うと、重要な特性において分布の偏りが大きい変数で調整するときに、調整によって回帰が突然バグる場合があるので注意しましょうという話です。
例を見ていこう
例として、ある環境汚染物質が健康被害を引き起こしている例を考えます。ここでは、それぞれの人の汚染物質への曝露量と、健康影響の程度（バイオマーカーの値で測定）のデータが得られているとします。
以下の話では、「環境曝露によって健康影響（バイオマーカーの値）の値が増加するという真の因果関係がある」という状況を考えていきます。ここで、曝露量とバイオマーカーの値の関係は、以下のシグモイド曲線で規定されているとします。（尚、毒性学的には曝露量と反応がシグモイド型になることはとても一般的な想定です）
具体的データとして、以下の散布図の例を考えていきます。このデータはRでシグモイド関数から生成した値に誤差を足して作成したものであり、ここで見られている相関そのものはシグモイド型の「曝露量→バイオマーカーの値」の用量反応関係を反映しています。ただし、誤差もまあまあ大きいため、全体をぱっと見すると線形の比例関係にも見えています。
このデータをRで線形単回帰（「バイオマーカー」が目的変数、「曝露量」が説明変数）をしてみると、以下のような解析結果になります：
## Call:
## glm(formula = BiomarkerScore ~ scale(ExposureDose), data = dataAB.df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -13.4270 -2.7530 0.0925 3.5405 9.1241
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.5716 0.7639 -5.985 5.98e-07 ***
## scale(ExposureDose) 3.6659 0.7736 4.739 2.99e-05 ***
散布図から容易に予想されることですが、曝露量（ExposureDose）との強い関連が推定されています（偏回帰係数=3.6659, p<0.0001）*1。結果の解釈としては、「曝露量と健康影響には強い関連がある」という結論になります。元々シグモイド型のデータに線形回帰をしているので、予めモデルは誤設定されていることになりますが、概ね正しい結論とメッセージが得られていることになります。線形モデルでも実用上はまずまずの推定はできるよね、というかんじです。
さて次は、データの値そのものは全く変えずに、設定を少し変えてみます。今回のデータは、異なる２つの地域（A村とB村）から得られたデータを含んでいる状況になっています。散布図は以下になります：
今回のポイントとして、この例ではA村とB村で曝露量の分布は重なっていることを覚えておいていただければと思います。さて、ここで説明変数として「曝露量」と「地域」を考慮した線形重回帰分析をしてみると、結果は以下のようになります：
## glm(formula = BiomarkerScore ~ scale(ExposureDose) + Area, data = dataAB.df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -13.1285 -2.7341 0.0121 3.3540 9.4185
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.2757 1.0931 -3.912 0.000378 ***
## scale(ExposureDose) 3.6783 0.7831 4.697 3.58e-05 ***
## AreaB -0.5918 1.5466 -0.383 0.704176
ここでは「曝露量（ExposureDose）」については前回の解析とほぼ同じような強い関連が推定されています（偏回帰係数:3.6783, p<0.0001）。一方、「地域（Area）」については特に関連がないと推定されました（偏回帰係数:-0.5918, p=0.701476）。全体的な結果の解釈は前回と同じで、「曝露量と健康影響には強い関連がある」という結論になります。穏当です。
さてここからが本番です。
データの値そのものは全く変えずに、設定をさらに変えます。今回は、異なる２つの地域（A村とB村）で曝露量の分布が全く重なっていない状況を考えます。たとえば、A村の中に汚染物質の特異的な発生源があるような状況ではこのようなパターンのデータになることが実際にありえます。散布図はこうです：
多重共線性も気になるので書いておきますが、この例では「地域（村）」と「曝露量」の相関係数は0.84くらいです。
ここで前回と全く同様に、説明変数として「曝露量」と「地域」を考慮した線形重回帰分析をしてみると、結果は以下のようになります：
## glm(formula = BiomarkerScore ~ scale(ExposureDose) + Area1, data = dataAB_mod1.df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.519 -2.712 1.443 3.054 5.755
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1971 1.1638 -0.169 0.866
## scale(ExposureDose) -0.6486 1.1568 -0.561 0.578
## Area1B -10.2929 2.3106 -4.455 7.5e-05 ***
はい！「曝露量（ExposureDose）」との関連は完全に消えてしまいました（偏回帰係数の符号は負に変化:-0.6489, p=0.578）。その代わりに、今回は「地域（Area1B）」との強い関連が推定されています（偏回帰係数:-10.2929 p<0.0001）。 念のため繰り返しますが、曝露量とバイオマーカーの値は前回までのデータと一切変わっていません。
この推定結果をそのまま解釈すると、「A村で起きている健康被害にはこの汚染物質の曝露は関係ないんだ！」「A村自体が問題なんだ、汚染物質とは別の原因があるんだ！」となります。これでいいはずがありません。このデータはたとえばA村内に特異的な発生源があるためにこうしたパターンになっているだけで、実際には健康被害は汚染物質の曝露で生じているわけです。この解析結果は、健康被害の原因を見つけるという観点からは、定量的にも定性的にも完全に間違った推定結果になっています。 こうした結果は、A村での環境汚染に対策における不作為を正当化しかねないような、かなり由々しき解析結果となってしまっています。
ついでとなりますが、設定を少し緩めて、部分的に分布が重なっている状況で全く同じように線形重回帰解析をした例は以下のようになります：
# glm(formula = BiomarkerScore ~ scale(ExposureDose) + Area2, data = dataAB_mod2.df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -10.9165 -3.0670 0.7647 3.4124 9.8028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.405 1.252 -1.921 0.0625 .
## scale(ExposureDose) 1.659 1.198 1.385 0.1744
## Area2B -5.098 2.393 -2.131 0.0398 *
分布に完全に重なりがないときほどひどくはありませんが、あいかわらず「曝露量（ExposureDose）」との関連は不明瞭であり（偏回帰係数:1.659, p=0.1744）、「地域（Area2B）」には有意差がついている（偏回帰係数-5.098, p=0.0398）という、いまだミスリードを起こしうる、なんともぼんやりとした推定結果になってしまっています。
まとめとして上記の現象を一枚のポンチ絵にしてみると、以下のようになります：
なんでこんな推定結果になるの？
隙間時間の中でここまで書いてきて、かなりお腹が空いてきてしまい、説明するのがだんだんと面倒くさくなってきていますが（私よりうまく説明できそうな方はぜひtwitterなどでご解説いただければと思います）、要因を簡単にまとめると以下のようになります：
- 「非線形応答に線形モデルを当てはめている」かつ「異なるカテゴリー（A村とB村）で曝露量の分布の重なりがない（小さい）」
この両方の状況が同時に生じると、下図のように、「各カテゴリー内での線形フィットのトレンド」と「全体でみたときの関係のトレンド」の乖離が大きくなってしまうことにより、関連の大きさの推定を大きく狂わせるバグ的な状況が成立してしまうことがあるわけです。
これは単回帰では大きな問題が無かったのに、調整のために変数を加えたことでかえって問題が生じるという、けっこう回避難度の高いバグかなあと思います。こうした状況はレアケースであるとは思いますが、推定が量的・質的に完全に狂ってしまう場合もありえるため、誤りの帰結は重大になりうるところが怖いところです。
ここまで読んできて、「調整のための変数の追加により生じるバグ」であることから、いわゆる合流点バイアスを思い浮かべた方もいるかもしれません。しかし、上記の話は本質的にはモデルの誤設定が根本にある話であり、（いわゆる）合流点バイアスとは少し違う話*2で、たとえば、適切なモデル（何らかの制約を課したシグモイド型のモデルなど）で曝露量と反応の関係を回帰していれば、こうした調整によるバグは生じない可能性があります*3。
そもそもなんでこんなこと書こうと思ったの？
『欧州の排外主義とナショナリズム』という本を読んでいたら、その中で「収入の多寡は移民忌避的な態度の原因でない」という解析結果が示されていて、でももしかしたら下図のようなケースもありうるのかなあ？とちょこっと思ったからです。
本当にこういうケースが生じていると思っているわけではないですが、ちょっと思ってしまったついでにこの機会にテクニカルな記事として話をまとめてみたという次第です。
最後に伝えたいこと
解析ソフトにかける前に色々なパターンで散布図を描いて眺めてみてデータに病的な部分（事前に想定していなかった極端な偏りの傾向とか）がないかを確かめるという作業は、地味ですが、本当に本当に大事な作業なので、必ずやりましょう！
以上、いかがでしたでしょうか。（空腹のため書く気力がもう切れました*4）
ではまたどこかで
引用文献
![[モンベル] mont-bell トレールワレット BK（ブラック）1133248](https://m.media-amazon.com/images/I/518fCBPUClL._SL500_.jpg "[モンベル] mont-bell トレールワレット BK（ブラック）1133248")
