なぜ無作為化なのか：『因果推論の根本問題』とその解法

こんにちは。林岳彦です。はてなジョシュ（バーネット）です。今回から「はてなブログ」へ引っ越しました。今後とも引きつづきよろしくお願いします。

さて。

前回までの記事では、実験データではない調査観察データを用いた因果効果の推定における注意すべきバイアスの類型について書いてきました。

ここでなぜわざわざ「実験データではない」という但し書きをつけているのかというと、適切なデザインに基づき行われた実験（もしくは介入を伴う調査）からのデータは、処理・条件の違いによる結果の差を素直に「因果効果」とみなして解釈できるので、余り細かいことを考えなくても大丈夫だからです*1。

はい。

では、そもそも、なぜそのような実験では「結果の差を素直に因果効果とみなせる」のでしょうか？

今回は、その背景となるロジックについて書いていきたいと思います。

（すみません今回もものすごく長いです。。。）

まずは「因果効果」を定義してみる：「職業訓練が将来の年収に与える効果」を例に

はい。そもそも「因果効果」って何でしょうか。

論理学者でも無ければ抽象的なままに考えていくのはなかなか難しいので、今回は「職業訓練が将来の年収に与える影響」の例を通して考えていきたいと思います。

あなたの親友Kが失業しており、行政による職業訓練プログラムを受講するかどうか迷っているとします。このとき、判断の材料としては色々ありうるかとは思われますが、「その職業訓練が親友Kの将来の年収を増やすと期待できるのか」という点は重要といえるでしょう。

このとき、「職業訓練→親友Kの将来の年収」の”因果効果”を言葉で定義すると以下のように表現できるかと思います：

職業訓練を受講した場合の親友Kの5年後の年収　ー　 職業訓練を受講しない場合の親友Kの5年後の年収

職業訓練を受講した場合に収入が伸びればこの"因果効果"はプラスになりますし、職業訓練を受講した場合に収入が減ればこの"因果効果”はマイナスになります。（ここでは、「将来の年収」については便宜的に”5年後の"年収を指標とすることにします）

ちなみに、上の例をより一般的な形で表現すると、対象Xに対する「処理A→結果Z」の"因果効果"は：

対象Xに処理Aを行った場合の結果Z　ー　対象Xに処理Aを行わない場合の結果Z

と表すことができるでしょう。

ここまでは、感覚的にも理解しやすいかと思います。

しかしながら。

上記の「因果効果」は、本質的に”不可知"なのです。

「因果効果」 は不可知である：『因果推論の根本問題』

なぜ不可知なのか。

それは：

職業訓練を受講した場合の親友Kの5年後の年収　ー　 職業訓練を受講しない場合の親友Kの5年後の年収

の、両者の「場合」を共に観察することはできないからです。

つまり、「親友Kが職業訓練を受講した場合」には「職業訓練を受講しない場合の親友Kの5年後の年収」を知ることはできませんし、同様に「親友Kが職業訓練を受講しない場合」には「職業訓練を受講した場合の親友Kの5年後の年収」は知ることはできないのです。

身も蓋もない話ですが、まあ、そうですよね。

なので、上記の「因果効果」はそもそも知ることができないシロモノなのです。

「因果推論」において横たわるこの本質的な問題は『因果推論の根本問題（The fundamental problem of causal inference）』と呼ばれており、まあ基本的には途方に暮れるしかない問題です。

だがしかし。いつかはゆかし。むかしむかし統計学の世界にはロナルド・フィッシャーという途方も無い天才がおりまして、この『根本問題』に対する"銀の弾丸"を遺してくれました。それが、「無作為化」というアイデアです。

以下では、そのアイデアについて（すこし回り道をしながら）説明していきます。

さてどうするか：「集団への因果効果」を考えてみる

まずは、『因果推論の根本問題』への対処として、「集団への因果効果」を考えてみることにします。

統計学の基本的な考え方は「多くの事例を見れば法則が見えてくるかも（帰納法）」なので、"親友K”ひとりではなく、「職業訓練を受講した人／しない人」の多くの事例を集めることにより、「職業訓練→将来の年収」の「因果効果」に迫ることを試みてみましょう。

まず、調査への協力者をn人集めたとします。ここで、各個人i（n=1,…,n）への"因果効果"は以下のように定義できるでしょう：

職業訓練を受講した場合の個人iの5年後の年収　ー　 職業訓練を受講しない場合の個人iの5年後の年収

はい。これは、平たく言うと「（個人iにおける）職業訓練の受講の有無による5年後の年収の差」になります。残念ながら、この"因果効果"は上記の『根本問題』で見たとおり、知ることのできないシロモノです。

では、各個人ではなく、「n人の5年後の”平均”年収」に着目し「集団への”平均"因果効果」を以下のように定義してみます：

職業訓練を受講した場合のn人の5年後の平均年収　ー　 職業訓練を受講しない場合のn人の5年後の平均年収

はい。こちらは、n人の"集団"で見た時の「職業訓練の受講の有無による5年後の平均年収の差」になっています。

さて。

実は、ここで、「集団への"平均"因果効果」を考えてみても、『根本問題』はちっとも解決されていません。（すみません。。）

なぜなら、「n人が職業訓練を受講した場合」には「職業訓練を受講しない場合のn人の5年後の平均年収 」は知ることができませんし、「n人が職業訓練を受講しない場合」には、「職業訓練を受講した場合のn人の5年後の平均年収 」は知ることができないからです。

さて。もう少し考えを進めてみましょう。

全員について知る必要はない：n人を2つの「グループ」に分けてみる

次善の策を練ってみましょう。

「職業訓練を受講した場合のn人の5年後の平均年収」と「職業訓練を受講しない場合のn人の5年後の平均年収」を同時に知ることができないのなら、n人を「職業訓練を受講した人のグループ」と「職業訓練を受講しない人のグループ」に分けて考えれば良いのかもしれません。

ここでn人の内訳として、「職業訓練を受講した」グループがm人、「職業訓練を受講しない」グループがk人いたとしましょう（n=m+kとする）。

このとき、「職業訓練→将来の年収」の因果効果を：

職業訓練を受講したm人の5年後の平均年収　ー　 職業訓練を受講しないk人の5年後の平均年収

という形で捉えることはできるでしょうか？

うむ。

これは、できそうな気もします。少なくとも、この場合には第1項と第2項の両方の量を（原理的には）知ることができるという点では現実的と言えます。

しかしながら、これはこれで、問題が生じます。

なぜなら、「職業訓練を受講したm人」と「職業訓練を受講しないk人」の両者のグループには、かなり性質の違う人々が含まれている可能性があるからです。

例えば、「行政による職業訓練プログラムを受講した人のグループ」の中には、そもそもの流れとして現在失業中の人々が相対的に多く含まれているかもしれませんし、「受講しない人のグループ」の中には、現状のところは安定した職を得ている、もしくは専門的スキルが既にある人々が相対的に多く含まれているかもしれません。

このように両者のグループの性質が異なりうる場合には：

職業訓練を受講したm人の5年後の平均年収　ー　 職業訓練を受講しないk人の5年後の平均年収

という量を計算したところで、それが「職業訓練→将来の年収」の因果効果に起因するものなのか、「グループ間のそもそもの性質の違い」に起因するものなのか、うまく判別がつきません。

いたしかゆしです。

ではここで一旦、ちょっと理屈っぽい話となりますが、この状況を俯瞰してみることにします。

問題を俯瞰する：「反事実」の世界へようこそ

上記の状況を、Rubinの反事実モデルの枠組みを用いて整理してみます。

今まで見てきた職業訓練の例では、処理（条件）として「職業訓練を受講した／受講しない」の２つの場合、帰結として「職業訓練を受講した場合の5年後の年収／受講しない場合の5年後の年収」の２つの場合がありました。

これら「条件」と「帰結」の組み合わせを表にまとめると、論理的には、以下の4種類の「平均年収データ」がありうることになります：

上の図の4種類の「平均年収データ」について書き下すと：

• (I) 職業訓練を受講したグループの「職業訓練を受講した場合の5年後の平均年収」
• (II) 職業訓練を受講しないグループの「職業訓練を受講した場合の5年後の平均年収」
• (III) 職業訓練を受講したグループの「職業訓練を受講しない場合の5年後の平均年収」
• (IV) 職業訓練を受講しないグループの「職業訓練を受講しない場合の5年後の平均年収」

となります。

一見とてもややこしいですが、(I)と(IV)の「平均年収データ」については、特に問題はないかと思います。文として素直に意味が通ってますし、実際にデータを得ることも（少なくとも原理的には）可能です。

一方、(II)(III)の「平均年収データ」は、原理的に「そもそも観測することができないデータ」となっています。「職業訓練を受講したグループ」についての「職業訓練を受講しない場合の5年後の年収」はそもそも観測することができませんし、同様に「職業訓練を受講しないグループ」についての「職業訓練を受講した場合の5年後の年収」も観測することができません。

(II)(III)のようなデータは「事実に反した条件（counterfactual condition）」の下でしか得られないデータとなります。つまり、ここでは：

「職業訓練を受講しないグループ」が”もし職業訓練を受けていたとしたとき”の、5年後の平均年収」

という形で、「（事実としては受講していないけれども）もし職業訓練を受けていたとしたとき」という、反事実的な状況が想定されていることになります。

はい。

ここで重要なポイントを2つまとめます：

• (1) 「因果効果」の推論においては、必然的に「反事実的状況下のデータ」の値を知る必要がある
• (2) われわれが現実的に得られるのは、「事実として起きた状況」のもとでのデータのみである

この(1)については、「職業訓練→年収」の因果効果というのが「(I) と (III)の差」、つまり

　(I) 職業訓練を受講したグループの「職業訓練を受講した場合の5年後の平均年収」
　ー　 (III) 職業訓練を受講したグループの「職業訓練を受講しない場合の5年後の平均年収」

で定義されることに対応します*2。この後者の項は、反事実的状況下でのみ得られるデータとなっています。

(2)については、いわずもがなですが、われわれが得ることができるのは(I)と(IV)の状況下のデータのみである、ということです。

これらのポイントを図中にまとめると：

という状況であると言えます。

このような状況を前にして、"事実的世界"に存在する我々が目指しうるゴールは：

反事実的世界においてしか観測できない(III)のデータを、現実に観測可能な(I)や(IV)のデータを用いて”観測”する

ことになるわけです。

どうやるのかって？

ここで、"無作為化”の出番です。

その銀の弾丸を放て：無作為化により『根本問題』を解く

さて。

反事実的世界においてしか観測できない(III)のデータを、現実に観測可能な(I)や(IV)のデータを用いて”観測”する

ためにはどうすれば良いでしょうか。

ここで、「（III)と(IV)の値が一致すると期待できるような条件」を考えてみます。

この(III)と(IV)で異なる部分は、「受講したグループ」と「受講しないグループ」の部分です。そこで、最初に職業訓練を受けてもらう前の段階で、「両者のグループの間に差がでないようなやり方」でn人のグループ分けを行うことにしましょう。

ここで、「両者のグループに差がでないようなやり方」として最も汎用的な方法が、「無作為化」です。

例えば、 集まってもらったn人にコインを投げてもらい、オモテが出た人は「受講するグループ」になり職業訓練を受講してもらい、ウラが出た人は「受講しないグループ」になり職業訓練を受講しないでいてもらいます。（このように各個体に対して無作為に条件や処理のグループを割り付けることを「無作為割付」と言います）

このとき、「受講するグループ」と「受講しないグループ」の間には本質的な差は無いと考えられるので、理屈として、5年後の年収データとしての

• (IV) 職業訓練を受講しないグループの「職業訓練を受講しない場合の5年後の平均年収」

の期待値は

• (III) 職業訓練を受講したグループの「職業訓練を受講しない場合の5年後の平均年収」

の期待値と一致すると考えられます。なぜなら、両グループは単にコイントスで分けただけなので、（もし値が測定可能であった場合に）両者の値（の期待値）に差が出ると考えるべき理由が何もないからです。（例えば、両者のグループの「身長の平均値」の期待値が異なるとは考えられないですよね。それと理屈的には同じです。また、数式を用いた補遺を記事末尾に載せましたので適宜ご参照ください。）

もちろん実際には(III)のデータは反事実的なので得ることはできません。しかし、無作為化（コイントス）によりグループを分けた場合には、(III)と(IV)の「平均年収」の期待値は理屈的に一致すると考えられるので、ここで「(III)のデータの身代わりとして(IV)のデータを使う」という"迂回路"が使えるようになるわけです。

この"迂回路"を使うことにより、

　(I) 職業訓練を受講したグループの「職業訓練を受講した場合の5年後の平均年収」
　ー　 (III) 職業訓練を受講したグループの「職業訓練を受講しない場合の5年後の平均年収」

という不可知な項を含む"因果効果”を

　(I) 職業訓練を受講したグループの「職業訓練を受講した場合の5年後の平均年収」
　ー　 (IV) 職業訓練を受講しないグループの「職業訓練を受講しない場合の5年後の平均年収」

と読み替えることが可能になるわけです。これらの(I)と(IV)の量は観測可能なデータなので、この読み替えを行うことにより、「職業訓練→将来の年収」の"因果効果”を推定することができることになります。

はい。

これで、『因果推論の根本問題』を統計学的に解くことができました。（これが、いわゆる『統計学的因果推論』というものの一例となります）

