こんにちは。林岳彦＠はてなジェシ（フロンターレの生命線）です。
さて。
先日、某勉強会で論文（Lim et al. 2012 in Lancet）の紹介をしたのでついでに発表スライドをアプしてみました。論文の内容は、簡単に言うと「世界の健康リスク要因のランキングを作ってみた」というものです。
詳しい説明は以下のスライドをご参照ください。図表が小さいですが、細かいところまで見たい方はリンク先よりスライドのPDFをダウンロードして拡大して見ていただければと思います。あくまで口頭での説明を前提とした発表ファイルなので、これだけを見てもわけが分からないかもしれませんが、例えばFig2とかFig5とかのリスク要因のランキングを見るだけでも「へえ」というかんじで面白いかもしれません。
ちなみに元論文が載ったGBD study特集号はこちら↓
The Lancet : Volume 380, Number 9859, 15 December 2012
お久しぶりです。林岳彦です。もうすぐ『愛なき世界』の日、いわゆる（マイブラッディ）バレンタインデーですね。何かと雑音が多いこの世界ですが、いつでも自分の足元を見つめて行きましょう。
さて。
今回は、以下の：
のあたりの皆様の良記事に触発されて「仮説検定」について何か書いてみようと思いました。で、書こうと思えば色々な側面から書ける気もするのですが、今回はちょっと斜めからのアプローチとして、「リスク分析の人の頭のなかで仮説検定はこんな感じに見えている」というところを書いていきたいと思います。
ここで、ひとくちに「リスク分析」と言っても思い浮かべるところは人それぞれかもしれませんが、今回は以下の：
では、つらつらと書いていきたいと思います。
（今回もとても長いです。いつもすみません。。。あと、今回とくに一部内容に余り自信がないので*2、もし間違ってたら適宜ツッコミおねがいいたします）
さて。
もしかしたら意外に思う方もいるのかもしれませんが、「リスク分析」においては仮説検定を使うことは殆どありません。私も一応リスク分析に類する論文を何報か書いていますが（論文へのリンク）、それらの論文の中で仮説検定を用いたことは一度もありません*3。
これは単に私が個人的に仮説検定を使わないという話ではなく、例えば、冒頭に挙げたDavid Vose著の『入門リスク分析』の中でも仮説検定の話はほんの僅かしか出てきません*4。これは、そもそも「リスク分析」のツールボックスの中に「仮説検定」が基本的に含まれていないことを意味しています。
では、なぜリスク分析においては仮説検定が殆ど使われないのでしょうか？
リスク分析において仮説検定が殆ど使われない理由としては、本来的なものと技術的なものの２通りの理由があります。
まずは、本来的な理由の方から説明してみましょう。
本来的に、「リスク分析」とは、リスク（＝影響と不確実性の在りよう）を推定し、意思決定者の参考となる情報を提供することが基本的な目的です。つまり、「意思決定支援」こそがリスク分析の本懐であるわけです。
ここで、「仮説検定」がリスク分析と本来的に相性が悪いところは、仮説検定が「意思決定の自動代行機能」の側面を持っているところです。つまり、リスク分析の目的が「意思決定支援」であるにもかかわらず、仮説検定は仮説の採択に関して論点を先取りして肝心の「意思決定」をオートマティックに代行してしまう*5側面があり、この性質は「リスク分析のツール」としてはとても悪質なものであると言えます（→参考：具体的な弊害についての過去記事）。
もう一つの、技術的な理由の方は、リスク分析では推定の不確実性を見るときに基本的に「不確実性の分布全体を見る」ので、仮説検定がそもそも不要であることです。
さて。
本記事では以下、この「不確実性の分布全体を見る」ということについて、例を交えつつ説明していきたいと思います。
では、「不確実性の分布全体を見る」とはどういうことなのかを見て行きます。
説明のための仮想例として、「ブリの養殖業」における、ブリをより大きく育てるための「新エサの開発」についてのケースを考えてみましょう。
あなたはブリの養殖業を営んでおり、ブリを大きく育てるための新エサを開発しているとします。そして実際に、その新エサの効果を調べるために、新エサを用いて1群のブリを育てているとします。
ある日、新エサで飼育した1群のブリの中から、30匹をランダムに選んでその体長を測ってみた結果、以下のヒストグラムが得られました（以下は対応するRスクリプト）：
# 30匹のサンプルの体長データ buri_length_30samples.v <- c(93.603192838865, 105.754649863331, 90.4655708138493, 126.929212032067, 107.94261657723, 90.6929742382298, 110.311435786427, 114.074870576938, 111.636720274802, 98.4191741926547, 125.676717526763, 108.847648546171, 93.6813912918729, 69.7795016923375, 119.873963772147, 102.325995864772, 102.757146053516, 117.157543160279, 115.318317926471, 111.908519818263, 116.784660574123, 114.732044510966, 104.118474750478, 73.1597245620494, 112.297386218421, 102.158068907065, 100.66306739942, 80.9387142415109, 95.8277491733707, 109.269123402996) # 散布図をプロット plot(buri_length_30samples.v,col="royalblue",lwd=3,ylim=c(60,140)) # ヒストグラムをプロット require(MASS) truehist(buri_length_30samples.v,h=5,col="royalblue",prob=F)
この図は、サンプルした30匹のブリの体長のヒストグラムを表しています。この図から、体長には70cm〜130cmほどのバラツキがあり、平均は100cm〜110cm程度であることが見て取れます。
ここで話を単純化するため、「従来のエサ」で育てた場合の平均体長が「100cm」であることが既に知られていると仮定します。（また、エサの違いは体長の分散に影響を与えないと仮定します）
このとき、「従来の餌と比べた場合」の新エサの体長への影響の指標は、「新エサで育てたブリの体長」から、従来の餌による平均体長である「100cm」を差し引くことにより得ることができるでしょう。（単に、上図の横軸が100cmぶん左側へずれるだけですが）一応、そのヒストグラムを描いてみると：
# 新エサの体長データから100引いたものをweight_diff.vに格納 weight_diff.v <- buri_length_30samples.v - 100 #「従来の餌との差」をプロット truehist(weight_diff.v,h=5,col="royalblue",prob=F) # 体長変化の平均値を算出 mean(weight_diff.v)
となります。図を見ると、分布の全体はゼロよりも若干右側に重心があるようにも見えます。「ゼロより右側が多い」ということは、従来のエサに比べて、「新エサは体長をプラス方向に変化させる傾向がある」ということになります。
ここで、上の30匹の体長の変化分の平均を計算すると「+4.23cm」という値が得られました。ここでもういっそのこと「新エサの方が従来の餌よりも体長を大きくする」ということでOK・・・と結論してしまいたいところですが、これはあくまでも30匹のサンプルから得られた「サンプルの平均値」であるため、この程度の差はたまたまの偶然により生じている可能性も否定できません。ここで、本当に新エサが「どのていど体長を増加させると期待できるのか」を判断する際には、限られたサンプルから「サンプル元の集団における平均値」を推定する際に生じる「不確実性」を考慮する必要があります。
（以下、煩雑さを避けるために「サンプル元の集団における平均値の推定」のことを単に「平均値の推定」と略記していきますが、推定している対象はあくまでもサンプル元の集団における平均値であることにご留意ください）
さて。このような、限られたサンプルからサンプル元の集団における値を推定する際の「不確実性」を取り扱うためのよく知られた方法の一つは仮説検定です。しかし、今回はいわゆる「リスク分析」っぽいアプローチを用いてこの「不確実性」を取り扱ってみたいと思います。
はい。
とは言っても、そんなに特別なことをするわけでもありません。
確率論的リスク分析においては、「推定に伴う不確実性」を可能なかぎり常に「分布の形」で取り扱っていきます。実務的には、そのための手法には大きく分けて「ブートストラップ法」と「ベイズ法」の２つのやり方があります。今回は、比較的とっつきやすい手法として「ブートストラップ法」の方を使っていきます*6。
ブートストラップ法というのは：
という方法です。（なんのこっちゃかよく分からない方は、基礎統計学講座＠ウィキやフリーソフトによるデータ解析・マイニングや奥村晴彦さんの記事のあたりの良記事をオススメいたします）
まあ案ずるより産むが易しということで、とりあえずやってみます*7。対応するRのコードは例えば：
### ブートストラップ（反復数50000）による平均値の推定分布の算出 # ベクトルの初期化 mean_bootstrap.v<- numeric(50000) # ブートストラップ（反復50000回） for(i in 1:50000){ # "新”サンプルを最初のサンプルからリサンプリング bs <- sample(weight_diff.v,30,replace = TRUE) # “新”サンプルの平均値を算出しmean_bootstrap.vに格納 mean_bootstrap.v[i] <- mean(bs) } # 平均値の推定分布を描く truehist(mean_bootstrap.v,h=0.5,col="royalblue")
となります*8。このブートストラップの計算では、「ブリの体長データそのもの」ではなく「従来のエサと比べた場合の変化量（つまり「体長データ」から100cmを差し引いたもの）」を用いています。
ブートストラップの結果として、以下のような「従来のエサと比べた場合の変化量」の推定された平均値の分布（ヒストグラム）が得られます：
はい。この分布が、「新エサが体長に与える変化の平均値の推定に伴う不確実性」を「分布の形」で表したものになります。確率論的リスク分析においては、推定に伴う不確実性をこのような分布の形で扱うことが一般的です。
（一般論として、このような分布が何を意味しているのかというと、例えば、この分布の「幅」は推定の「精度」を表していることになります。もし30匹よりも少ないサンプルに基づきこのような分布を描けば「より幅の広い分布」が得られますし、30匹よりもずっと多数のサンプルに基づけば「より幅の狭い分布」が得られることになります）
さてさて。では、リスク分析においてこのようなやり方で不確実性を扱うことのメリットとは一体何でしょうか。その主なメリットは、「推定の不確実性」の二次利用のしやすさにあります。
はい。
推定における不確実性を上記のような「分布の形」で保持しておくと、その分布データを意思決定支援のための新たな解析に簡単に供することができるというメリットが生じます。
・・・と言われてもにわかには何のことか分からないとは思いますので、上記の解析を用いて、意思決定において非常に重要な項目であると思われる「ブリの売価」がどう変化するかについての新たな解析例を展開してみましょう。
解析例としての仮定として、生産したブリの売値（100匹当り）は『「ブリの平均体長の3/2乗」の1000倍』、つまり以下の式*9：
「売値」 ＝ 1000 × （（ブリの平均体長） ^ （3/2））
で表せることが知られているとします*10。
従来のエサで育てた場合には、仮定により平均体長は「100cm」なので、この式より「売価」はちょうど 1000 × (100 ^ (3/2)) = 1,000,000円（100匹当り）になると計算できます。
同様に、新エサを与えた場合の「売値」は：
「売値」 ＝ 1000 ×（（100 + 新エサによる体長変化の平均値) ^（3/2））
と計算することができるでしょう。
ここで、「新エサによる体長変化の平均値」の推定の不確実性を前述のような「分布の形」で保持している場合には、その分布を表す数値データ群を「新エサによる体長変化」の値としてそのまま代入してやることにより、体長変化に関する推定の不確実性の影響込みでの「新エサで育てた場合の売値」の予測分布をすぐに算出することができます。
試しにRで計算してみましょう：
# ブートストラップによる推測分布を代入しつつ売値の予測分布を算出する price.v <- 1000 * ((100 + mean_bootstrap.v)^(3/2)) # 売値の予測分布のヒストグラムを描く truehist(price.v,h=5000,col="orange") # 売値の予測の平均値 mean(price.v) # 売値が1,000,000円となるパーセンタイルの算出 ecdf(price.v)(1000000)
計算の結果として、上図のような「新エサで育てた場合の売値」の予想分布が得られます。
この予想分布がどのように「意思決定を支援」しうるかを見てみましょう。
まず、上図の予想売価の平均値を算出すると「1,064,406円」となり、従来の売価が「1,000,000円」なので、新しいエサによる100匹当りの追加コストが「＋64,406円」のところが（新エサの効果推定に伴う不確実性の下での）新エサ導入の"損益分岐点"になると計算できます。
また、売価が従来のエサと同じ「1,000,000円」となるのは分布の4.8パーセンタイルのところであるため、新エサの効果推定に伴う「不確実性」を考慮に入れた場合には、売価が従来のエサのものを下回ってしまうリスクも依然5％程度はあることがわかります。
はい。
こんなかんじで統計的推定に伴う不確実性の分布を「二次利用」していくのがいわゆる「リスク分析」っぽい解析アプローチとなります。このように、もともとの「データ（今回の場合は新エサによる体長データ）」のコンテクストを、より意思決定に直接的に関連するコンテクスト（今回の場合は「売値」の予測分布）へとできるだけ定量的に繋ぐための解析を行うことがリスク分析の基本的なミッションになります*11。
さて。
とは言っても、場合によってはリスク分析においても（主に説明のための便宜として）、「仮説検定」や「信頼区間」のような計算が求められることもあります。次は、不確実性の「分布」がそれらの古典的な統計的概念とどういう関係にあるのかを見ていきます。
不確実性の「分布」と「仮説検定」「信頼区間」の関連性は、図で見ると分かりやすいので図に描いてみます：
この図は上述のブートストラップ法で算出した「新エサによる体調変化の平均値」の推定値の分布に、「仮説検定」や「信頼区間」との関係の説明を書き込んだものです。この図を使って説明をしていきます。
まず、上図(B)の場合である「信頼区間」を見ていきます。推定値の分布と信頼区間との関連性はシンプルで、いわゆる95%信頼区間というのは推定値の分布の2.5パーセンタイル*12と97.5パーセンタイルの値に挟まれる範囲を意味します。
実際にRで推定分布から計算すると：
# 95%信頼区間の計算 quantile(mean_bs.v,p=c(0.025,0.975))
2.5% 97.5%
- 0.8022501 8.9183979
として95%信頼区間 [-0.80 〜 +8.9 ]を簡単に計算することができます。（このようにブートストラップ法で計算された信頼区間を特に「ブートストラップ信頼区間」と呼ぶこともあります→参考：基礎統計学講座＠ウィキ, hoxo_mさんの記事）
また、上図(A)の仮説検定との関連に関しては、上記の95%信頼区間が「ゼロをまたぐ」ことから、「新エサによる変化の平均値はゼロと異ならない」という仮説を棄却することができない、という5%有意水準を念頭においた仮説検定的な判断を行うことも可能です。
（＊補足：ここでなぜ仮説検定”的”という言い方をわざわざしているかと言うと、いわゆる「仮説検定」は帰無仮説が正しいという仮定のもとで生成された統計量の分布のもとで、実際に得られているサンプルからの統計量を評価するというロジックで行うので、上記の例のように（帰無仮説ではなく）実際に得られているサンプルの統計量の分布をもとにした判断を「仮説検定」と呼ぶのは適切ではないと思われるからです*13。ただし、帰無仮説のもとでの統計量分布に基づき考えるのも、得られているサンプルからの信頼区間に基づき考えるのも、見ようとしている先のものは実質的には殆ど同じです（が、比較対象先の値自体もまた推測値である場合には、「2つの推測値の分布」の重なりを考える必要があるので注意が必要です→参考：「馬車馬のように」の記事, kamedo2さんの記事）
同様に、もし「新エサによる変化の平均値がゼロより大きいと言えるかどうか」の判断を5%有意水準を念頭においた（片側）仮説検定的に行う場合には、5%信頼下限値（＝分布の5パーセンタイル値）より上の区間がゼロを含むかどうかにより判断を行うことも可能です。（今回の例では推定値の分布の5パーセンタイル値は「0.039」というゼロより大きい数値になり、5%信頼下限値より上の区間にはゼロを含まないことになるので「新エサによる変化の平均値はゼロよりも大きい」と判断することができます）
はい。
このように、解析を通じて不確実性を「分布の形」で保持してさえいれば、信頼区間や仮説検定的な解析もニーズに応じて簡単に行うことができます。このため、リスク分析では不確実性を「分布の形」で取り扱うやり方がデフォルトとして好まれる一方で、仮説検定に関するエトセトラには余り本来的に興味が持たれない、ということになります。
（喩えて言うと、不確実性を「分布の形」で取り扱うことを志向することは、魚を「一本まるごと」取り扱うことを志向することに似ていると言えます。魚を「一本まるごと」保持しておけば、必要に応じて「信頼区間（さく切り）」や「仮説検定（切り身）」の形に加工することは容易ですが、不確実性を最初から「信頼区間（さく切り）」や「仮説検定（切り身）」の状態にしてしまうと、それらに向いた用途以外への二次加工（あら煮とか）を行うのが難しくなり、意思決定支援のコンテクストに沿った解析の可能性の幅をあらかじめ大きく狭めてしまう危険性があるのです）
では、今回の記事の内容をまとめます：
はい。こんなかんじになるかと思います。
で、少し補足しておきますと、この記事を通して私が言いたいのは「何がなんでも不確実性を扱うときは分布で扱え」ということではありません。
どちらかと言うと、「不確実性をどのように扱ってもいいけれど、アタマの中のイメージとしては常にどっかで「分布の形」を想像しておくべし」ということが伝えたいことです。なんというか、不確実性について「分布のイメージ」を持っている人がそれを「仮説検定」や「信頼区間」の言葉に変換するのは容易ですが、仮説検定のp値しか気にしてない人がそれを分布のイメージに翻訳するのは難しいので、実際の解析法として何を使うかはともかくとしても、「分布のイメージ」を常にアタマの中にもっておくことが大事かと思います。
実務的にも、例えば、リスク分析においては「分布のかなり端っこの部分が特に重要」だったり「推定分布が多峰型」だったりすることもありますので、そんなときに「分布の形」を無視して信頼区間だけで解析しようとすると痛い目にあうので注意が必要です。このような場合も不確実性についての「分布のイメージ」がアタマに入っていると地雷を避けることができるかと思います。
で、ついでにもう一件補足しておきますと、「意思決定支援」が本懐であるのは何もリスク分析だけではなく、ビジネス畑の「データサイエンティスト」の方々の中にもそのような方々が多くおられるのでないかと推測しております。さいきんわたくしも
![データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)](http://ecx.images-amazon.com/images/I/51jDJwFEw8L._SL160_.jpg "データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)")
データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)
という良書を拝読させていただき大変勉強になったのですが、この中に今回示したような「リスク分析」、つまり「データ解析」の結果を二次利用してさらなる解析を行い、意思決定により密接に関連するコンテクストにおける指標への影響を提示しうるような、そんな解析アプローチも載っていても良いのではなかろうか、という感想を抱きました。
スーパーで出来合いの「切り身」を買ってくるようなデータサイエンティストではなく、クライアントのニーズに沿って魚を一本まるごといかようにでも捌ける「寿司屋の大将」のようなデータサイエンティストを目指すのならば、そんな「リスク分析という包丁」も一つ持ってると身を助くこともあるのかも、とか思ったりする今日このごろです（←無責任）。
で、もし、本記事がそのような「寿司屋の大将」への入り口となりうればとても嬉しいと思っております。
# 「リスク分析」をもっと知りたい方にはこちらをオススメいたします↓
仮説検定関連：
ブートストラップ関連：
信頼区間関連：
Rによるベイズ解析（MCMC）の実装関連：
あんまり関係ないかもですがリスク分析における「変動性」と「不確実性」について自分で昔に書いた記事：
*1:なので、もし読んでいて違和感がある場合には、「リスク分析」を「確率論的リスク分析」に適宜置き換えて読んでいただければと思います
*2:特に古典的な「統計学」に関するところが苦手です
*3:文献からの既往データとして仮説検定からの結果を（仕方なく）使うことはあります
*4:本書の全567ページ中で仮説検定の話が出てくるのは僅か10ページほどで、しかも内容は分布モデルの選択の文脈での適合度検定の話なので、本来なら仮説検定よりもAICを使うべき部分になります
*5:あるいは、代行しうると多くの人にみなされている
*7:まあ30匹程度でブートストラップやるのは精度が低いかもですが今回はその辺はあまり気にせずにやっていきます
*8:実装のコードはいろいろ方法はありますが、今回はsample関数を使いました。分布の形状が綺麗になるように多めに50000回の反復をさせています
*9:すみません。この式はテキトーに考えたもので特に意味はありません
*10:説明の単純化のため、ブリの体長の分散の影響と、この関係式自体に関する不確実性は無視できるものと仮定します
*11:今回の例は説明のためかなり単純化していますが、実際には多数のコンテクストからの多数のパラメータが絡み合う複雑な解析になることももちろんあります
*12:xパーセンタイル = 小さい方から数えて全体のx%目の順番に相当する値
*13:あんまりこのへん自信ない
こんにちは。林岳彦です。はてなジョシュ（バーネット）です。今回から「はてなブログ」へ引っ越しました。今後とも引きつづきよろしくお願いします。
さて。
前回までの記事では、実験データではない調査観察データを用いた因果効果の推定における注意すべきバイアスの類型について書いてきました。
ここでなぜわざわざ「実験データではない」という但し書きをつけているのかというと、適切なデザインに基づき行われた実験（もしくは介入を伴う調査）からのデータは、処理・条件の違いによる結果の差を素直に「因果効果」とみなして解釈できるので、余り細かいことを考えなくても大丈夫だからです*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人の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)と(IV)の「平均年収データ」については、特に問題はないかと思います。文として素直に意味が通ってますし、実際にデータを得ることも（少なくとも原理的には）可能です。
一方、(II)(III)の「平均年収データ」は、原理的に「そもそも観測することができないデータ」となっています。「職業訓練を受講したグループ」についての「職業訓練を受講しない場合の5年後の年収」はそもそも観測することができませんし、同様に「職業訓練を受講しないグループ」についての「職業訓練を受講した場合の5年後の年収」も観測することができません。
(II)(III)のようなデータは「事実に反した条件（counterfactual condition）」の下でしか得られないデータとなります。つまり、ここでは：
「職業訓練を受講しないグループ」が”もし職業訓練を受けていたとしたとき”の、5年後の平均年収」
という形で、「（事実としては受講していないけれども）もし職業訓練を受けていたとしたとき」という、反事実的な状況が想定されていることになります。
はい。
ここで重要なポイントを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年後の年収データとしての
の期待値は
の期待値と一致すると考えられます。なぜなら、両グループは単にコイントスで分けただけなので、（もし値が測定可能であった場合に）両者の値（の期待値）に差が出ると考えるべき理由が何もないからです。（例えば、両者のグループの「身長の平均値」の期待値が異なるとは考えられないですよね。それと理屈的には同じです。また、数式を用いた補遺を記事末尾に載せましたので適宜ご参照ください。）
もちろん実際には(III)のデータは反事実的なので得ることはできません。しかし、無作為化（コイントス）によりグループを分けた場合には、(III)と(IV)の「平均年収」の期待値は理屈的に一致すると考えられるので、ここで「(III)のデータの身代わりとして(IV)のデータを使う」という"迂回路"が使えるようになるわけです。
この"迂回路"を使うことにより、
(I) 職業訓練を受講したグループの「職業訓練を受講した場合の5年後の平均年収」
ー (III) 職業訓練を受講したグループの「職業訓練を受講しない場合の5年後の平均年収」
という不可知な項を含む"因果効果”を
(I) 職業訓練を受講したグループの「職業訓練を受講した場合の5年後の平均年収」
ー (IV) 職業訓練を受講しないグループの「職業訓練を受講しない場合の5年後の平均年収」
と読み替えることが可能になるわけです。これらの(I)と(IV)の量は観測可能なデータなので、この読み替えを行うことにより、「職業訓練→将来の年収」の"因果効果”を推定することができることになります。
はい。
これで、『因果推論の根本問題』を統計学的に解くことができました。（これが、いわゆる『統計学的因果推論』というものの一例となります）
（ただし、ここで推定できたのはあくまで集団への「"平均的"因果効果」であることに注意してください。そもそもの話である、職業訓練が「親友K」に与える因果効果は、やはり不可知なままなのです。たとえ職業訓練の"平均"効果がプラスであったとしても、親友Kにはそれがマイナスに働く可能性は常に残ります。個体に対する因果効果は、本質的にやはり不可知なのです。）
ではまとめます。
今回は以下のようなことを見てきました：
はい。
一見地味に感じるかもしれませんが、「ランダムネス」の意図的な利用により「法則」を浮かび上がらせるというこの「無作為化というアイデア」はなかなかに凄いものであるようにも思います。かのリチャード・ドーキンスが「ダーウィン以降で最も偉大な生物学者」としてフィッシャーの名を挙げていましたが、この「無作為化というアイデア」はその「ダーウィン以降で最も偉大な生物学者」が遺したものの中でも、「最も偉大な発明」の一つと言えるかもしれません。
そんなわけで、「統計的因果推論」の分野で有名な研究者といえば、RubinやPearlの名が挙がるのかもしれませんが、「統計的因果推論のキング」を一人選ぶとしたらやはりフィッシャーになるのではないかと思います。
はい。
というわけで、『統計学の道を歩むかぎり、お天道様とフィッシャーはどこへいってもついてまわる』という格言を再認識したところで、今回の記事も終わろうかと思います*3。
＃次回からは、そろそろ「確率概念」についてガッツリ書いていこうかなあ、と思ってます。
上記の説明では、「コイントス」を用いて無作為化を行っていますが、別にそこは必ずしも「コイントス」でなくとも構いません。大事なことは、割付（＝各人をどのグループに分けるのか）を、「受講した場合／しない場合の5年後の平均年収」と本来的に全く関係ない（＝独立な）要素にもとづき決めるということです。
数式で書くと、上記の表の
の期待値は、それぞれ
と書けます。ここでE(A|B)という表記は、「Bという条件」の下での「Aの期待値」を表しています。
ここで、「B」の条件がAと全く関係ない（＝独立な）要素によって決められている場合には、E(A|B)=E(A)となるため：
となり、「(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。
正直どうにもならないことも多いですが、とりあえず「傾向スコア」とか「バックドア基準」とかを参照してみてください。
統計的因果推論（傾向スコア）の勉強会資料をアプしてみた - Take a Risk：林岳彦の研究メモ
発表資料アプ『相関と因果について考える：統計的因果推論、その(不)可能性の中心』 - Take a Risk：林岳彦の研究メモ
===【20180806追記ここから】======================
久しぶりに読み返したところ、上の記述に対しては「サッパリすぎるだろ」と思ったので追記しておきます。無作為化が成立しているか怪しいときは「まずバックドアパスを閉じれるように工夫してみて、それでも閉じれているかどうか分からないときには想定されるバックドアパスの「太さ」について定量的に考察してみましょう（感度分析など）」となるかなあと思います。簡単にあきらめないで！
＊とりあえず以下の良記事を参照すると吉：
www.krsk-phs.com
www.krsk-phs.com
===【20180806追記ここまで】======================
詳しくは、以下の参考文献を読むのをオススメいたします。
今回の記事は以下の二つの名著の中の説明をmixtureしたものとなっております。ぜひこれらの書籍もご参照いただければと思います。
統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)
このブログのデザインは、はてな提供の"epic"のデザイン（背景を白に変更）を利用したものですが、「見出し」のデザインとサイドバーの「人気記事」の部分だけカスタマイズを加えています。
「見出し」のスタイル変更は以下のサイトを参考にさせていただきました：
【CSS初心者向け】ブログで良く見かけるシンプルな見出しの作り方 | delaymania
はてなブログのカスタマイズの方法についてまとめてみた【ま】 - はてブのまとめ
「人気記事」のサイドバーは以下の記事のものを利用させていただきました：
【修正済】コピペで簡単！ はてなブログの人気記事を画像付きで表示させる方法。 - #ChiroruLab
大変ありがとうございました。リスペクト＆感謝申し上げます。
どもっす。林岳彦です。先日、某所で統計解析の講師役をしました。その際に解析環境の準備の手間を省こうと思って、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、媒介変数法（フロントドア基準）を用いて因果効果を推定できたりするという利点もあります）
では、今回のシリーズ全体を軽くまとめてみたいと思います：
今回の一連の記事では「因果関係があるのに相関関係が見られない」以下のケースを見てきました：
はい。いろいろ書きました。お付き合いいただき誠にありがとうございました＞読者さまがた
今、自分はジャワでリゾートを満喫しながらこの記事を書いているのですが*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。
はい。
現場からは以上です。
＃あ、次回くらいから「はてなブログ」に引っ越そうかな、と思ってます（なんとなく）。
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)
統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)
多変量解析の展開―隠れた構造と因果を推理する (統計科学のフロンティア 5)
*1:結局みなさまにRStudioのデスクトップ版をインストールしていただきました。後で見なおしたらAmazon EC2のマイクロインスタンスではそら落ちるわ、って話ですわ。想定が甘すぎ＞俺
*2:統数研あたりが全国の迷える統計講師役のために登録制で使えるぶっといRStudioサーバー版を立ててくれたりしないかなあ（妄想）。。統計実習って実習時間のかなりの部分が皆の解析環境のセットアップに割かれがちなんすよ。。
*3:10000人のデータ全部を描くとゴチャゴチャしすぎるので最初の500データだけプロットしてます
*4:あるいは、めんどくさいオッサンに絡まれがちな
*5:もちろん、日常的な用法はTPOなどに応じてかなり幅があるものなので、これとは違う言い方も同様に成り立ちうるかと思います。
*6:パス分析における「塩分摂取量→超過死亡リスク」の「総合効果」を見ていることになります
*7:パス分析における「塩分摂取量→超過死亡リスク」の「単独効果」を見ていることになります
*8:例えば、今回の解析例は、「血圧」自体も超過死亡リスクを下げるための効果的な介入ポイントであることを示しています
*9:嘘です。本当は午前5:45のマクドナルドで書いてます
*10:実験計画法に基づいた方法により得られたものではないデータ
*11:あるいは、対処不可能な場合には、その対処不可能性が明晰に理解できるようになる
*12:その意味で私は、一ノ瀬正樹氏のベイジアンネットワーク関連の理論に対する評価にはかなり不満があります。Pearlの理論のキモをちゃんと理解してるのかなあ？
どもです。林岳彦＆オメガトライブです。きみは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最小のモデル」を見つけたとしても、その"ベストモデル"の偏回帰係数の値が「各変数の因果効果を最も適切に反映している」とナイーブに前提するのはやめておきましょう。（特に、解析のそもそもの目的が「介入」にある場合はなおさらです！→過去記事）
はい。今回の「中編」をまとめます。
今回は「因果関係があるのに相関関係が見られないケース」として：
というケースについて見てきました。また、その中で：
ということにも触れました。
さて。
今回の合流点バイアスの件ですが、バイアスロン／トライアスロンの例のように変数間の質的な関係性（＝背景の因果構造）が分かっている場合には、合流点バイアスを避けることは比較的簡単かもしれません。しかし、背景の因果構造が不明な場合にそれを避けることはなかなかに難しいものです。
良かれと思って加えたもので逆に墓穴を掘る、というのは人生でもよくあることですよね。
一方、何かをやっているうちに「そもそもの目的がよく分からなくなってくる」というのも人生ではよくあることです。
次回の「後編」ではそのような例として、「中間変量によるマスク」のケースについて見ていきたいと思います。
今回の話をformalな形で理解したい方は以下の本をどうぞ：
統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)
*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:前者は本質的に「確率」のレイヤーの話であり、後者は本質的に「因果」のレイヤーの話になります
どもお久しぶりです。林岳彦です。ローソンなどで売ってるいなばのタイカレーはそうめんのつけ汁として使ってもマジうまいのでオススメです。
さて。
今回は前々回の記事：
因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） - 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とか）にして考えていく方法が一般的となっています。
この辺で、いったんまとめてみます：
というかんじになるかと思います。（ここまではよろしいでしょうか？）
さて。
ここまで、「検定力」は「サンプルサイズ」「効果量」「有意水準」に依存することを見てきました。仮説検定の枠組みにおいて、これらの4つの量は常にお互いに依存し合う"カルテット"的な量になっています。そして、これらの中でも、多くの場合に実験や調査の前にわれわれが主体的に大きく変えられる余地があるものは「サンプルサイズ」です。そのため、検定力と聞くと「とにかくサンプルの数を増せばいいんでしょ？」と思われる方もいるかもしれません。
まあそれはそうなのです。が、ただ、そうは言っても「やみくもにサンプルの数を増せばいい」というわけでもありません。
そもそも、検定における「検定力」「サンプルサイズ」「効果量」「有意水準」の中で、私たちが知りたい「調査対象において実際に起きていること」を最も直接的に反映しているのは何といっても「効果量」なのです。その一方、「有意水準」「検定力」というのは「実験・調査デザイン（のクオリティ）」の方をより直接的に反映する量といえます。
その意味で、特に論文の査読という「実験・調査デザインのクオリティ」を重視するプロセスにおいて「有意水準」や「検定力」が重視されてきたのは理解できますが、そのせいで「効果量」の方が軽視されるようになってしまったら、それはもう本末転倒と言えるでしょう（例えばこちら）。
（殆どの場合において）大事なのは「調査対象において実際に起きていること」を理解することであり、そのためには「効果量」を中心に考えを進めていくことが重要です。
サンプルサイズを決定する際にも、そもそもの本来的な目的に照らして「どの程度の"差"が実質的な差と言えるのか？」という点をまず明確にし、その「実質的な差（実質的な効果量）を検出するためにはどの程度のサンプルサイズが必要なのか？」という順番で考えていくことが王道といえます。（ここでちゃんと考えることをサボると、あとで「nが足りな〜い...」という番町皿屋敷的な状態になるので気をつけましょう）
逆に言うと、サンプルサイズをやみくもに多くして、「瑣末な差」に統計的有意性を見いだしたとしても、それは「現実を知る／現実に対処する」ためには何の役にも立ちません*13。
（ついつい忘れがちですが）統計解析の際にはいつでもその本来的な目的を見失わないように注意しましょう。
では、今回の内容をまとめます。（殆どさっきのまとめと重複しますが）
はい。
ちなみに今回の記事では仮説検定の枠組みを基にして書きましたが、またそもそも論を言うと、「効果量」に着目するのであればそもそも仮説検定の枠組みよりも「効果量の信頼区間」や「効果量の推定分布」を求めたほうが断然スジが良いと言えます。この辺りを書きはじめるととても長くなるので、ぜひ以下の本をご一読いただければ幸いです（すごく良い本だと思うので）：
（もしくは、みんな仮説検定のことなんか忘れてベイジアンになっちゃえばよいのにと思います）
はい。
では、次回の後編では、「因果関係があるのに相関関係が見られないケース」の中でも、交絡要因や中間変量が関連するものについて書いていきたいと思います。
（今後はコンスタントに更新していきます）
*1:counterfactualさんのブコメを見てここに量的な表現がくるのはおかしいかなと思ったので改変
*2:ちなみに、私は因果関係に対しては「ある／なし」、相関に対しては「見られる／見られない」「現れる／現れない」の語を使うことが多いですが、これは因果は本質的に存在論的なものであるのに対し、相関は本質的に現象論的なものである、という理解の仕方が背景にあります
*3:ブコメで「サンプル数」の語は誤りとの指摘あり。まじどもサンクスです。勉強になります。以下の文では瑣末になるのでこっそり修正しました
*4:つまり等分散性を仮定
*5:各群のデータは正規分布すると仮定します
*6:体重の「増加」に興味があるので片側検定を用います
*7:その都度サンプルをランダムに発生させているので、コードを実行するごとに値は変わります
*8:一般には検定力は0.8くらいあることが望ましいと言われていたりします→詳しくは大久保・岡田(2012)
*9:厳密に言えば片側検定なので「B群の体重の方が大きい場合に」
*10:ここはちょっと大雑把な言い方になっています。実際には効果量の定義にはいろいろあります。詳しくは大久保・岡田(2012)などをご参照ください
*11:あるいは、何も考えずに単に慣習に従って
*12:「効果量」の定義はいろんな形があり得ますが、ここではとりあえずこの形で
*13:もちろんそこで得られた「効果量」には意味があります／まあそれがサンプルサイズの増加に費やした努力量に見合ったものかは分かりませんが
どもです。林岳彦です。息子の3DSにバーチャルコンソールの「ソロモンの鍵」を密かに入れました（まだ３面）。
さて。
前回の記事：
因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） - Take a Risk：林岳彦の研究メモ
につきましては沢山ブクマ等をいただき大変ありがとうございました*1。大変感謝しております。
さて。上記記事について、ublftboさんから「相関関係の定義が書かれていないのでは」（相関と因果 - Interdisciplinary）とのご指摘をいただいたきました。
ご指摘は確かにごもっともですので、今回は「相関」概念についてと、そのついでに近年に開発された"21世紀の相関（MIC）"の話について私なりに書いてみたいと思います。
（以下、ややマニア向けの話になるかもしれません。あと前回ほどではないですが、それなりに長いです。）
まずは、前回の記事の補足的なところから話をはじめたいと思います。
前回の記事において「相関」という語を私がどういう意味で用いていたかというと、あまり深くは考えてはいませんでしたが、改めて文字にすると「得られたデータにおける項目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 より引用
ここでそれぞれの図はデータの散布図を表していて、上の数字はそれぞれの相関係数を表しています。
上の図から：
という「相関係数」の特徴がわかると思います。
はい。いわゆる「相関係数」とは、こういうものなんです。
さて。上記の特徴を理解して使えば相関係数は大変便利なものです。しかし、世の中の全てが「直線的関係」だと思うなよ、と言われたらそれはまさしくそうであります。
そこで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では非線形的な関連が捉えられているようです。
正直ちょっと「ほんまかいな」と思わないでもないですが、研究の「とっかかり」を得る分には十分なのかなあとも思います。
まとめます。
今回は：
てな話でした。
次回は、「因果関係があるにもかかわらず相関関係（＊より正確には"統計的関連"）が生じないケース」についてまとめたいと思います。
私も寄稿させていただいた「線と情事」が、ニコニコ超会議2と併せて行われる「超文学フリマ」において以下の要領において販売される予定です。
「超文学フリマ」in ニコニコ超会議2:
2013年04月28日（日）10:00〜17:00幕張メッセ（「ニコニコ超会議2」内）ブース: 【イ-03】
ニコニコ学会βなどのついでにでもお立ち寄りいただければ幸いです。何卒よろしくお願いいたします。
＊通販でも買えます
甘茶茂、しまおまほ、島田虎之介、関根美有ほか「線と情事 創刊号」 - タコシェオンラインショップ
創作・オリジナル 同人誌専門店 サッシ / 「線と情事」 / NECOfan
.
