どもです。林岳彦＆オメガトライブです。きみは1005%（消費税込）
さて。
今回は、前回の記事：
今回は因果関係があるのに相関関係が見られない4つのケースをまとめてみた（前編：検定力が低い） - Take a Risk: 林岳彦の研究メモ
のつづきの”中編”になります。本記事では「因果関係があるのに相関関係が見られないケース」の中でも、「交絡・合流点」が関わるケースについて書いていきます*1。
扱う内容の範囲としては、最初の記事：
因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） - Take a Risk: 林岳彦の研究メモ
と重複する部分がかなりありますが、今回の記事では、「仮想例のデータ生成」の段階からRでの計算を交えて説明していきたいと思います。（今回はちょっと「R実習」のような趣になるので、Rの読み書きができないと分かりにくい部分が多々あるかもしれませんが、申し訳ありません）
（最初に断っておきますが、今回もけっこう長いです。本当にすみません）
では書いていきたいと思います。
因果関係があるのに相関が見られないケース（２）：交絡によるバイアスにより打ち消されてる
まずは、「交絡によるバイアスにより、因果効果が見かけの上で打ち消されている」ケースについて見ていきます。
仮想例として、以下のような因果構造のもとでの「河川において農薬が生物種数に与えている影響」を考えていきます：
ここで、「農薬濃度」は調査地点（5000地点*2）における河川中の農薬濃度を、「生物種数」は各地点での河川中の生物（動物プランクトン）の種数を、「近傍の水田面積」は各調査地点の近傍の水田面積を表すとします。
上図の因果構造は、各地点において「河川中の農薬濃度の増加は河川中の生物種数を減少させる」一方で、「近傍の水田面積」の増加は「農薬濃度」と「生物種数」をともに増加させる*3効果を持つことを意味しています。
はい。
では、このような因果構造のケースに対して回帰分析を行うとどうなるか、Rで解析を行っていきましょう（以下のRのコードファイルはConfounding_bias.R 
）。
今回は、そもそもの仮想データを生成するところから始めたいと思います。まずは、近傍の「水田面積」のデータ（サンプルサイズn=5000地点）を以下のように生成します：
n.sample <- 5000
SuidenMenseki_scaled.v <- rnorm(n.sample,mean=0,sd=1)
ここでは、正規分布からの乱数として5000個のデータが生成され、"SuidenMenseki_scaled.v"に格納されています。解析の単純化のために、データの数値は「平均0, 標準偏差1」となるように既に標準化されている状況を想定しています。（今後、本記事ではデータの数値は基本的に標準化されているものとして考えていきますので、面積や濃度の値が「マイナス」になっていたりしますが、そこは本記事の本質的な議論とは関係してこないので余り気にしないでください）
次は、この「水田面積」のデータから、対応する「農薬濃度」のデータを生成します：
error1.v <- rnorm(n.sample,mean=0,sd=1)
NouyakuNoudo.v <- 1.0*SuidenMenseki_scaled.v + error1.v
NouyakuNoudo_scaled.v <- scale(NouyakuNoudo.v)
ここでは「農薬濃度＝1.0×水田面積＋誤差」という関係式をもとに「農薬濃度」のデータが作りだされています*4。
さらに、この「農薬濃度」と「水田面積」から、対応する「生物種数」のデータを生成します：
error2.v <- rnorm(n.sample,mean=0,sd=1)
SeibutsuSyusuu.v <- 0.5*SuidenMenseki_scaled.v - 0.3*NouyakuNoudo_scaled.v + error2.v
ここでは「生物種数＝0.5×水田面積 - 0.3×農薬濃度＋誤差」という関係式をもとに「生物種数」のデータが作りだされています。
この関係式は、「農薬濃度を1増加させる→生物種数は0.3減少する」ことを意味しています。つまり、ここで「農薬濃度→生物種数」の因果関係の存在が規定されていることになります。ここで、「農薬濃度」が「生物種数」に対してマイナスの影響（係数-0.3）を与えていることを覚えておいてください。
最後に、各データをデータフレームにまとめて最初の数行の中身をチラ見してみましょう：
>
> syusuu_nouyaku_suiden.df <- data.frame(SYUSUU = SeibutsuSyusuu.v, NOUYAKU = NouyakuNoudo_scaled.v, SUIDEN = SuidenMenseki_scaled.v)
>
> head(syusuu_nouyaku_suiden.df)
SYUSUU NOUYAKU SUIDEN
1 1.2260034 1.45111399 1.0752558
2 1.3868468 0.05863758 0.6772032
3 0.7522993 1.00719354 2.2926979
4 -0.4177403 0.76182684 1.9777108
5 0.6866292 0.13305951 -0.3588315
6 0.6610084 -0.46573089 0.4718610
それぞれの列は生物種数（SYUSUU）、農薬濃度（NOUYAKU）、水田面積（SUIDEN）のデータの値を表しています。それぞれの行は各調査地点を表しており、このデータは全部で5000行（=5000地点）のデータを含んでいます。
さて。これでデータは準備できました。（おつかれさまでした）
ではまず、データの様子を直感的に理解するために、ペアプロットを描いてみます*5：
pairs(syusuu_nouyaku_suiden.df[1:500,], col="green", pch=20, ps=0.1)
これらの散布図を見ると、「水田面積」と「農薬濃度」の間には強い相関が見られることが分かります。しかし、そもそも「生物種数＝0.5×水田面積 - 0.3×農薬濃度＋誤差」という関係式をもとに生物種数のデータを生成したにもかかわらず、「生物種数」はそのどちらとも相関が見られません。
これはどういうことでしょうか？
ここで、回帰分析を適用してみましょう。
もともとの目的は「河川において農薬が生物種数に与えている影響」を見ることでしたので、まずは「生物種数」を目的変数、「農薬濃度」を説明変数にして単回帰分析をしてみます（モデル式：生物種数＝A×農薬濃度＋B）：
>
> lm_syusuu_by_nouyaku <- lm(SYUSUU ~ NOUYAKU, data=syusuu_nouyaku_suiden.df)
> summary(lm_syusuu_by_nouyaku)
（中略...）
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.02252 0.02887 0.780 0.435
NOUYAKU 0.03156 0.02887 1.093 0.274
上記の単回帰の結果からは、農薬濃度（NOUYAKU）と生物種数（SYUSUU）の関連は小さく（回帰係数が0.03）、p値も0.27であることが示されています。
この例においては、農薬濃度と生物種数の間に明白な因果関係がある（＝そもそも生物種数データ自体が「生物種数＝0.5×水田面積 - 0.3×農薬濃度＋誤差」という式を用いて生成されている）にもかかわらず、単回帰による解析では農薬濃度の回帰係数の値はそれよりも1ケタ小さく、しかも負号も逆になっています。
このような解析結果のみをもとに「因果効果」についての結論を出すと、「農薬濃度は生物種数に影響を与えない」あるいは「農薬濃度を低下させても生物種数は増加しない」という誤った結論を導きかねません。全くもって困った事態です。
（既にご承知の方々も多いと思いますが）この例の状況は「交絡要因」の存在により生じています。
交絡要因とは「説明変数と目的変数の両方に影響を与える上流側の要因」のことであり、今回の例では以下の因果グラフ（再掲）より、「農薬濃度→生物種数」の影響を見る上で「水田面積」が交絡要因になっていることが分かります（交絡についてのより詳しい解説はこちら）。：
この例では、「農薬濃度」と「生物種数」の間には因果関係*6が確かにあります（つまり、「農薬濃度」を変化させると「生物種数」も変化する）。
しかし、「農薬濃度が高い地点では近傍の水田面積も大きい」傾向があり、さらに「近傍の水田面積が大きいと生物種数が増える」ため、「農薬濃度→生物種数」の影響が、「水田面積→生物種数」の影響により見かけの上で打ち消されやすい（相関関係が見えにくい）形になっています。（実際問題として交絡の影響によりどの程度「打ち消される」かは、各影響の大きさのバランスによって決まります）
このような交絡の影響を除去するためには、説明変数として「交絡要因」となる変数を加えた重回帰分析を行うのが定石です。
実際に、説明変数に交絡要因である「水田面積（SUIDEN）」も加えた重回帰分析をしてみましょう（モデル式：生物種数＝A×農薬濃度＋B×水田面積 + C）：
>
> lm_syusuu_by_nouyaku_suiden <- lm(SYUSUU ~ NOUYAKU + SUIDEN, data=syusuu_nouyaku_suiden.df)
> summary(lm_syusuu_by_nouyaku_suiden)
（...中略...）
Coefficien
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.01912 0.02848 0.671 0.502
NOUYAKU -0.30682 0.04054 -7.568 4.5e-14 ***
SUIDEN 0.47934 0.04087 11.727 < 2e-16 ***
結果として得られた「農薬濃度」の偏回帰係数は「-0.31」となっています。これは、データ生成時に用いた関係式の係数（-0.3）と近い値になっており、「水田面積」を加えた重回帰により「農薬濃度→生物種数」の因果効果*7が適切に推定されていることが分かります。
この解析結果をもとに結論を出せば、「農薬濃度の増加は生物種数を減少させる」あるいは「農薬濃度を低下させると生物種数は増加する」という正しい（背景の因果構造を適切に反映した）結論に至ることができ、河川生態系保全のための適切な施策を立てることができそうです。
このように、「交絡要因の影響」については、交絡要因を重回帰の説明変数として加えることにより除去することができます*8。そのため、このような事例からの実践的なアドバイスとしてとりあえず「説明変数と関連のありそうな変数はモデルに入れておく」 というものが、ありうるのかもしれません（例えばここ）。
が。
だが。
しかし。
実は、とっても悩ましいことに、「説明変数と関連のありそうな変数」を加えたがために却ってバイアスが生じてしまう（その結果、偏回帰係数の解釈でやらかしてしまう）、というケースもあるのです。
以下、そんな２つのケースを見ていきます。
因果関係があるのに相関が見られないケース（３）：合流点バイアスにより打ち消されている
まずは、「合流点バイアス（collider bias）」のケースについて見ていきます。
合流点バイアスとは、データの収集時または解析時において、データが「因果の合流点において選抜／層別／調整されている」ために生じるバイアスです。
ここで"合流点（collider）"とは、例えば、以下の因果グラフでは、要因Bが「合流点」と呼ばれるものになります。
前々々回の記事では、データの"収集時"における合流点バイアス（選択バイアス）の例を検討しましたが、今回はデータの"解析時"における合流点バイアスについて見ていきたいと思います。
今回の仮想例は、「デュアスロンのタイム記録に対する自転車能力の影響」になります*9。
「デュアスロン」というのは、「トライアスロン（水泳・自転車・長距離走）」から水泳を除いたものに相当します（Wikipediaのデュアスロンの項）。今回の例では、「デュアスロン」のタイム記録に対して、以下のような因果構造を考えます。
ここでは、「長距離走・自転車・水泳」の個々の能力を示す値である「走力」「自転車力」「泳力」というものが独立に存在し、しかも何らかの形で測定されているとします（この辺りは現実的に考えるとちょっとアレですが、仮想例のための便宜として色々とご容赦ください）。
上の因果グラフは、「自転車力・走力が増加すると、デュアスロンのタイム記録は短縮される」という因果関係の存在を示しています。同様に、「自転車力・走力・泳力が増加すると、トライアスロンのタイム記録は短縮される」という因果関係も存在します。
はい。
果たして、このようなケースに対して回帰分析を行うとどうなるでしょうか？
今回もRでデータ生成の段階から解析を試行してみましょう（RのコードはCollider_duathron.R ）。
まずは、「自転車力」「走力」「泳力」のデータ（サンプルサイズn=5000人）を生成します：
n.sample <- 5000
bike_ability.v <- rnorm(n.sample,mean=0,sd=1)
run_ability.v <- rnorm(n.sample,mean=0,sd=1)
swim_ability.v <- rnorm(n.sample,mean=0,sd=1)
次に、「自転車力」と「走力」から、「デュアスロンのタイム記録」を生成します。
error1.v <- rnorm(n.sample,mean=0,sd=1)
duathlon_time.v <- -0.5*bike_ability.v - 1.5*run_ability.v + error1.v
ここでは、「デュアスロンタイム＝-0.5×自転車力 - 1.5×走力 ＋ 誤差」の関係式*10に基づき5000人分のデータを作成しています。
最後に、「自転車力」「走力」「泳力」から、「トライアスロンのタイム記録」を生成します。
error2.v <- rnorm(n.sample,mean=0,sd=1)
triathlon_time.v <- -1*bike_ability.v - 0.5*run_ability.v - 0.5*swim_ability.v + error2.v
ここでは、「トライアスロンタイム＝-1×自転車力 - 0.5×走力 - 0.5×泳力＋ 誤差」の関係式*11に基づき5000人分のデータを作成しています。
最後に、各データをデータフレームにまとめて最初の数行の中身をチラ見してみましょう：
>
> du_tri_data.df <- data.frame(BIKE_ABILITY = bike_ability.v, RUN_ABILITY = run_ability.v, SWIM_ABILITY = swim_ability.v, DUATHLON_TIME = duathlon_time.v, TRIATHLON_TIME = triathlon_time.v)
>
> head(du_tri_data.df)
BIKE_ABILITY RUN_ABILITY SWIM_ABILITY DUATHLON_TIME TRIATHLON_TIME
1 -0.37212056 -1.19281058 -1.37855174 1.1314993 -0.4740669
2 1.09813338 2.46914712 0.31185433 -3.8039360 -3.1358702
3 -0.05780920 -0.14392736 -0.80598730 -0.6573044 0.1752532
4 0.61386369 -0.56583193 0.03463905 0.7448045 -0.1604868
5 -0.67924502 -0.33796826 0.15632568 1.8537339 2.5391096
6 0.05139414 -0.04468604 -0.37751680 1.4944796 0.3739054
それぞれの列は自転車力（BIKE_ABILITY）、走力（RUN_ABILITY）、泳力（SWIM_ABILITY）、デュアスロンタイム記録（DUATHLON_TIME）、トライアスロンタイム記録（TRIATHLON_TIME）、のデータの値を表しています*12。
それぞれの行は各個人を表しており、このデータは全部で5000行（=5000人）のデータを含んでいます。
これでデータは準備完了です。（ふう）
データの様子を直感的に理解するために、まずはペアプロットを描いてみます*13：
pairs(du_tri_data.df[1:500,], col="pink", pch=20, ps=0.1)
これらの散布図を見ると、「自転車力」「走力」と「デュアスロンタイム」「トライアスロンタイム」との間には負の相関が見られます。一方、「自転車力」と「走力」の間には相関が見られません。これは各データの生成における仮定／過程と整合しており、何ら不思議のないものです。
はい。
では、回帰分析を適用していきます。
もともとの目的は「自転車力がデュアスロンタイムに与える影響」を見ることでしたので、まずは「デュアスロンタイム」を目的変数、「自転車力」を説明変数にして単回帰分析をしてみます（モデル式：デュアスロンタイム＝A×自転車力＋B）：
>
> lm_du_by_bike.res <- lm(DUATHLON_TIME ~ BIKE_ABILITY, data=du_tri_data.df)
> summary(lm_du_by_bike.res)
（...中略…）
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.005374 0.025155 0.214 0.831
BIKE_ABILITY -0.529530 0.025121 -21.079 <2e-16 ***
解析結果から、「自転車力」の回帰係数は「-0.53」と推定されており、これはデータ生成時の関係式における係数「-0.5」と近い値になっています。
つまり、ここでは「自転車力」のみを説明変数とした単回帰により、「自転車力→デュアスロンタイム」の因果効果は適切に推定されています。この例のような因果構造の場合には、単回帰が「必要にして充分」なのです*14。
では、このケースに対して重回帰分析を行ってみましょう。
ここでは、変数として「走力」が測定されておらず、その代わりに「トライアスロンタイム」のデータだけがある状況を考えてみます*15。
そして、とりあえず「説明変数（自転車力）と関連のありそうな変数はモデルに入れておく」という方針に従い、「トライアスロンタイム」を加えた重回帰分析を行います（モデル式：デュアスロンタイム＝A×自転車力 ＋ B×トライアスロンタイム ＋ C）：
>
> lm_du_by_bike_triathlon.res <- lm(DUATHLON_TIME ~ BIKE_ABILITY + TRIATHLON_TIME, data=du_tri_data.df)
> summary(lm_du_by_bike_triathlon.res)
（...中略…）
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.003051 0.024115 -0.127 0.899
BIKE_ABILITY 0.021665 0.030586 0.708 0.479
TRIATHLON_TIME 0.466000 0.019454 23.954 <2e-16 ***
この重回帰の結果からは、「自転車力」の偏回帰係数は「0.021」と推定されました。もともとのデータ生成時の関係式における係数「-0.5」とは値の大きさも負号も異なり、また有意差もすっかり消え失せています。
つまり、「トライアスロンタイム」を変数として加えたことにより、「自転車力→デュアスロンタイム」の因果関係はすっかり見えなくなってしまったわけです。。。
実は、この因果グラフの形では、合流点である「トライアスロンタイム」を加えることにより、「走力」が"交絡要因"へと変化してしまうのです*16。このように、合流点における変数を加える*17ことにより却って生じるバイアスってのが「合流点バイアス（collider bias）」というやつなのです。
さてさてさてさて。
もしここで、AICによるモデル選択を行うとどうなるでしょうか？
線形の重回帰で「自転車力」を含む各モデルのAICを計算すると以下のようになります：
最もAICが低いモデルは、説明変数に「自転車力」と「走力」を用いた重回帰モデルになっています。これは、このデータの生成におけるもともとの関係式（デュアスロンタイム＝ -0.5×自転車力 - 1.5×走力 ＋ 誤差）を適切に反映した妥当な結果と言えます。
ではここで、何らかの事情で「走力」のデータが利用できないようなケースを考えてみましょう。このときの最小AICのモデルは、「自転車力」と「トライアスロンタイム」を用いた重回帰モデルとなります。
ここで悩ましいのは、（上記の解析で示されているように）この「自転車力＋トライアスロンタイム」がもたらす「自転車力」の偏回帰係数（0.021）は、「自転車力」の本当の因果効果*18（関係式の係数＝-0.5）とはかなりかけ離れた値になっていることです。
このケースでは、「自転車力」のみを説明変数とした単回帰モデルにより適切に「自転車力」の因果効果を推定できるにもかかわらず、そのモデルよりも「自転車力＋トライアスロンタイム」のモデルの方がAICが低くなっているのです。
この例は、「モデルのAICが低いこと」と「そのモデルの偏回帰係数が因果効果をより適切に反映していること」は本質的には別モノの話であることを示しています*19。AICはあくまでもモデル全体による「予測能力」を最大化することを目指した指標であり、各変数の因果効果（介入による効果; cf. 過去記事）の推定を目指したものではないのです。
なので、一連のモデル群の中から「AIC最小のモデル」を見つけたとしても、その"ベストモデル"の偏回帰係数の値が「各変数の因果効果を最も適切に反映している」とナイーブに前提するのはやめておきましょう。（特に、解析のそもそもの目的が「介入」にある場合はなおさらです！→過去記事）
中編のまとめ
はい。今回の「中編」をまとめます。
今回は「因果関係があるのに相関関係が見られないケース」として：
- 交絡要因となる変数を考慮していないために、因果効果が見かけの上で打ち消されている（＝相関関係が見られない）
- 合流点にある変数を加えたがために、因果効果が見かけの上で打ち消されている（＝相関関係が見られない）
というケースについて見てきました。また、その中で：
- モデルにおける「AICが低いこと」と「各変数の偏回帰係数が因果効果を適切に反映していること」は本質的に別モノの話
ということにも触れました。
さて。
今回の合流点バイアスの件ですが、バイアスロン／トライアスロンの例のように変数間の質的な関係性（＝背景の因果構造）が分かっている場合には、合流点バイアスを避けることは比較的簡単かもしれません。しかし、背景の因果構造が不明な場合にそれを避けることはなかなかに難しいものです。
良かれと思って加えたもので逆に墓穴を掘る、というのは人生でもよくあることですよね。
一方、何かをやっているうちに「そもそもの目的がよく分からなくなってくる」というのも人生ではよくあることです。
次回の「後編」ではそのような例として、「中間変量によるマスク」のケースについて見ていきたいと思います。
参考文献：
今回の話をformalな形で理解したい方は以下の本をどうぞ：
既にある程度の知識のある方には以下の論文もオススメです：
構造的因果モデルについて
（単なる宣伝）
私が詩を寄稿させていただいた自主制作文芸誌が こちらやこちらで買えます（まだまだ赤字らしいのでどぞ4649です）。
.
![データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)](http://ecx.images-amazon.com/images/I/51jDJwFEw8L._SL160_.jpg "データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)")