（ただし、ここで推定できたのはあくまで集団への「"平均的"因果効果」であることに注意してください。そもそもの話である、職業訓練が「親友K」に与える因果効果は、やはり不可知なままなのです。たとえ職業訓練の"平均"効果がプラスであったとしても、親友Kにはそれがマイナスに働く可能性は常に残ります。個体に対する因果効果は、本質的にやはり不可知なのです。）

まとめ

ではまとめます。

今回は以下のようなことを見てきました：

• そもそも：因果推論を行うためには、「反事実的状況下におけるデータ」を観測する必要がある
• ゆえに：因果推論は不可能である（『因果推論の根本問題』）
• しかし：無作為化により、「反事実的データ」の身代わりに観測可能なデータを使うことができる（『根本問題』の統計学的解決）
• なかにし：「無作為化ハンパないって　反事実的なデータめっちゃトラップするもん」
• かんとく：「おれ握手してもらったぞ」

はい。

一見地味に感じるかもしれませんが、「ランダムネス」の意図的な利用により「法則」を浮かび上がらせるというこの「無作為化というアイデア」はなかなかに凄いものであるようにも思います。かのリチャード・ドーキンスが「ダーウィン以降で最も偉大な生物学者」としてフィッシャーの名を挙げていましたが、この「無作為化というアイデア」はその「ダーウィン以降で最も偉大な生物学者」が遺したものの中でも、「最も偉大な発明」の一つと言えるかもしれません。

そんなわけで、「統計的因果推論」の分野で有名な研究者といえば、RubinやPearlの名が挙がるのかもしれませんが、「統計的因果推論のキング」を一人選ぶとしたらやはりフィッシャーになるのではないかと思います。

はい。

というわけで、『統計学の道を歩むかぎり、お天道様とフィッシャーはどこへいってもついてまわる』という格言を再認識したところで、今回の記事も終わろうかと思います*3。

＃次回からは、そろそろ「確率概念」についてガッツリ書いていこうかなあ、と思ってます。

補遺(1)：無作為化のキモについて

上記の説明では、「コイントス」を用いて無作為化を行っていますが、別にそこは必ずしも「コイントス」でなくとも構いません。大事なことは、割付（＝各人をどのグループに分けるのか）を、「受講した場合／しない場合の5年後の平均年収」と本来的に全く関係ない（＝独立な）要素にもとづき決めるということです。

数式で書くと、上記の表の

• (III) 職業訓練を受講したグループの「職業訓練を受講しない場合の5年後の平均年収」
• (IV) 職業訓練を受講しないグループの「職業訓練を受講しない場合の3年後の平均年収」

の期待値は、それぞれ

• (III) E(職業訓練を受講しない場合の5年後の平均年収｜ 職業訓練を受講したグループに割付された）
• (IV) E(職業訓練を受講しない場合の5年後の平均年収｜ 職業訓練を受講しないグループに割付された）

と書けます。ここでE(A|B)という表記は、「Bという条件」の下での「Aの期待値」を表しています。

ここで、「B」の条件がAと全く関係ない（＝独立な）要素によって決められている場合には、E(A|B)=E(A)となるため：

• (III) E(職業訓練を受講しない場合の5年後の平均年収｜ 職業訓練を受講したグループに割付された）＝ E(職業訓練を受講しない場合の5年後の平均年収）
• (IV) E(職業訓練を受講しない場合の5年後の平均年収｜ 職業訓練を受講しないグループに割付された） ＝E(職業訓練を受講しない場合の5年後の平均年収）

となり、「(III)の期待値＝(IV)の期待値」となるわけです。

ここのキモは「Bの条件がAと全く関係ない（＝独立な）要素によって決められている」ことなので、コイントスやサイコロを使わずとも、例えば「各人の携帯電話番号の末尾が偶数か奇数か」にもとづき職業訓練受講の有無（B）の条件を割りつけても、E(A|B)=E(A)を満たすと考えられるので問題ありません。

逆に、職業訓練を受講していない場合の5年後の平均年収（A）と潜在的に関係のある要素に基づき割付を行うと、E(A|B)=E(A)が満たされないので因果推論がうまくいきません*4。例えば、「現在求職中かどうか」にもとづき職業訓練受講の有無（B）の条件を割りつけると、BとAは独立でない可能性が高く、E(A|B)=E(A)は満たされないでしょう。

また、細かい話になりますが、条件Bを無作為に割付けたとしても、「条件を割付けたこと自体」がAに影響を及ぼす場合は、E(A|B)=E(A)が満たされないのでやはりダメです。例えば、「訓練を受講するグループに割り付けられた」こと自体がモチベーションを変化させ、（実際の訓練受講とは関係ない部分で）年収に影響を与えてくる場合には、BとAは独立とならず、E(A|B)=E(A)は満たされません。このような場合はいわゆるSUTVA違反のケースに相当します*5。

補遺(2)：無作為化ができないときにはその悲しみをどうすりゃいいの？誰がぼくを救ってくれるの？

正直どうにもならないことも多いですが、とりあえず「傾向スコア」とか「バックドア基準」とかを参照してみてください。

統計的因果推論（傾向スコア）の勉強会資料をアプしてみた - Take a Risk：林岳彦の研究メモ
発表資料アプ『相関と因果について考える：統計的因果推論、その(不)可能性の中心』 - Take a Risk：林岳彦の研究メモ

===【20180806追記ここから】======================
久しぶりに読み返したところ、上の記述に対しては「サッパリすぎるだろ」と思ったので追記しておきます。無作為化が成立しているか怪しいときは「まずバックドアパスを閉じれるように工夫してみて、それでも閉じれているかどうか分からないときには想定されるバックドアパスの「太さ」について定量的に考察してみましょう（感度分析など）」となるかなあと思います。簡単にあきらめないで！

＊とりあえず以下の良記事を参照すると吉：
www.krsk-phs.com
www.krsk-phs.com
===【20180806追記ここまで】======================

詳しくは、以下の参考文献を読むのをオススメいたします。

【今回の参考文献】

今回の記事は以下の二つの名著の中の説明をmixtureしたものとなっております。ぜひこれらの書籍もご参照いただければと思います。

統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

【ついでに謝辞：今回のブログの引っ越しに際して参考にさせていただいた記事】

このブログのデザインは、はてな提供の"epic"のデザイン（背景を白に変更）を利用したものですが、「見出し」のデザインとサイドバーの「人気記事」の部分だけカスタマイズを加えています。

「見出し」のスタイル変更は以下のサイトを参考にさせていただきました：
【CSS初心者向け】ブログで良く見かけるシンプルな見出しの作り方 | delaymania
はてなブログのカスタマイズの方法についてまとめてみた【ま】 - はてブのまとめ

「人気記事」のサイドバーは以下の記事のものを利用させていただきました：
【修正済】コピペで簡単！ はてなブログの人気記事を画像付きで表示させる方法。 - #ChiroruLab

大変ありがとうございました。リスペクト＆感謝申し上げます。

（ついでのついでに単なる宣伝）

私が詩を寄稿させていただいた自主制作文芸誌が こちらやこちらで買えます。どうぞよろしくおねがいいたします。

.

*1:交互作用があると必ずしもそうでもない

*2:因果効果の定義としては特に「職業訓練を受講したグループ」に限定する必要はなく、同じグループにおける差異であればOKです。今回の場合は「職業訓練」の文脈上、「職業訓練を受講したグループ」に着目するのが自然と思われるのでこの形で書いています

*3:嗚呼、ここでwhat_a_dudeさんのフィッシャー関連記事を引用しようと思ったのにもう全て消されてるんだよね。。。what_a_dudeさんのいけず。。。

*4:バイアスがかかる

*5:「薬を処方するグループに割り付けられた」ことによるプラセボ効果がある場合などもこのケースになるかと思います

（後編）今回は因果関係があるのに相関関係が見られない4つのケースについてまとめてみた：中間変量の影響

どもっす。林岳彦です。先日、某所で統計解析の講師役をしました。その際に解析環境の準備の手間を省こうと思って、Amazon EC2上にRStudioのサーバー版を立てて、聴講者にそこに繋いでもらって実習をしようとしたのですが、いざ皆が繋いだらサーバーがクラッシュしまくって実習が全く進みませんでした*1 。。。すみませんでした(泣)*2。。

さて。

良かれと思ったもので逆に墓穴を掘る、というのは人生ではよくあることですよね！

前回の「合流点の追加によるバイアス」はそんな例の一つでしたが、今回の「後編」ではそのようなもう一つの例として、「中間変量の追加によるマスク」のケースについて見ていきます。

因果関係があるのに相関が見られないケース（４）：中間変量によってマスクされている

はい。では、中間変量によって因果効果がマスクされてしまうケースを見ていきます。

ここで「中間変量」というのは、「A→Z」のような因果関係がある場合に、「A→M→Z」における"M"のような因果の流れの中間に位置づけられるような変量のことを指します。分野によっては、媒介変数とか中間変数とか呼ばれることもあります。（本記事内での「因果」の定義については過去エントリーをご参照ください）

今回の記事では、「A→Z」の因果効果、つまり「Aを変化"させた"ときのZの変化量」に興味がある場合の、中間変量の影響について見ていきます。

今回は仮想例として、「塩分摂取量」が「死亡リスク（死亡率）」に与える影響の例を考えていきます。

具体的には、最も単純な形の一つとして、「塩分摂取量」「死亡リスク」「血圧」が以下のような因果構造を持っている場合を想定していきます（医学的には全く正確ではないかとは思いますが、説明のための便宜としてご容赦ください）：

上の図は、「塩分摂取量」が増加すると、「血圧」が増加し、そのこと（のみ）を通じて「死亡リスク」が上昇する因果構造を意味しています。

このような因果構造の場合に、回帰分析を行うとどうなるのか、Rで解析を試行してみましょう（以下のコードのファイルはintermediate.R ）。

今回も、そもそもの仮想データを生成するところから始めたいと思います。まずは、「塩分摂取量」のデータ（サンプルサイズn=10000人）を以下のように生成します：

# サンプルサイズを10000に設定
n.sample <- 10000
# 平均0、標準偏差1の正規分布から10000個の<span class="deco" style="font-weight:bold;">「塩分摂取量」</span>の値をランダムに生成
Enbun.v <- rnorm(n.sample,mean=0,sd=1)

ここでは、解析の単純化のために塩分摂取量データの数値は「平均0, 標準偏差1」となるよう既に標準化されているものと想定します。

次は、この「塩分摂取量」のデータから、対応する「血圧」のデータを生成します：

# 誤差項の作成
error1.v <- rnorm(n.sample,mean=0,sd=1)
# 「血圧＝5×塩分 ＋誤差＋平均値」の関係式により血圧データを作成
ketsuatsu.v <- 5*Enbun.v + error1.v + 70

ここでは「血圧＝5×塩分摂取量＋誤差＋平均値」という関係式をもとに血圧データが作りだされています。この関係式は、「塩分摂取量を1増加させる→血圧は5増加する」という因果関係を意味しています。

さらに、この「血圧」から、対応する「超過死亡リスク」のデータを生成します：

# 誤差項の作成
error2.v <- rnorm(n.sample,mean=0,sd=1)
# 「死亡リスク＝0.5×血圧 ＋誤差」の関係式により死亡リスクデータを作成
tyouka_shibou_caseA.v <- 0.5*ketsuatsu.v + error2.v

ここでは「超過死亡リスク＝0.5×血圧 ＋誤差」という関係式をもと超過死亡リスクデータが作りだされています。ここでは便宜的に、「超過死亡リスク」は「1万人あたりの死亡数の（ベースラインに対する）上昇分」の値であるとします。上記の関係式は、「血圧を1増加させる→死亡リスクは0.5増加する」という因果関係を意味します。

最後に、各データをデータフレームにまとめて中身をみると：

> # 生成したデータをデータフレームの形でまとめる
> Enbun_Shibou.df <- data.frame(ENBUN = Enbun.v, KETSUATSU = ketsuatsu.v, SHIBOU = tyouka_shibou_caseA.v)
> # データフレームの最初の数行の中身を見る
> head(Enbun_Shibou.df)
        ENBUN   KETSUATSU   SHIBOU
1 -0.19187249  70.56252       34.63921
2  0.23711386  72.01089       36.00139
3  1.65857552  79.01415       40.53588
4  0.01811544  70.96663       34.97887
5 -0.14860591  68.18043       33.54352
6  0.59624079  74.54500       36.02395

それぞれの列は塩分摂取量（ENBUN）、血圧（KETSUATSU）、超過死亡リスク（SHIBOU）のデータの値を表しています。それぞれの行は各個人を表しており、このデータは10000行（=10000人）のデータを含んでいます。

さて、これでデータは準備完了です。

ではまず、データの様子を直感的に理解するために、ペアプロットの散布図を描いてみます*3：

# 項目の対ごとに散布図を描く
> pairs(Enbun_Shibou.df[1:500,], col="red", pch=20, ps=0.1)

「塩分摂取量」「血圧」「超過死亡リスク」の全ての間に強い相関があることが分かります。（まあ、そういうふうにデータを作ったので当たり前です）

では、回帰分析を行ってみましょう。

もともとの目的は「塩分摂取量が死亡リスクに与える因果的影響」を見ることでしたので、まずは素直に「超過死亡リスク」を目的変数、「塩分摂取量」を説明変数にして単回帰分析（モデル式：超過死亡リスク＝ A×塩分摂取量＋B）をしてみます：

> # 目的変数「塩分摂取量」、説明変数「超過死亡リスク」の単回帰
> lm_tyouka_shibou_by_Enbun <- lm(tyouka_shibou_caseA.v ~ Enbun.v)
> summary(lm_tyouka_shibou_by_Enbun)
（...中略…）
Coefficients:
            Estimate Std.    Error        t value     Pr(>|t|)   
(Intercept) 34.99540       0.01105   3167.2     <2e-16 ***
Enbun.v      2.49680      0.01108    225.4      <2e-16 ***

解析結果から、「塩分摂取量（Enbun.v）」の回帰係数（『超過死亡リスク＝ A×塩分摂取量＋B』の"A"の値）は「2.49」と推定されています。ここは少し難しいかもですが、データ生成時に用いた関係式から「塩分摂取量1単位の増加→血圧の5単位増加」「血圧の5単位増加→釣果死亡リスクの2.5単位増加」という関係があり、これらをまとめると「塩分摂取量1単位の増加→釣果死亡リスクの2.5単位増加」となるので、上記の回帰式における回帰係数の理論値は「2.5」となります。で、ここの「2.49」という推定値は非常に理論値に近い値をとっており、「塩分摂取量→超過死亡リスク」の因果効果は適切に推定されていることが分かります。

今回の例のような因果構造の場合には、「塩分摂取量→超過死亡リスク」のシンプルな単回帰が、その因果効果の推定において十分な解析となっているわけです。

さて。

ここに「血圧」も変数に加えた重回帰分析をしてみましょう。「超過死亡リスク」を目的変数、説明変数として「塩分摂取量」と「血圧」の２変数を用います（モデル式：超過死亡リスク＝A×塩分摂取量＋B×血圧＋C）：

> lm_tyouka_shibou_by_Enbun_kestuatsu <- lm(tyouka_shibou_caseA.v ~ Enbun.v + ketsuatsu.v)
> summary(lm_tyouka_shibou_by_Enbun_kestuatsu)
（...中略…）
Coefficients:
                  Estimate   Std. Error     t value    Pr(>|t|)   
(Intercept)   1.069235   0.697701      1.533      0.125   
Enbun.v      0.075569   0.050775      1.488      0.137   
ketsuatsu.v 0.484687   0.009967    48.631      <2e-16 ***

はい。重回帰の結果を見ると、「超過死亡リスク」への「塩分摂取量（Enbun.v）」の偏回帰係数（『超過死亡リスク＝A×塩分摂取量＋B×血圧＋C』の”A”の値）は「0.076」になっています。単回帰での回帰係数「2.49」と比べると実質的にはゼロと言える値にまで小さくなっています。また、有意差も無くなっています（p=0.137）。

これは、そもそも重回帰の偏回帰係数というものが、「他の変数の影響を除去（＝他の変数は一定に固定）したとき」の「その説明変数が目的変数に与える影響の大きさ」を示しているからです。ある説明変数の影響が中間変量を介してのみ伝わる場合には、「他の変数を一定に固定した」ときにはその影響が「中間変数の固定」によりマスクされてしまうので、偏回帰係数はゼロになってしまうわけです。

さてさて。上記の重回帰の結果だけを見ると、あたかも「塩分摂取量」は「超過死亡リスク」に影響しないように見えます。はたして、この重回帰の結果から「塩分摂取量→超過死亡リスクの因果関係はない」と結論づけてしまって良いものでしょうか？

この辺りで、因果の「定義」の問題が絡んでくる*4ちょっとヤヤコシイ話になります。

今回の私の一連の記事では、「因果の定義」を：

「要因Aを変化させた（介入した）とき、要因Zも変化する」ときに、「要因A→要因Zの因果関係がある」と呼ぶ

としています。いわゆる「介入効果」を基にした「因果」の定義です（詳しくはここ）。このような「因果」の定義を採る場合には、この重回帰の結果から「塩分摂取量→超過死亡リスクの因果関係はない」と結論づけてしまうのは明らかに誤りと言えます。なぜなら、介入により「塩分摂取量」を1単位変化させると、（因果構造から明らかなように）結果的に「超過死亡率」も変化するからです。

一方、我々は日常的な用法においては、因果のメカニズム（特に、因果の直近の上流-下流関係）に基づいて「因果」という言葉を用いていることも多いように思います。

「因果」の日常的な用法の範囲においては、上記の重回帰の結果は、「塩分摂取量は血圧に対して直接的な因果関係を持ち」、「血圧は超過死亡リスクに対して直接的な因果関係を持つ」という言い方も成り立つように思います。このとき、「塩分摂取量は超過死亡リスクに対して（直接的な）因果効果はない」という言い方も（日常的な用法の範囲としては）おそらく成り立つのかもしれません*5。

因果関係の「有無」を巡っては、このように「介入効果に基づく定義」と「メカニズムに基づく定義」のどちらを採るかによって、言い方が変わってくることがあります。

もっとコントラストがわかりやすい例として、今までの図に「塩分摂取量」から「超過死亡リスク」への小さなマイナスの直接的影響の効果を追加した、以下のような因果グラフの場合を考えてみましょう：

この因果グラフにもとづき生成したデータを単回帰および重回帰で解析すると、結果として以下のような回帰式が得られたとします（コードはintermediate.R ）：

単回帰式：超過死亡リスク ＝ 2.3×塩分摂取量 + 定数
重回帰式：超過死亡リスク ＝ -0.12×塩分摂取量 ＋ 0.48×血圧 ＋ 定数

このとき、「塩分摂取量（の1単位増加）が死亡リスクに与える”因果的"影響」は、「介入効果」に基づく定義を採用すれば、単回帰の回帰係数によって示された「+2.3」となります*6。一方、「直近のメカニズム的関係に基づく因果の定義」を採用すれば、重回帰の偏回帰係数よって示された「-0.12」を因果的影響として採る場合もあるかもしれません*7。

このような状況において、ある人は端的に「塩分摂取量の増加は死亡リスクを増加させる」と言い、一方で、他の人は端的に「塩分摂取量の増加は死亡リスクを減少させる」という言い方をするかもしれません。

これは、どちらかが本質的に正しいという問題ではなく、（暗黙の前提として）「因果効果」をどう定義しているか、による違いといえるでしょう。

このような定義によるすれ違いを解消するためには、どのような意味で「因果効果」という語を用いているか、また、解析や調査デザインにおいてどのような変数で調整／層別されているのかを常に明確に意識することが大切です。また、可能であれば、上のような因果グラフを実際に描きながら議論をすると非常に捗るでしょう。

そして、このようなケースを前にして一番大事なことは、解析の途中で「そもそも何がしたいのか（介入効果が見たいのか／要因分析がしたいのか／予測したいのか／マイニングしたいのかetc）」を見失わないことです。

（まあ、「迷子になったがゆえにできる発見」というのも時折あるんですけどね...）

今回のまとめ

ではあっさりまとめます。

今回の内容をまとめると：

• 中間変量を加えることにより「因果効果」が見かけ上マスクされることがあるので注意
• 「因果」の（しばしば暗黙の）定義によって、同じ現象を見ていても言い方が変わることがあるので注意
• 解析結果を解釈する際には、どのような変数で調整／層別された上での解析結果なのか常に明確に意識しよう
• 解析の「そもそもの目的」を見失わないようにしよう

こんなかんじですかね。はい。

あ。一応付け加えておきますが、本記事の主旨は「中間変量を加えたときの偏回帰係数の解釈には気をつけろ」ということであって、決して「どんなときでも中間変量を加えるな」という中間変量dis話ではないので、その旨よろしくお願いいたします。

（例えば、中間変量を付け加えてみることにより新たな介入可能ポイントが見つかったり*8、媒介変数法（フロントドア基準）を用いて因果効果を推定できたりするという利点もあります）

今回のシリーズの全体まとめ

では、今回のシリーズ全体を軽くまとめてみたいと思います：

今回の一連の記事では「因果関係があるのに相関関係が見られない」以下のケースを見てきました：

• （１）検定力が低い（→過去記事1）
• （２）交絡によるバイアスにより打ち消されてる（→過去記事2）
• （３）合流点バイアスにより打ち消されている（→過去記事3）
• （４）中間変量によりマスクされている（本記事）
• （補）（ピアソンの相関係数を尺度とした場合に）関連の仕方が線形でない（→過去記事4）

はい。いろいろ書きました。お付き合いいただき誠にありがとうございました＞読者さまがた

今、自分はジャワでリゾートを満喫しながらこの記事を書いているのですが*9、一連の記事を通して言いたいことは何だったかなぁと思い返してみると、要するに：

調査観察データ*10における因果効果の解析においては、最も重要なのは先ず「検定力・交絡・合流点バイアス・中間変量」の問題がクリアされているかどうかであって、信頼区間とかp値とかは「それらの問題が既にクリアされている」という前提の上で初めて意味をもつ二次的なものにすぎない

ということかもしれません。（今回のシリーズで示されたいろんな例を思い返してみてくださいね）

なんというか、調査観察データにおける「信頼区間とかp値とか」というのは「牛丼の上に載った”紅ショウガ"」みたいなもので、まあ全然重要じゃないとまでは言わないけれども、解析の本体たる「牛丼本体」の具合の議論をしないで”紅ショウガ"の話ばっかりしてるってのもどうなのよ、と思っちゃうような場面を統計結果の議論に関してけっこう見かけるんですよね。

今回の一連のシリーズを通して、統計解析における「信頼区間とかp値とか」以外の部分の重要さが少しでも伝わったらといいなぁ、と思います。

今回のシリーズのリンク：

いちおう、関連する一連の過去エントリーのリンクを貼っておきます：
因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） - Take a Risk: 林岳彦の研究メモ
"相関"の話＆そのついでに"21世紀の相関（MIC）"の話（ややマニア向け） - Take a Risk: 林岳彦の研究メモ
今回は因果関係があるのに相関関係が見られない4つのケースをまとめてみた（前編：検定力が低い） - Take a Risk: 林岳彦の研究メモ
（中編）今回は因果関係があるのに相関関係が見られない4つのケースについてまとめてみた：交絡・合流点の影響 - Take a Risk: 林岳彦の研究メモ

おまけ：因果グラフの理論がもたらす"明晰さ"について

おまけです。

上記の一連のエントリーでは、統計解析における「因果関係」と「相関関係」の顕れ方の関係について色々と書いてきました。この辺りの話は、今までは各研究分野の中で職人芸的に伝承されていた部分が多く、分野によって様々な語られ方をしてきたものだと思われます（なので、この辺りの術語については、しばしば分野間で統一性がありません）。

手前味噌とはなりますが、上記の一連のエントリーの中では、その辺りの話をなるべく統一的な形で整理し、説明できたのではないかと考えております。で、それを可能にしたのは取りも直さず、一連のエントリーの背景にあるJudea Pearlの因果グラフの理論体系のおかげなんです。

Pearlの因果グラフの理論体系は、今まで様々な分野で経験的に試行錯誤的に語られてきたものに対して、背後にある「共通の論理的構造」に着目することにより、非常に統一的かつ明晰な理論的見通しを与えてくれるんですよね。（その様は感動的ですらあると思う）

例えば、シンプソンのパラドックスなんかはPearlの理論を学ぶとかなり明晰な形で理解することができて、Wikipedia英語版のSimpson’s paradoxの項にも書いてあるとおり（強調引用者）：

As to why and how a story, not data, should dictate choices, the answer is that it is the story which encodes the causal relationships among the variables. Once we extract these relationships and represent them in a graph called a causal Bayesian network we can test algorithmically whether a given partition, representing confounding variables, gives the correct answer. The test, called "back-door," requires that we check whether the nodes corresponding to the confounding variables intercept certain paths in the graph. This reduces Simpson's Paradox to an exercise in graph theory.

Pearlの理論からもたらされる「バックドア基準」を理解することにより、「シンプソンのパラドックス」なんかはグラフ理論における”an exercise”として何の不思議もなく対処することができるようになります*11。

まあ、Pearlの理論を学んだところで私たちが直面している統計学的諸問題が必ずしも解決するとは限らないのですが（逆に「解決不能な問題」であることが判明したりする）、私たちが直面している統計学的諸問題が『何であって何でないのか』が非常に明晰になるというご利益があるかと思います。

そんなこんなで個人的には、「統計的諸問題におけるPearlのグラフ理論の効用」は「哲学的諸問題における分析哲学の効用」にとても似ているかなあ、と思ってます*12。

はい。

現場からは以上です。

＃あ、次回くらいから「はてなブログ」に引っ越そうかな、と思ってます（なんとなく）。

参考文献

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

統計的因果推論の鉄板本その１（Rubin系）

統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)

統計的因果推論の鉄板本その２（Pearl系）

Excelで学ぶ共分散構造分析とグラフィカルモデリング

上の宮川本が難しい方は、こちらからどうぞ。

多変量解析の展開―隠れた構造と因果を推理する (統計科学のフロンティア 5)

統計的因果推論の入門としてはこの本もオススメできます。

統計的因果推論 -モデル・推論・推測-

Pearl師匠の本。いろんな意味で激ムズ。

Modern Epidemiology

この本（3rd ed）の12章の「Causal Diagrams」の章が（それなりに難しいが）英語が読めて、疫学がベースの人には一番有益な因果グラフ入門かもしれない。

（単なる宣伝）

私が詩を寄稿させていただいた自主制作文芸誌が こちらやこちらで買えます（まだまだ赤字らしいのでどぞ4649です）。

*1:結局みなさまにRStudioのデスクトップ版をインストールしていただきました。後で見なおしたらAmazon EC2のマイクロインスタンスではそら落ちるわ、って話ですわ。想定が甘すぎ＞俺

*2:統数研あたりが全国の迷える統計講師役のために登録制で使えるぶっといRStudioサーバー版を立ててくれたりしないかなあ（妄想）。。統計実習って実習時間のかなりの部分が皆の解析環境のセットアップに割かれがちなんすよ。。

*3:10000人のデータ全部を描くとゴチャゴチャしすぎるので最初の500データだけプロットしてます

*4:あるいは、めんどくさいオッサンに絡まれがちな

*5:もちろん、日常的な用法はTPOなどに応じてかなり幅があるものなので、これとは違う言い方も同様に成り立ちうるかと思います。

*6:パス分析における「塩分摂取量→超過死亡リスク」の「総合効果」を見ていることになります

*7:パス分析における「塩分摂取量→超過死亡リスク」の「単独効果」を見ていることになります

*8:例えば、今回の解析例は、「血圧」自体も超過死亡リスクを下げるための効果的な介入ポイントであることを示しています

*9:嘘です。本当は午前5:45のマクドナルドで書いてます

*10:実験計画法に基づいた方法により得られたものではないデータ

*11:あるいは、対処不可能な場合には、その対処不可能性が明晰に理解できるようになる

*12:その意味で私は、一ノ瀬正樹氏のベイジアンネットワーク関連の理論に対する評価にはかなり不満があります。Pearlの理論のキモをちゃんと理解してるのかなあ？

（中編）今回は因果関係があるのに相関関係が見られない4つのケースについてまとめてみた：交絡・合流点の影響

どもです。林岳彦＆オメガトライブです。きみは1005%（消費税込）

さて。

今回は、前回の記事：

今回は因果関係があるのに相関関係が見られない4つのケースをまとめてみた（前編：検定力が低い） - Take a Risk: 林岳彦の研究メモ

のつづきの”中編”になります。本記事では「因果関係があるのに相関関係が見られないケース」の中でも、「交絡・合流点」が関わるケースについて書いていきます*1。

扱う内容の範囲としては、最初の記事：

因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） - Take a Risk: 林岳彦の研究メモ

と重複する部分がかなりありますが、今回の記事では、「仮想例のデータ生成」の段階からRでの計算を交えて説明していきたいと思います。（今回はちょっと「R実習」のような趣になるので、Rの読み書きができないと分かりにくい部分が多々あるかもしれませんが、申し訳ありません）

（最初に断っておきますが、今回もけっこう長いです。本当にすみません）

では書いていきたいと思います。

因果関係があるのに相関が見られないケース（２）：交絡によるバイアスにより打ち消されてる

まずは、「交絡によるバイアスにより、因果効果が見かけの上で打ち消されている」ケースについて見ていきます。

仮想例として、以下のような因果構造のもとでの「河川において農薬が生物種数に与えている影響」を考えていきます：

ここで、「農薬濃度」は調査地点（5000地点*2）における河川中の農薬濃度を、「生物種数」は各地点での河川中の生物（動物プランクトン）の種数を、「近傍の水田面積」は各調査地点の近傍の水田面積を表すとします。

上図の因果構造は、各地点において「河川中の農薬濃度の増加は河川中の生物種数を減少させる」一方で、「近傍の水田面積」の増加は「農薬濃度」と「生物種数」をともに増加させる*3効果を持つことを意味しています。

はい。

では、このような因果構造のケースに対して回帰分析を行うとどうなるか、Rで解析を行っていきましょう（以下のRのコードファイルはConfounding_bias.R ）。

今回は、そもそもの仮想データを生成するところから始めたいと思います。まずは、近傍の「水田面積」のデータ（サンプルサイズn=5000地点）を以下のように生成します：

# サンプルサイズを5000に設定
n.sample <- 5000
# 平均0、標準偏差1の正規分布から5000個の「水田面積」の値をランダムに生成
SuidenMenseki_scaled.v <- rnorm(n.sample,mean=0,sd=1)

ここでは、正規分布からの乱数として5000個のデータが生成され、"SuidenMenseki_scaled.v"に格納されています。解析の単純化のために、データの数値は「平均0, 標準偏差1」となるように既に標準化されている状況を想定しています。（今後、本記事ではデータの数値は基本的に標準化されているものとして考えていきますので、面積や濃度の値が「マイナス」になっていたりしますが、そこは本記事の本質的な議論とは関係してこないので余り気にしないでください）

次は、この「水田面積」のデータから、対応する「農薬濃度」のデータを生成します：

# 平均0、標準偏差1の正規分布から誤差項を作成
error1.v <- rnorm(n.sample,mean=0,sd=1)
# 「農薬濃度＝1.0×水田面積＋誤差」の関係式から「農薬濃度」の値（5000個）を作成
NouyakuNoudo.v <- 1.0*SuidenMenseki_scaled.v + error1.v
# データの標準化
NouyakuNoudo_scaled.v <- scale(NouyakuNoudo.v)

ここでは「農薬濃度＝1.0×水田面積＋誤差」という関係式をもとに「農薬濃度」のデータが作りだされています*4。

さらに、この「農薬濃度」と「水田面積」から、対応する「生物種数」のデータを生成します：

# 平均0、標準偏差1の正規分布から誤差項を作成
error2.v <- rnorm(n.sample,mean=0,sd=1)
# 「生物種数＝0.5×水田面積 - 0.3×農薬濃度＋誤差」の関係式から「生物種数」の値（5000個）を作成
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人）を生成します：

# サンプルサイズを5000に設定
n.sample <- 5000
# それぞれの能力値を平均ゼロ、標準偏差1の正規分布に従いランダムに生成
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)
# デュアスロンタイムデータを「デュアスロンタイム＝ -0.5×自転車力 - 1.5×走力 ＋ 誤差」の関係式により作成
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)
# トライアスロンタイムデータを「トライアスロンタイム＝ -1×自転車力 - 0.5×走力 - 0.5×泳力＋ 誤差」の関係式により作成
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です）。

.

*1:本当は中間変量の話題まで書いていたのですが、長くなりすぎたのでとりあえず「中編」としていったん切り上げました

*2:実際にはこんなに大規模の実測データが手に入ることは殆どありませんが、推定結果が安定するように多めのデータサイズを仮定していきます

*3:近傍の水田環境が河川中の動物プランクトンの供給源となっている状況を想定

*4:濃度がマイナスになるデータがありますが、ここの「農薬濃度」は[http://c4s.blog72.fc2.com/blog-entry-99.html:title=標準化]済みの値であると考えてください

*5:5000地点の全部を描くとゴチャゴチャしすぎるので最初の500データだけプロットしてます

*6:今回の記事における因果関係の定義は[http://d.hatena.ne.jp/takehiko-i-hayashi/20130418/1366232166:title=こちら]

*7:定義により、「農薬濃度を変化"させた"ときの生物種数へ与える影響」

*8:そもそも交絡要因となる変数が測定されてさえいなかったり、あるいは因果構造自体が全く分からないような状況ではお手上げですが

*9:この例を思いついたときは個人的にはとても嬉しかったのです

*10:現実のデュアスロンにおける関係とはかなり違うと思いますがすみません

*11:現実のトライアスロンにおける関係とはかなり違うと思いますがすみません

*12:それぞれのデータの値は平均がゼロになるように標準化済みのものとご解釈ください

*13:スペースの節約のため示した図中では「泳力」は省いています。また、煩雑さを避けるため最初の500データだけプロットしてます

*14:この件において単回帰が「必要にして十分」であることを判断するためのformalな基準として、[http://d.hatena.ne.jp/takehiko-i-hayashi/20120625/1340611310:title=バックドア基準]というものがあります

*15:現実にはそんなかんじの状況も多々あります

*16:なので、ここに「走力」をさらに加えることにより交絡が除去できます。興味のある方はRで実際に試してみてくだされ

*17:あるいは、合流点において"調整"されることにより。データ選別の段階でバイアスがかかるような場合もその選別基準に関わる「変数」上における一種の「調整」とみなすことが可能である

*18:「自転車力」を変化"させた"ときのデュアスロンタイムの変化

*19:前者は本質的に「確率」のレイヤーの話であり、後者は本質的に「因果」のレイヤーの話になります

今回は因果関係があるのに相関関係が見られない4つのケースをまとめてみた（前編：検定力が低い）

どもお久しぶりです。林岳彦です。ローソンなどで売ってるいなばのタイカレーはそうめんのつけ汁として使ってもマジうまいのでオススメです。

さて。

今回は前々回の記事：

因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） - Take a Risk：林岳彦の研究メモ

の続編として、逆のケースとなる「因果関係があるのに相関関係が見られない」ケースについて見ていきたいと思います。あんまり長いと読むのも書くのも大変なので、今回はまずは前編として「検定力の問題」に絞って書いていきます。

（＊今回は上記の前々回の記事での記述を下敷きに書いていきますので、分からないところがあったら適宜前々回の記事をご参照ください）

まずは（今回の記事における）用語の定義：「相関」と「因果」

今回も少しややこしい話になると思うので、まずは用語の定義をしておきたいと思います。（＊細かいところはあまり気にしない方はこの節は飛ばしてもらってもたぶん大丈夫です）

前々回の記事と同じく、今回の記事内では「因果」の定義については：

「要因Aを変化させた（介入）とき、要因Zも変化する」ときに、「要因A→要因Zの因果関係がある」と呼ぶ

ことにしたいと思います（詳しくは、前々回のエントリーをご参照ください）。今回の記事の後編では、この「因果」定義自体の難点にも触れることにもなっていきます。

また、今回の記事内では、「相関」については：

得られたデータにおける項目Aと項目Zのあいだに何らかの関連がみられる（＝独立でない）ときに、両者間に「相関がある」と呼ぶ

ことにします。より一般には、「相関」というと、 ピアソンの相関係数で示されるような項目間での「直線的な関連性」を指すことが多いですが、今回の記事で言及するのは基本的に「直線的な関連」に限定した話ではないので、上記のようなかなりユルい定義を使っていきます（より詳しい議論は前回の記事をご参照ください）。

そして「相関が見られる／見られない」という表現については、「データにおける項目間にある関連性が、偶然に生じたものと区別できるほど大きい*1／区別できないほど小さい」という意味で用いていきます。（そもそも相関は本質的に連続量的なものなので、「見られる／見られない」という二分法自体がスジが悪いともいえますが、ここでは日常的用語法に従い「見られる／見られない」という表現を使っていきます*2）

また、今回の記事では、あくまでも狭義の相関（ ピアソンの相関係数r）ではなく、データにおける項目間に「何らかの関連性」が見られるかどうかを問題としていくので、例えばt検定で「2群間の平均値に差がある（＝2群間での「処理の違い」と「結果」の間に関連性がある）」ことをもって「相関がある」という言い方もしていきます（まあ基本的に、2群間で平均値に差が見られれば相関も見られるので同じことになります）。

因果関係があるのに相関が見られないケース：（１）検定力が低い

はい。では「因果関係があるのに相関が見られない」ケースを見て行きましょう。

まずは、なんといっても、そもそも「検定力が低い」ケースです。

検定力とは、「（偶然によるものではない）差や関連性がある場合に、それを統計的に有意な差や関連性として検出できる確率」です。そして、結論から言うと、その「検定力」の大きさは、「サンプル数サンプルサイズ*3」「影響（関連性）の大きさ」「データにおけるバラつきの大きさ」「有意水準」に本質的に依存します。

では、具体的な理解のために、仮想例を見ていきましょう。

今回は、「新開発のキンギョのエサが、キンギョの体重に与える影響」のケースを考えていきます。

あなたはキンギョのエサの開発者で、「新しく開発したエサが（従来のエサに対して）どのくらいキンギョの体重を増加させるのか」を知りたいとします。ここで、生まれたばかりのキンギョの稚魚をA群/B群へとランダムに分け、A群には「従来のエサ」を、B群には「新しく開発したエサ」を与えつづけて1年間育ててみるという実験をしてみます：

ここで、「従来のエサで育てたときの平均体重は30g、標準偏差は5gで、新開発のエサはキンギョの1年後の体重を平均+3g増加させる（標準偏差は変化させない*4）」という状況を考えていきましょう。

果たして、各群に何匹のキンギョを用いればこの2群の「差」を検出できるでしょうか？

では、Rで色々と計算していきたいと思います。まずは、試しに散布図の例から描いてみます*5。

この図は各群のキンギョが10匹（各群がn=10）のケースについて、横軸の左側は従来のエサ（A群）、右側は新開発のエサ（B群）、縦軸はそれぞれの群における「1年後の体重（g）」を表しています（RコードはKingyo_Ex.R ）。このプロットを見ると、新開発エサのB群の体重の方が大きいように見えますが、その差の程度はビミョウなかんじですね。。

では、実際にサンプルサイズがA群、B群ともにn=10匹の場合のp値を計算してみましょう。今回は「正規分布する2群のサンプルにおける平均値の差を見る」のが目的なのでt検定になります*6：

Rで解析を行うと次のようになります：

> n <- 10; controlmean <- 30;
> truedifference <- 3; truesd <- 5;
>
> weightA.data <- rnorm(n,mean=controlmean,sd=truesd)
> weightB.data <- rnorm(n,mean=controlmean+truedifference,sd=truesd)
>
> ttest.res <- t.test(weightA.data, weightB.data, var.equal=T, alternative="less")
> ttest.res

     Two Sample t-test

data:  weightA.data and weightB.data
t = -0.6752, df = 18, p-value = 0.2541
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 2.454056
sample estimates:
mean of x mean of y
 31.15390  32.71873

今回の各群n=10匹の例では、p値（p-value）はp=0.25になりました*7。有意水準をp=0.05とすると、この場合では"統計的に有意な差"を検出できなかったことになります。

さて。では、このような各群n=10匹の場合のシミュレーションを100回繰り返してみましょう（A群、B群のデータは各回のシミュレーションごとに新たなものを生成しています：RコードはKingyo_ttest100sim.R ）：

各群n=10のシミュレーションを100回繰り返すと、以上のような結果になりました。横軸は100回のシミュレーションを、縦軸はそれぞれのシミュレーションで得られたp値を示しています。赤色の破線はp=0.05の水準を示しています。この例では、有意水準をp=0.05とすると100回中の37回において（のみ）統計的に有意な差（p<0.05）を検出できているという結果になりました。

この「（偶然によるものではない）差がある場合に、それを統計的に有意な差として検出できる確率」が、「検定力」と呼ばれるものになります。
この「検定力」がなぜ重要かというと、要するに：

「検定力」が低い場合に、「因果関係があるのに相関関係が見られない」というケースが数多く生じる

からです。

例えば、上記の図の例では、本当は「新開発のエサによって体重が+3g変化している（＝因果関係がある）」にもかかわらず、「100回中37回」しかそこに「統計的有意差が見られない（＝相関関係が見られない）」という状況になっています。

もし、こんな検定力の低い（10回に4回弱しか差を検出できない）デザインの実験結果に基づいて「差がない」という結論になり、あなたが心血を注いで新開発したエサがお蔵入りになってしまったら悲しいですよね。。。

検定力が低い（あるいは、検定力に関する意識が低い）というのは、とても由々しき問題なのです。

どのような場合に検定力が低くなるのか

では、どのような場合に検定力が低くなるのかを見ていきましょう。

まずは、上記と同じキンギョの実験の例を念頭に、「コントロール群の平均体重30g、標準偏差5g、新エサ投与時の増加体重+3g、各群n=10、有意水準p=0.05（片側検定）」の条件における検定力をRを用いて計算してみます。t検定を対象とした検定力の計算には、power.t.test()という関数を用います：

> power.t.test(n=10,delta=3,sd=5,sig.level=0.05,type="two.sample",alternative=("one.sided"))

     Two-sample t test power calculation

              n = 10
          delta = 3
             sd = 5
      sig.level = 0.05
          power = 0.3617837
    alternative = one.sided

検定力（power）は0.36と計算されました。同じパラメータ値を用いた上述のシミュレーションで得られている結果（100回中37回で有意差あり）と良く整合してますね。この検定力の「0.36」という値は低いのか高いのかというと、「よく分からないからもうコイントスで決めようぜ（←やけくそ）」というのでも検定力は0.5あるわけなので、まあ少なくともとても高いとは言えない値です*8。

次は、検定力に対するサンプルサイズの影響をみていきましょう。上記と同じ「コントロール群の平均体重30g、標準偏差5g、新エサ投与時の増加体重+3g、有意水準p=0.05（片側検定）」という状況でサンプルサイズn（1群あたりのキンギョの匹数）だけを変化させた検出力の算出を行います：

はい。サンプルサイズが増すに従って検定力が上昇するのが見て取れます（Rコードはkensyutsuryoku.R ）。もし8割の確率で「群間で体重に差がある場合に*9、それを統計的に有意なものとして検出できる」こと（検定力=0.8）を目指すのなら、サンプルサイズとして少なくとも各群40匹程度は必要となることが分かります。

一方、サンプルサイズが各群20匹程度の場合には、検定力が0.5くらいであり、おおよそ半々の確率で「新エサと体重の間に、因果関係があるのに（統計的に有意な）相関がみられない」ことが予測されます。

では、次は体重における「影響の大きさ」「ばらつきの大きさ」の影響を見てみます。今までと同じくデフォルトは「コントロール群の平均体重30g、標準偏差5g、新エサ投与時の増加体重+3g、有意水準p=0.05（片側検定）」で、体重の変化量と標準偏差のみを変えてみましょう（Rコードはkensyutsuryoku.R ）：

他の要因が同じ場合、「変化量（影響の大きさ）」が大きくなるにつれて検定力は増加します（左図）。上図の例では、新エサによる体重の変化量が+6gを超えると検定力も0.8以上になっていますが、体重の変化が+3gくらいだと検定力は0.4以下に留まります。つまり、「変化量」が小さいほど検定力が低下するため、例えば、変化量が小さい場合にも同じ検定力を維持するためには「サンプルサイズを増す」などの対処が必要になってくる、というわけです。

また、他の要因が同じ場合、「体重のバラツキ」が大きくなると検定力は減少します（右図）。上図の例では、体重の標準偏差が2gくらいであれば検定力は0.8を超えますが、標準偏差が5gくらいになると検定力は0.4を下回ってきます。これは例えば、「同じ実験デザインの実験を行っても、実験の下手な人がやると（結果におけるバラツキが増えるため）有意差が出にくくなる」という現象に相当するものです。逆の言い方をすると、検定力を上げるための方法として「バラツキ」を抑える工夫をする、という方向性もあることになります。

で、実は、上の「変化量の大きさ」と「バラツキの大きさ」は相対的なもので、検定力は「変化量／標準偏差」に依存する、とよりシンプルに両者の影響をまとめることができます。例えば、左図ではデフォルトの標準偏差として5gが用いられていますが、「変化量／標準偏差=1」（＝変化量5g）のとき、検定力は0.6強になっています。一方、右の図ではデフォルトの変化量として+3gが用いられていますが、こちらも「変化量／標準偏差=1」（＝標準偏差3g）のとき、検定力は同じ値（0.6強）になっているのが分かるかと思います。

この「変化量／標準偏差」の値（つまり、変化量のスケールを標準偏差で標準化した値）は一般に「効果量（effect size）」と呼ばれるものに対応します*10。

最後に、有意水準を変化させてみましょう。今までと同じく「コントロール群の平均体重30g、標準偏差5g、新エサ投与時の増加体重+3g、各群のサンプルサイズn=10匹」で、有意水準だけを変化させてみます（Rコードはkensyutsuryoku.R ）：

（あたりまえなのですが）有意水準が大きくなるにつれ、検定力が上がっていきます。有意水準が大きく（緩く）なると「本当は差がないのに、差があると判定」する確率（第一種の誤り）が大きくなる一方で、検定力は大きくなり、「本当は差があるのに、差がないと判定」（第二種の誤り）する確率が小さくなります。このように、有意水準と検定力は本来的にトレードオフの関係にあるわけです。

ここで、「有意水準」と「検定力」のどちらを重視する必要があるのかは本来はケースバイケースなのですが、多くの場合、「第一種の誤り」の方を制御することを重視して*11、有意水準の方を先ず固定（p=0.05やp=0.01とか）にして考えていく方法が一般的となっています。

ここまでの小まとめ

この辺で、いったんまとめてみます：

• 検定力が低いと、「因果関係があるのに相関関係が見られない」ケースが多く生じる
• 検定力が低くなるのは以下のケースである：
  • (1) サンプルサイズが小さい
  • (2) 効果量（＝変化量／標準偏差*12）が小さい
  • (3) 有意水準が小さい

というかんじになるかと思います。（ここまではよろしいでしょうか？）

でも「やみくもにサンプルサイズを増せばいい」というわけでもないの：本質は効果量にあり

さて。

ここまで、「検定力」は「サンプルサイズ」「効果量」「有意水準」に依存することを見てきました。仮説検定の枠組みにおいて、これらの4つの量は常にお互いに依存し合う"カルテット"的な量になっています。そして、これらの中でも、多くの場合に実験や調査の前にわれわれが主体的に大きく変えられる余地があるものは「サンプルサイズ」です。そのため、検定力と聞くと「とにかくサンプルの数を増せばいいんでしょ？」と思われる方もいるかもしれません。

まあそれはそうなのです。が、ただ、そうは言っても「やみくもにサンプルの数を増せばいい」というわけでもありません。

そもそも、検定における「検定力」「サンプルサイズ」「効果量」「有意水準」の中で、私たちが知りたい「調査対象において実際に起きていること」を最も直接的に反映しているのは何といっても「効果量」なのです。その一方、「有意水準」「検定力」というのは「実験・調査デザイン（のクオリティ）」の方をより直接的に反映する量といえます。

その意味で、特に論文の査読という「実験・調査デザインのクオリティ」を重視するプロセスにおいて「有意水準」や「検定力」が重視されてきたのは理解できますが、そのせいで「効果量」の方が軽視されるようになってしまったら、それはもう本末転倒と言えるでしょう（例えばこちら）。

（殆どの場合において）大事なのは「調査対象において実際に起きていること」を理解することであり、そのためには「効果量」を中心に考えを進めていくことが重要です。

サンプルサイズを決定する際にも、そもそもの本来的な目的に照らして「どの程度の"差"が実質的な差と言えるのか？」という点をまず明確にし、その「実質的な差（実質的な効果量）を検出するためにはどの程度のサンプルサイズが必要なのか？」という順番で考えていくことが王道といえます。（ここでちゃんと考えることをサボると、あとで「nが足りな〜い...」という番町皿屋敷的な状態になるので気をつけましょう）

逆に言うと、サンプルサイズをやみくもに多くして、「瑣末な差」に統計的有意性を見いだしたとしても、それは「現実を知る／現実に対処する」ためには何の役にも立ちません*13。

（ついつい忘れがちですが）統計解析の際にはいつでもその本来的な目的を見失わないように注意しましょう。

今回のまとめ：

では、今回の内容をまとめます。（殆どさっきのまとめと重複しますが）

• 検定力が低いと、「因果関係があるのに相関関係が見られない」ケースが多く生じる
• 検定力が低くなるのは以下のケースである：
  • (1) サンプルサイズが小さい
  • (2) 効果量（＝変化量／標準偏差）が小さい
  • (3) 有意水準が小さい
• 検出力はやみくもに高ければ良いというものでもない（本質は「効果量（effect size）」にあり）
• わたしの検定力は530000です

はい。

ちなみに今回の記事では仮説検定の枠組みを基にして書きましたが、またそもそも論を言うと、「効果量」に着目するのであればそもそも仮説検定の枠組みよりも「効果量の信頼区間」や「効果量の推定分布」を求めたほうが断然スジが良いと言えます。この辺りを書きはじめるととても長くなるので、ぜひ以下の本をご一読いただければ幸いです（すごく良い本だと思うので）：

伝えるための心理統計: 効果量・信頼区間・検定力

（もしくは、みんな仮説検定のことなんか忘れてベイジアンになっちゃえばよいのにと思います）

はい。

では、次回の後編では、「因果関係があるのに相関関係が見られないケース」の中でも、交絡要因や中間変量が関連するものについて書いていきたいと思います。

（今後はコンスタントに更新していきます）

今回の参考文献：

検定力分析入門

検定力分析の入門書。さすがの分かりやすさ。事例（心理学中心）も多くオススメです。

伝えるための心理統計: 効果量・信頼区間・検定力

上記の豊田本より高レベルですが、読み応えはあります。仮説検定が廃れてきた理由と経緯についても詳しく書かれています（あまりその辺りのところを詳しく書いてくれている本はあまりないので大変ありがたい）。こちらも大変オススメです。

（単なる宣伝）

私が寄稿させていただいた自主制作文芸誌が こちらやこちらで買えます（まだまだ赤字らしいのでどぞ4649です）。

.

*1:counterfactualさんのブコメを見てここに量的な表現がくるのはおかしいかなと思ったので改変

*2:ちなみに、私は因果関係に対しては「ある／なし」、相関に対しては「見られる／見られない」「現れる／現れない」の語を使うことが多いですが、これは因果は本質的に存在論的なものであるのに対し、相関は本質的に現象論的なものである、という理解の仕方が背景にあります

*3:ブコメで「サンプル数」の語は誤りとの指摘あり。まじどもサンクスです。勉強になります。以下の文では瑣末になるのでこっそり修正しました

*4:つまり等分散性を仮定

*5:各群のデータは正規分布すると仮定します

*6:体重の「増加」に興味があるので片側検定を用います

*7:その都度サンプルをランダムに発生させているので、コードを実行するごとに値は変わります

*8:一般には検定力は0.8くらいあることが望ましいと言われていたりします→詳しくは大久保・岡田(2012)

*9:厳密に言えば片側検定なので「B群の体重の方が大きい場合に」

*10:ここはちょっと大雑把な言い方になっています。実際には効果量の定義にはいろいろあります。詳しくは大久保・岡田(2012)などをご参照ください

*11:あるいは、何も考えずに単に慣習に従って

*12:「効果量」の定義はいろんな形があり得ますが、ここではとりあえずこの形で

*13:もちろんそこで得られた「効果量」には意味があります／まあそれがサンプルサイズの増加に費やした努力量に見合ったものかは分かりませんが

"相関"の話＆そのついでに"21世紀の相関（MIC）"の話（ややマニア向け）

どもです。林岳彦です。息子の3DSにバーチャルコンソールの「ソロモンの鍵」を密かに入れました（まだ３面）。

さて。

前回の記事：

因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） - Take a Risk：林岳彦の研究メモ

につきましては沢山ブクマ等をいただき大変ありがとうございました*1。大変感謝しております。

さて。上記記事について、ublftboさんから「相関関係の定義が書かれていないのでは」（相関と因果 - Interdisciplinary）とのご指摘をいただいたきました。

ご指摘は確かにごもっともですので、今回は「相関」概念についてと、そのついでに近年に開発された"21世紀の相関（MIC）"の話について私なりに書いてみたいと思います。

（以下、ややマニア向けの話になるかもしれません。あと前回ほどではないですが、それなりに長いです。）

広義の「相関」（より適切には「関連／association」と呼ばれるもの）

まずは、前回の記事の補足的なところから話をはじめたいと思います。

前回の記事において「相関」という語を私がどういう意味で用いていたかというと、あまり深くは考えてはいませんでしたが、改めて文字にすると「得られたデータにおける項目Aと項目Bの間に何らかの関連が見られる＝相関がある」という意味で用いていたと思います。

これはかなりざっくりした用法で、かなり広義の意味での「相関」という語の使い方と言えます（日常言語における「相関」の語のニュアンスの方を強く反映した用法とも言えるかもしれません）。

Wikipediaの"correlation"の項を見てみると：

In loose usage, correlation can refer to any departure of two or more random variables from independence, but technically it refers to any of several more specialized types of relationship between mean values.

という記述がありましたが、この前半の"loose usage"における"any departure of two or more random variables from independence"というところが、前回の記事において私が念頭においていた「相関」の語の意図するところになります。

ここで、"departure of two or more random variables from independence"というところの意味が良く分からない方も多いかもしれないので、「関連／association」と「独立／independence」の関係についてちょっと説明を加えてみます：

さて。そもそも「得られたデータにおける項目Aと項目Bの間に何らかの"関連"が見られる／見られない」というのは、統計学的にどう表現されうるでしょうか？

一般的には、Aが起こる確率P(A)と、Bが起こる確率P(B)、および、AとBが同時に起こる確率P(A,B)の関係が：

P(A,B) = P(A)P(B)

のときに、AとBは「独立／independence」と呼ばれます。また、AとBが独立であるとき、それらの間には「（統計的に）関連がない」とされます。

例えば、２つのコインA, Bを投げたときに、それらがオモテとなる確率がそれぞれ「 P(A=オモテ)=0.5、P(B=オモテ) = 0.5 」であるとしましょう。

ここで「AがオモテでBもオモテ」となる確率 P(A=オモテ, B=オモテ）が、「 P(A=オモテ, B=オモテ）= P(A=オモテ) x P(B=オモテ) = 0.5 x 0.5 = 0.25 」であるときには、「コインAを投げてオモテになる確率」と「コインBを投げてオモテになる確率」は「独立」であり「関連がない」ということになります。

逆に、 P(A=オモテ, B=オモテ) が0.25 から逸脱する場合には、コインAとコインBのあいだに「何らかの関連性」が推測されることになります。

ちなみに上の式を別の形で書くと：

P(A) = P(A|B)

とも書くことができます。これは、「Bがどうであるか」はAの確率に影響を及ぼさない、ということを示しています。また別の言い方をすると、「Bに関する情報は、Aに関する予測の役に立たない（＝何の情報ももたらさない）」という言い方もできます。

つまり「AとBが独立である」というのは、AとBの関係が「関連がない」「影響を及ぼさない」「予測の役に立たない」「情報をもたらさない」などなどの意味に対応することになります。

はい。

というわけで、「AとBの間に広義の意味で相関がある」は数学的／統計学的にいうと「AとBは独立でない」という状況に対応するものになります。

一般には、この意味での「広義の相関」については、「相関correlation」よりも「関連association」という語が使われるのが普通ですので、細かい議論をする際には区別して使うことが望ましいでしょう。（その意味で、前回の私の記事の用語法はあまり良くないと言えます*2）

では次は、ご存知の方も多いと思いますが、狭義（というかより一般的な用法における）「相関」について説明してみたいと思います。

いわゆる「相関」は直線的関係の指標である

一般に、統計学の文脈において「相関」と言った場合には、「ピアソンの相関係数」に基づくものを意味することが殆んどかと思われます。

Wikipediaの「相関係数」の項からピアソンの相関係数の数学的定義を引用すると：

となります。

さて。では、ピアソンの相関係数の直感的な性質を見てみましょう。

（ピアソンの）相関係数の特徴は、データ間の直線的関係のみを見ていることにあります。

直線的な関係を見ているということは、AとBの間に「明らかな関連／association」がありそうな場合にも、相関係数（correlation coefficient）は低い値になりうる、ということを意味しています。

「百聞は一見にしかず」なので、Wikipediaの図を見てみましょう：

http://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A1%E3%82%A4%E3%83%AB:Correlation_examples2.svg より引用

ここでそれぞれの図はデータの散布図を表していて、上の数字はそれぞれの相関係数を表しています。

上の図から：

• 直線的関係でばらつきが全くない場合は相関係数は1（または-1）になる
• 直線的関係*3の傾きの大きさは相関係数の大きさには関係ない
• ただし傾きが完全にフラットのときは相関係数はゼロ*4
• 直線的関係においてばらつきが大きいと相関係数はその分小さくなる
• 非直線的な関連を評価したい場合には相関係数はあまり役にたたない

という「相関係数」の特徴がわかると思います。

はい。いわゆる「相関係数」とは、こういうものなんです。

さて。上記の特徴を理解して使えば相関係数は大変便利なものです。しかし、世の中の全てが「直線的関係」だと思うなよ、と言われたらそれはまさしくそうであります。

そこでMICですよ（ドヤ）。

21世紀の"相関"：MICとは？

MIC（Maximum Information Coefficient : MIC@Wikipedia）とは、すごくざっくり言うと「どんな形でも対応可能な"相関"係数」です。

2011年のScienceで出版されたReshef et al. 2011において発表されたもので、そのときのScience誌の解説文では「A correlation for the 21st century」なんて書かれたりしています（スゴイネ!）。

MICはピアソンの相関係数のようには単純な数式では表せず、コンピュータによってゴリゴリと計算されます。基本的なアルゴリズムとしては、データ散布図を様々な数のグリッド（＝解像度）で区切っていきながら、様々な解像度の値において相互情報量が最大（＝各グリッド内に含まれるデータ密度のコントラストが最大となるようなイメージ）となるような区切り方を決定し、それらを規格化したのちの最大の情報量をMICの値として選択しているようです*5。

ちなみに相互情報量を数式で表すと：

となっており、独立（ P(x,y)=P(x)p(y) ）のときには相互情報量はゼロとなることがわかります。

MICの特徴は、どんな形のassociationでも定量化できるできるところにあります。例えばこんな感じです：

http://lectures.molgen.mpg.de/algsysbio12/MINEPresentation.pdf より引用・改変

これは、一番左の列のタイプのデータの場合に、MICとピアソンの相関係数がそれぞれどのような値をとるかを示しているものです。ピアソンの相関係数では検出できていないような非線形の場合においても、MICでは高い値を示すことがわかります。（ちなみにMICがとる値の範囲は0から1までになっています）

MICの特徴として、データがばらつくに従ってその値が低下することも挙げられます：

http://lectures.molgen.mpg.de/algsysbio12/MINEPresentation.pdf より引用・改変

この辺りの性質は、ピアソンの相関係数の性質を良く受け継いでおり、直感的にも違和感のないものです。

MICはかなりgeneralなものなので、基本的にはどんな対象にでも適応できます。その中でも有望な応用例として、遺伝子発現における「非線形的関連の検出」が挙げられているようです。

幸い、最近RでMICを簡単に計算するための"minerva"というパッケージが出たようなので、それを使って計算も試してみたいと思います：

install.packages("minerva")
library(minerva)
data(Spellman)
Spellman <- as.matrix(Spellman)
res <- mine(Spellman,master=1,n.cores=1)

ここで"Spellman"というデータセットには、CDC15 Yeast Geneの4382個の転写産物の量を時系列（23 time points）で計測したデータが入っています（詳しくはこちら）。mine関数においてMICの値が計算されており、"res"にはその結果が格納されています。

"res"の中身を見て、MICが高い値になっていた2つの転写産物の例をピックアップしてみると以下のようなパターンになっていました：

いずれも、ピアソンの相関係数では低い値となるケースですが、MICでは非線形的な関連が捉えられているようです。

正直ちょっと「ほんまかいな」と思わないでもないですが、研究の「とっかかり」を得る分には十分なのかなあとも思います。

まとめ

まとめます。

今回は：

• 「（統計的に）関連がない」とは「（統計的に）独立である」ということ
• （統計学の文脈で）「相関」といえば一般には「ピアソンの相関（係数）」のことを指す
• （ピアソンの）相関係数は直線的関係しか評価していない
• 直線的関係に限らない「関連度」の一般的な指標としてMICなんてのがあります

てな話でした。

次回は、「因果関係があるにもかかわらず相関関係（＊より正確には"統計的関連"）が生じないケース」についてまとめたいと思います。

関連情報などまとめ

• Correlation and dependence @Wikipedia（link）
• 確率論的独立性 @Wikipedia （link）
• MIC @Wikipedia (link）
• MICの元論文（Reshef et al. 2011 in Science） （link）
• 同号に載っていたMICの解説（link）
• Reshef et al. 2011の補遺（link)
• MICの解説資料（オススメ）link
• RのMICが使えるパッケージ"minerva" （link）

＃関連しない単なる宣伝：超文学フリマにおいて私も寄稿させていただいた文芸誌「線と情事」買えます！

私も寄稿させていただいた「線と情事」が、ニコニコ超会議2と併せて行われる「超文学フリマ」において以下の要領において販売される予定です。

「超文学フリマ」in ニコニコ超会議2:
2013年04月28日（日）10:00〜17:00幕張メッセ（「ニコニコ超会議2」内）ブース: 【イ-03】

ニコニコ学会βなどのついでにでもお立ち寄りいただければ幸いです。何卒よろしくお願いいたします。

＊通販でも買えます
甘茶茂、しまおまほ、島田虎之介、関根美有ほか「線と情事 創刊号」 - タコシェオンラインショップ
創作・オリジナル 同人誌専門店 サッシ / 「線と情事」 / NECOfan
.

*1:小心者なのでブクマ数が伸びすぎてるプレッシャーで風邪を引きました

*2:まあ、ざっくり「相関と因果」という場合には広義の意味での「相関」が含意されていると考えたほうが良い気もしますが、それにしてもどこかで一度定義しておいたほうが良かったですね

*3:相関の話がメインなので、ここでは回帰直線を意図していません。回帰と相関の違いについては→ここなど

*4:と思ったけど、ばらつきがない場合に限っては厳密に言うと相関係数の式の分母がゼロになるから計算できないってのが正解かも／←コメント欄でのご指摘に従い修正

*5:ここあんまり自信ないです。詳しくはこの解説資料および原論文の補遺を参照推奨

因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き）

どもっす。林岳彦です。ファミコンソフトの中で一番好きなのは『ソロモンの鍵』です*1。

さて。

今回は、因果関係と相関関係について書いていきたいと思います。「因果関係と相関関係は違う」というのはみなさまご存知かと思われますが、そこをまともに論じていくとけっこう入り組んだ議論となります。

「そもそも因果とは」とか「因果は不可知なのか」のような点について論じるとヒュームから分析哲学（様相論理）へと語る流れ（ここのスライド前半参照）になりますし、統計学的に因果をフォーマルに扱おうとするとRubinの潜在反応モデルやPearlのdo演算子やバックドア基準（ここのスライド後半参照）の説明が必要になってきます。

その辺りのガッツリした説明も徐々に書いていきたいとは考えておりますが（予告）、まあ、その辺りをいちどきに説明しようというのは正直なかなか大変です。

なので今回は、あまり細かくて遭難しそうな話には立ち入らずに、「因果関係がないのに相関関係があらわれる4つのケース」についてあくまでもざっくりとした説明を試みてみようと思います。

(＊最初に断っておきますが、以下、長いです)

まずは（とりあえずの）因果関係の定義から

まずは、（とりあえずの）「因果関係」の定義をしておきたいと思います。

因果の定義には色々なスタイルがありうるのですが、今回の記事内では：

「要因Aを変化させた（介入した）とき、要因Zも変化する」ときに、「要因A→要因Zの因果関係がある」と呼ぶ

ことにしたいと思います。

このような「介入に基づく因果の定義」は非常に日常的／普遍的なものです。たとえば、初めて訪れた大きな教室において「壁にあるスイッチと照明電灯の因果関係」を把握したい場合を想像してみましょう。

このとき、私たちはおそらくスイッチを適当に押し、「どのスイッチを変化させると、どの電灯の状態が変化するか」の対応をみて、「各スイッチと各電灯」についての「因果」を了解します。このような「介入による変化にもとづく因果関係の判断」は私たちが普段の生活の中で日常的に行なっていることです。

因果関係について本質論的に考えていくと遭難しがちなので、とりあえず今回はこのようなカジュアルな形で「因果関係」という概念を捉えておきます。

因果関係なしの相関関係：（１）偶然によるケース

では、因果関係がないのに相関関係があらわれるケースを見て行きましょう。

まず第一に挙げられるのは、単なる偶然によるケースでしょう*2。

独立した2つの現象に関連性が現れることは、単なる偶然によっても生じえます。もっとも単純には、「二つのサイコロを転がしたら同じ目がでた」なんてのがそれにあたります。ここではもちろん「同じ目がでた」と言っても、２つのサイコロの目の間に「因果関係」はありません。つまり、「サイコロAの目を変化させると、サイコロBが同じ目に変化する」という関係があるわけではありません。

もう少し複雑な例として、25人の幼稚園児に「うまい棒」と「チロルチョコ」を、10枚のコインを投げたときのオモテの数の個数だけ与えるような場合を考えてみましょう*3。

まず、25人の園児に10枚のコインを投げたときのオモテの数の個数だけ「うまい棒」を与えます。つまり、1人目が10枚のコインを投げたときのオモテの数が「5」なら1人目には「5個のうまい棒」を与え、2人目が10枚投げたときのオモテの数が「3」なら2人目には「3個のうまい棒」を与え…という調子で、25人に対してそれぞれオモテが出た数だけ「うまい棒」を与えていきます。

その後、最初に与えたうまい棒の数とは全く関係なく、25人それぞれに新たにコインを10枚投げたときのオモテの数の個数だけ「チロルチョコ」を同様に与えていきます。このようにして、25人全員に「うまい棒」と「チロルチョコ」の配布を終えたとします。

この場合、もちろん、各園児が持つ「うまい棒」と「チロルチョコ」の個数の間に因果関係はありません。つまり、「ある園児が持つうまい棒の個数を変化させると、その園児が持つチロルチョコの個数が変化する」という関係はありません。

さて。

では、このとき、各園児が持つ「うまい棒」と「チロルチョコ」の個数に、どのくらいの頻度で「有意な相関」が検出されると思いますか？ Rによるシミュレーションをやってみます（ソースコードは Umaibou-Tiroru.R ）：

この図は、25人の園児のそれぞれに「うまい棒」と「チロルチョコ」を、それぞれ10枚のコインを投げたときのオモテの数の個数だけ与えるシミュレーションを100回くり返したときの結果です。

横軸はそれぞれのシミュレーション（1〜100回目）を、縦軸はそれぞれのシミュレーションにおいて得られた「各園児が持つうまい棒とチロルチョコの個数」の相関係数のp値を示しています。赤線はp=0.05となるところを示しており、赤線を下回ると有意差（有意水準 p=0.05）があるということになります。

上の図からは、100回中5回くらい、有意水準を下回る場合があり、相関係数に「有意差」が現れていることが分かります。つまり、100回に5回くらいは、各園児が持つ「うまい棒」と「チロルチョコ」の個数に偶然による「有意な相関」が現れるわけです*4。

ちなみに、上の図は園児の数が25人でしたが、園児の数を25,000人に変えても同じような結果が得られます：

これは「統計的有意差」の意味を正しく理解している人々には極めてあたりまえのことなので、これらの図を見て「えっ！？」と思ってしまった方々は、この機会に"有意差"というものを復習しておきましょう。（例えばここ）

また、このような件に関する特に注意が必要なケースとしては、同一のサンプルに対して多数の仮説検定を適応（＝多重比較）しているときが挙げられます。このような場合には、上記のような「偶然によって生じる相関関係」を有意なものとして特に拾いやすくなることが知られています。この多重比較の問題については以前に書きましたので、ぜひ以下をご参照ください：

無から有（意差）を生む：多重比較でウソをつく方法 - Take a Risk：林岳彦の研究メモ

調査観察データにおける「相関関係」が偶然によるものかどうかを知りたい場合には*5、可能な場合には再度の調査を行ったり、あるいは独立に行われた類似の研究間で一貫した結果が得られているかを調べると良いでしょう*6。多数の研究で一貫した結果が出ている場合には、その「相関関係」が偶然によるものとは考えられません。

因果関係なしの相関関係：（２）因果の流れが「逆」のケース

次は、因果の流れが「逆」のケースについて見ていきたいと思います。

例として、「事故多発注意の看板」について考えていきましょう。「事故多発注意」の看板がある幾つかの場所と、そのような看板がない幾つかの場所で、事故の発生率を調べてみたとします。そのような調査の結果、おそらく「事故多発注意」の看板のある場所のほうが、ない場所よりも、事故の発生率が高いかもしれません。図や数式で表すと：

P(事故の発生|看板あり) > P(事故の発生|看板なし)

と書けるような状況です。

このような場合に、「看板あり→事故発生率高い」という因果関係はあると考えられるでしょうか？

これは、普通に考えると、因果の流れが「逆」ですよね。このような場合には、おそらく「事故の発生率が高い→看板の設置」というのが順当な因果の順序であり、その結果として、「事故の発生率」と「看板の設置状況」に相関関係が生み出されているものと考えられます。

このような場合には、「事故の発生率」と「看板の設置状況」の相関関係は、「看板あり→事故発生率上昇」の因果関係の存在を意味しません。また、もちろん「看板設置への介入（看板を外す）→事故の発生率の減少」という変化も期待できません。

看板の例はあまりにも安直すぎるかもしれないので、もうちょっとリアルに難しい例も考えてみましょう。

サンゴとその捕食者の例を採り上げてみます（この話はこのtogetterの内容を基にしていますが、本記事ではあくまでも説明のための仮想例としてディテールは無視して取り扱います）。

サンゴの保全のための調査から、サンゴの生存率とサンゴの捕食者Oの個体数に以下の相関関係が示されているとしましょう。また、捕食者Oは実際にサンゴを捕食していることがフィールドでの観察から分かっているとします。

このとき、「捕食者Oの増加→サンゴの生存率の減少」という因果関係を想起するのは自然なことかもしれません。

もしこのような因果関係が存在するならば、「捕食者O」を減少させることにより「サンゴの生存率」を増加させることができそうです。

はてさてしかしながら：

より詳細な調査から、「捕食者Oは死にかけのサンゴしか食べない」ことが分かってきたとします。このとき、捕食者Oは実は生態系の中でスカベンジャー的役割を果たしていたということになります。

こうなると、「サンゴの生存率が低下→スカベンジャーである捕食者Oが増加」という逆の因果が真である可能性も出てきます。もしこの形の因果が真ならば、「捕食者Oを減少させることによりサンゴの生存率を増加させる」という保全施策は全く効果を及ぼさないことになります。（むしろ、スカベンジャーを排除することによりサンゴの健全な新陳代謝が妨げられる可能性さえあるかもしれません）

そして、このどちらの「因果の向き」がより真に近いのかは、基本的には現場での観察 and/or 介入によってしか明らかにすることはできません*7。

はい。

このサンゴと捕食者の例は、「因果の向き」を正しく認識し、適切な統計的因果推論を行うためには、対象とする現象についての適切な背景知識を持っていることが本質的に重要であることを示しています。（因果の森に一歩足を踏み入れたからには、stepwise AICみたいなオートメーションな形で話を進めることはできないのです！）

対象とする現象に対する理解が（その現象の複雑さに比して）乏しい場合には、「因果の向き」を正しく理解することも案外と難しいものです。油断して自分の思い込み／予断で進みすぎないように、注意しましょう。

あと、以下のような、複数の要因を経た「逆向きの因果」も非常に気づきにくくやっかいなので、注意が必要です。

因果関係なしの相関関係：（３）因果の上流側に共通の要因が存在するケース

では次は、「因果の上流側に共通の要因が存在するケース」について見ていきたいと思います。

これはいわゆる「交絡」と表現されるものに対応するケースです*8。因果グラフで表すとその構造がわかりやすいので、因果グラフを用いて説明していきたいと思います。

「交絡」の状況を因果グラフで書くと以下のようになります：

コトバで書くと「説明変数Aの"上流側"に、説明変数Aと結果変数Zの両者に影響をもたらす要因がある」 というかんじですかね。そのような要因のことを、「交絡要因」とか「交絡因子」と呼びます。上の図の場合、要因B,Cが交絡要因となります。

（ここで要因Cを「交絡要因」と呼ぶかどうかは微妙ですが、要因Cで調整すれば交絡は消える状況にあるのは確かなので一応「交絡要因」として呼んでおきます。ここで要因Cも含めるためのニュアンスとして単に"上流"ではなく"上流側"という表現を使っています）

では、具体的な例で考えてみましょう。河川中の「亜鉛濃度」と「底生生物の種数」の仮想例について見ていきます（この仮想例の生成と解析に用いたRスクリプトは zinc-BOD.R ）。

河川の底生生物の保全のために「底生生物の種数」と河川中の環境汚染物質について調査を行う場合を考えていきます。まず、手始めに重金属濃度データについて調査したところ、河川中の「底生生物の種数」と「亜鉛濃度」に以下のような関係が見られました：

このとき、「亜鉛は底生生物に対して毒性がある」という背景知識を考慮に入れると、「亜鉛濃度の増加→底生生物の種数の減少」の因果関係を想起するのは自然なことかもしれません。ここで「底生生物の種数」に対して「亜鉛濃度」を説明変数とした単回帰をしてみると、その回帰係数は「-1.0」で、傾きの有意差は「p<2.6 x (10の-10乗)」となっています。

はてさてしかしながら：

さらに調査を進めていくと、河川の有機汚濁（BOD）に関しても同様の関係が見られることがわかりました：

はて。これを見る限りは、「BODの増加→底生生物の種数の減少」という因果関係もありえそうです。「底生生物の種数」に対して「BOD」を説明変数とした単回帰を行うと、その回帰係数は「-1.9」で、傾きの有意差は「p < 2x(10の-16乗)」となっています。

はてはて。どちらが「真」の因果関係なのでしょうか？

こういう場合の定石として、結果変数を「底生生物の種数」、説明変数を「亜鉛濃度」と「BOD」とした重回帰分析を行なってみましょう。結果として得られたのは：

> res.num_BOD_zinc <- lm(num_species ~ BOD + zinc)
> summary(res.num_BOD_zinc)
...略
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 18.32240    2.75233   6.657 1.68e-09 ***
BOD         -2.02024    0.16300 -12.394  < 2e-16 ***
zinc         0.08852    0.12811   0.691    0.491    

という結果でした。この結果を読むと、「亜鉛濃度」「BOD」の2つの説明変数による重回帰を用いたところ、「亜鉛濃度」の影響が消えてしまっていることが分かります（対応する偏回帰係数が0.088, p値は0.49）。この結果から、「亜鉛濃度」と「底生生物の種数」の相関関係は交絡によるものであり、「因果関係がないのに相関関係が現れている」ケースであることが分かります。

一方、「BOD」の影響は殆ど変わらず（対応する偏回帰係数が-2.0, p値は < 2x(10の-16乗)）、「BODの増加→底生生物の種数の低下」の因果関係の方は真であると推測できることになります。

上記のような場合には、因果グラフは以下のような構造であると推測できます：

このとき結果変数として「底生生物の種数」を、説明変数として「亜鉛濃度」をとった場合には、「BOD」が交絡要因となっているのが分かります。このような場合には、「BOD」の値で重回帰において説明変数として追加する等による調整を行わないと「亜鉛濃度→底生生物の種数」の因果効果を適切に推測することはできません。

（逆に、上記の因果グラフは、結果変数として「底生生物の種数」を、説明変数として「BOD」をとった場合には、「亜鉛濃度」は交絡要因ではないことを示しています）

一般に、上記の例のように「説明変数Aの"上流側"に、説明変数Aと結果変数Zの両者に影響をもたらす要因」がある場合には、簡単に「因果関係なしの相関関係」のグラフを作ることができます。例えば、説明変数Aと結果変数Zの両者が「年代につれて増加」するような傾向を持つときには、両者の変数の間に因果関係がなくとも、相関関係は簡単に現れます。例えば、このあたりの有機食品と自閉症の例やテレビの台数と死亡率の例などはその典型例といえるでしょう。

このような「交絡」の影響を除去するためにはどうしたら良いでしょうか？

定石としては、交絡の原因となる要因（＝「説明変数Aの"上流側"に、説明変数Aと結果変数Zの両者に影響をもたらす要因」）に基づきデータを層別化して解析を行うか、重回帰においてその交絡要因を説明変数に追加した解析を行うことになります。また、興味のある説明変数以外をまとめてエイヤっと交絡を調整する方法として傾向スコア法なんてのもあります。

ある要因が交絡要因かどうかを判断したい場合には、（厳密なものでなくてもよいので）因果グラフを描いて、その要因が「説明変数Aと結果変数Zの両者に影響をもたらす要因」 かどうか検討してみましょう。

（重回帰における変数選択の際に、「関連のありそうな変数はとりあえず入れとけ」みたいなアドバイスもあるようですが*9、変数を追加することにより逆にバイアスが生じたり*10因果効果の過小推定を引き起こす*11可能性もあり、また余計な変数を入れると推定精度が低下しがちなので注意が必要です。交絡を調整するためにどの説明変数を選択するべきかを判別するformalな基準としてバックドア基準というものがありますので、それについてもいつか書きたいと思います）

また、いわゆる「シンプソンのパラドックス」というのもこのタイプの交絡の一種です。以下の筒井淳也さんの素晴らしい記事を適宜ご参照いただければと思います：
シンプソンのパラドックスの図解 - 社会学者の研究メモ

シンプソンのパラドックスにおいての「どの変数において層別化（重回帰の説明変数として追加）するべきか」という決定は、確率論的考察からは議論することすらできない問題ですが*12、因果グラフさえ描くことができれば、バックドア基準に基づきどの変数が調整すべき変数（＝交絡をもたらす変数）なのかを明確に判別することができます。

＊交絡の問題に関しては、以下の記事でもまとめて書いておりますので適宜ご参照ください：
因果グラフからみる交絡問題：「遺伝統計学における因果問題の特殊性」について考えてみた - Take a Risk：林岳彦の研究メモ

（４）因果関係なしの相関関係：因果の合流点において選抜／層別／調整されてしまっているケース

では最後に、「因果の合流点において選抜／層別／調整されてしまっているケース」について説明していきたいと思います。

このケースは、いわゆる「選択バイアス」と呼ばれることが多いものです*13。因果グラフ系の用語では、「合流点バイアス（collider bias）」と呼ばれます。

因果グラフで見てみましょう：

因果グラフがこのような形のとき、要因Bが「合流点（collider）」と呼ばれるものになります。このような合流点で選抜／層別／調整が起きていると、「説明変数A→結果変数Z」の因果関係がない場合でも、両者のあいだに相関が生まれることがあります。

たぶん合流点バイアスは具体的な例を通して理解したほうがてっとり早いと思いますので、芸術系大学の入学試験の仮想例を考えてみましょう。

ある芸術系大学で入学試験があり、その内訳は「実技試験（500点満点）」と「学力試験（500点満点）」の二つの科目からなるとします。仮想例のための仮定として、実技試験と学力試験の間には本来的に相関関係も因果関係もないものとします。プロットをしてみると、こんなかんじの状態ですね（用いたRスクリプトは exam.R ）：

さて。

ここで、合格基準が「二つの科目の総得点が700点以上」だとします。このとき、諸事情によりデータ解析者には「合格した生徒のデータしか渡されていない」状況を考えてみます。このような状況で「渡されたデータ」からグラフを描き、相関係数を求めてみると：

という関係が現れました。ここでは、実技試験と学力試験の間には、因果関係がないにもかかわらず有意な相関が生じています。これは、データが「二つの科目の総得点が700点以上」という条件であらかじめ選別されていることに起因するバイアスです。図で説明すると：

この上記の青色の線が「二つの科目の総得点が700点以上」を満たすラインになっており、その部分の合格者だけのデータが選抜されていることにより相関が生じていることが分かるかと思います。こういうのがいわゆる「選択バイアス」というやつです。

上記の例のケースを因果グラフで書くと：

のような形になり、「合格の可否（試験の総合点）」が合流点となります。このような合流点の値に基づく選別があらかじめ行われていると、そもそも因果関係のない「実技試験の点数」と「学力試験の点数」の間に相関関係が生じてしまうのです。

（ここで一つ関連して注意しておきたいのは、たとえデータを「ランダムにサンプリングする」といっても、サンプリングの対象となる集団において既にバイアスが形成されている場合には、その"ランダムサンプリング"によってはそのバイアスを除くことはできません。例えば、上記の例の「合格した生徒のグループ」からランダムサンプリングにより50人選んだところで、その合流点バイアス自体が消えることはありません）

上記の例では、得られているデータの形成過程において「合流点での値に基づくデータの選択」が起きていましたが、合流点バイアスは、統計解析の時点における「合流点での値に基づく層別化」や「重回帰の変数として合流点にある要因を追加する」ことによっても同様に生成してしまいます。とっても気をつけましょう。

合流点バイアスを避けるためには、（厳密なものでなくてもよいので）因果グラフを描いてみて、上記のような「因果の合流点」において選抜／層別／調整が行われていないかチェックしてみると良いでしょう。

おまけ：建設的な質問をしよう！：統計解析者への「質問テンプレート」

はい。ここからは「おまけ」です。

上記の４つのケースを踏まえて、因果関係と相関関係に関して統計解析者（プロの研究者を想定）へ建設的な質問をするための「質問テンプレート」を幾つか考えてみたので解説してみます。

(A) 「説明変数Aと結果変数Z、および関連する主要なその他の説明変数（共変量）に関して、あなたが想定している因果グラフを描いてみてください（厳密なものでなくとも構いません）」

これは、上記ケースの２[逆]・３[交絡]・４[合流点]が当てはまっているかどうかをチェックするためのもっとも包括的な質問になります。ここで大事なのは、とりあえずの議論のためには必ずしも厳密な因果グラフを描く必要はなく、変数間の関係がおおざっぱに（上流か／下流か／独立か）議論できれば、多くの場合、間に合うということです。

相手が因果グラフを描いてくれた場合には、上記２・３・４のケースが当てはまりそうか吟味し、また、その因果グラフの理論的／メカニズム的妥当性について議論を進めることにより、建設的なやりとりが期待できます。（ひょっとすると、新たなリサーチクエスチョンが見つかる、なんてこともあるかもしれません）

(B) 「説明変数Aと結果変数Zについて、説明変数Aにおいて異なるサンプル群について、説明変数A以外の要因については同質であると言えますか？その論拠は？」

これも、上記の２[逆]・３[交絡]・４[合流点]のケースに対応する質問となります。この質問をより具体的に発展させると「質問(A)」の形になる、というかんじですかね。

(C) 「説明変数Aと結果変数Zの両者に影響を与えている、共通の原因が存在する可能性については考慮しましたか？」

これは上記３[交絡]のケースに絞って質問する場合になります。要するに「交絡要因についてちゃんと考慮していますか？」という質問です。具体的に頭に浮かんでいる「交絡要因」があるのならば、より具体的に「◯◯が交絡要因となっているのではないでしょうか？」と質問してみましょう。

(D) 「サンプリングもとの集団において既にバイアスがかかっている可能性はありませんか？」

これは上記４[合流点]のケースについて、「合流点バイアス（選択バイアス）があるのではないでしょうか？」という意味ですね。

(E) 「説明変数Aと結果変数Zについて、独立に行われた同様の研究において結果は一貫していますか？」

これは上記１の「偶然によるケース」に対する質問です。

(F) 「あなたの研究デザインは潜在的に多重比較になっていませんか？」

これも上記１の「偶然によるケース」に対する質問です。

あなたが、もし質問される側（統計解析者の側）ならば、質問者から蜂の巣にされないように、これらの点についてはちゃんと研究計画の段階からあらかじめ自問自答しておきましょうね（ニッコリ）。

まとめ：全ては地に足のついた丁寧な考察のために

まとめます。

今回説明してきた「因果関係がないのに相関関係があらわれるケース」は、以下の4つになります：

• (1) 偶然によるケース
• (2) 因果の流れが「逆」のケース
• (3) 因果の上流側に共通の要因があるケース
• (4) 因果の合流点において選抜／層別／調整されてしまっているケース

もしかしたら、上記の4つの他にも（測定や実験デザインの不備等により）「因果なしの相関」が生じるケースもあるのかもしれませんが、それらも結局は上記1-4の形（やその組み合わせで）表すことができるかと思います*14。

で、ここで逆から言うと、われわれがある特定のデータから「相関関係が因果関係を意味するのか」を判断する際には、実際には「上記の4つのケースを除外できるかどうか」を判断しているとも言えるわけです。つまり、様々な角度からの検討により、上記4つのケースを高い確度で除外できると判断できるような場合には、「因果関係があって、相関関係が生じている」と推論できる、ということになります。

はい。

で、

まあ。なんといいますか。詰まるところをいえば（ヒューム的な意味で）我々は因果関係を知ることができない、というのは真だとは思うのです。完璧な文章や完璧な絶望が存在しないように、完璧な因果推論などといったものは存在しないのです。

ですが、でもそこで「因果なんてけっきょく分かんないんだから統計データから因果の考察するのなんて無意味っつうか無粋だよね〜」とか「因果なんてけっきょく分かんないんだからなし崩し的に因果があるって解釈しちゃってOKですよね〜」みたいな雑なかんじになるのではなく、そのつどそのつど「データ」と「理論／メカニズム」と「現場知」を照らし合わせながら「因果の有無と程度」について地に足をつけて丁寧に考察していくのが大事かな、って思う今日このごろであります。

ほんで

今回の記事がもしそのような「地に足のついた丁寧な考察」の一助となれば本当に嬉しいです。

＃次回は、今回の記事とは逆の「因果関係があるのに相関関係があらわれないケース」についてまとめたいと思います。
＃あといくら何でも記事が長すぎですよね。本当にすみません。。。 ＞ 読者さまがた
＃尚、今回の記事は以下のtogetterに触発されて書かれたものです。まとめ人および登場人物の方々に感謝＆リスペクト申し上げます。
相関関係と因果関係をごっちゃにしないために - Togetterまとめ

関連図書など

統計的多重比較法の基礎

多重比較についてはこれで勉強しました。

メタ・アナリシス入門―エビデンスの統合をめざす統計手法 (医学統計学シリーズ)

メタアナリシスについてはこれから勉強したいです。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

統計的因果推論の鉄板本その１（Rubin系）

統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)

統計的因果推論の鉄板本その２（Pearl系）

Excelで学ぶ共分散構造分析とグラフィカルモデリング

上の宮川本が難しい方は、こちらからどうぞ

統計的因果推論 -モデル・推論・推測-

Pearl師匠の本。いろんな意味で激ムズ。

Modern Epidemiology

この本（3rd ed）の12章の「Causal Diagrams」の章が（それなりに難しいが）英語が読めて、疫学がベースの人には一番有益な因果グラフ入門かもしれない。

関連する過去記事

「美人ほど女の子を産む」はウソ？：A. Gelmanによる統計的欠陥の指摘のメモ - Take a Risk：林岳彦の研究メモ
無から有（意差）を生む：多重比較でウソをつく方法 - Take a Risk：林岳彦の研究メモ
確率と因果を革命的に架橋する：Judea Pearlのdo演算子 - Take a Risk：林岳彦の研究メモ
因果グラフからみる交絡問題：「遺伝統計学における因果問題の特殊性」について考えてみた - Take a Risk：林岳彦の研究メモ
統計的因果推論（傾向スコア）の勉強会資料をアプしてみた - Take a Risk：林岳彦の研究メモ
発表資料アプ『相関と因果について考える：統計的因果推論、その(不)可能性の中心』 - Take a Risk：林岳彦の研究メモ

単なる宣伝

私が寄稿させていただいた文芸誌が中野のタコシェで買えるようになったようです→ http://t.co/FzHMI9siWJ

*1:本当に名作だと思う

*2:"単なる偶然"とは何か、と問い詰められるとちょっと困るけど

*3:どんな場合やねん

*4:暇なひとはRのスクリプトで実際に試してみてください

*5:実験データの場合には再試験してください（キッパリ）

*6:いわゆるメタ・アナリシス

*7:フロントドア基準的な方法で明らかにできる可能性はある、かも？

*8:「交絡」の定義については、研究分野間で微妙にちがったりするので、この表現がピンとこない方々ももしかしたらいるかもしれませんがご容赦を

*9:まあ応対している相手の統計リテラシーのレベルにもよる話なので一概には否定はしませんですよ

*10:合流点によるバイアスなど

*11:中間変量による調整の場合など

*12:do演算子のような確率と因果を架橋する概念が必要となる

*13:この用語法も分野によるらしいので「そんなの選択バイアスじゃないやい！」と思われる方々もいらっしゃるかと思いますがご容赦くださいね

*14:たぶん

宣伝（文学フリマ）：自主制作文芸誌『線と情事』に詩を書かせていただきました

あまりここにはわたくし事を書かないようにしているのですが、今日は諸事情により宣伝させていただきます。

このたび、何の因果か、甘茶茂＠NECOfanさんが責任編集を務める文芸誌『線と情事』へと詩（12ページ、7篇）を寄稿させていただきました。

＊文芸誌『線と情事』の内容は以下をご覧ください：
コミックナタリー - 島田虎之介、しまおまほ、関根美有ら自主制作文芸誌に

シマトラさんやしまおまほさんなどの豪華執筆陣に紛れて私の名前がエントリーされているのは、なんというか「21世紀枠」的なかんじは否めませんが、文学フリマなどへお出かけの方がおりましたら、ぜひご購入のほど何卒よろしくお願いいたします。

（個人的に私はタマフルTHE MOVIEのDVDを買うくらいにはコアなタマフルファンなので、しまおさんと名前が並んで載っているだけで鼻血が出るくらいうれしいのです。あとザ・キャプテンズの傷彦氏を初めとする個人的にリスペクトしてやまない仙台時代の知人たちとこんな形でまた「対バン」できるのは感慨深いものがあります）

とりあえず大阪の文学フリマ（4/14）を皮切りに、東京の文学フリマおよびweb通販などで販売される予定です（詳しくはこちら）。好事家の皆様、ぜひよろしくお願いいたします。

関連図書：

東京命日

ガールフレンド (P‐Vine BOOKs)

THE CAPTAINS ANTHOLOGY