なぜリスク分析のプロは仮説検定を使わないのか（ややマニア向け）

お久しぶりです。林岳彦です。もうすぐ『愛なき世界』の日、いわゆる（マイブラッディ）バレンタインデーですね。何かと雑音が多いこの世界ですが、いつでも自分の足元を見つめて行きましょう。

さて。

今回は、以下の：

• そもそもビジネスの現場ではどういう「レベル」の統計学を使うべきなのか - 銀座で働くデータサイエンティストのブログ
• 統計学的検定に対するある拒絶反応: ニュースの社会科学的な裏側
• A/Bテストのガイドライン：仮説検定はいらない（Request for Comments｜ご意見求む） - 廿TT

のあたりの皆様の良記事に触発されて「仮説検定」について何か書いてみようと思いました。で、書こうと思えば色々な側面から書ける気もするのですが、今回はちょっと斜めからのアプローチとして、「リスク分析の人の頭のなかで仮説検定はこんな感じに見えている」というところを書いていきたいと思います。

ここで、ひとくちに「リスク分析」と言っても思い浮かべるところは人それぞれかもしれませんが、今回は以下の：

入門リスク分析―基礎から実践

に載っているような「確率論的リスク分析手法」を用いるアプローチ一般を指して「リスク分析」という言葉を使っていきます*1。

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

（今回もとても長いです。いつもすみません。。。あと、今回とくに一部内容に余り自信がないので*2、もし間違ってたら適宜ツッコミおねがいいたします）

前置き：リスク分析では仮説検定は殆ど使わなかったりします

さて。

もしかしたら意外に思う方もいるのかもしれませんが、「リスク分析」においては仮説検定を使うことは殆どありません。私も一応リスク分析に類する論文を何報か書いていますが（論文へのリンク）、それらの論文の中で仮説検定を用いたことは一度もありません*3。

これは単に私が個人的に仮説検定を使わないという話ではなく、例えば、冒頭に挙げたDavid Vose著の『入門リスク分析』の中でも仮説検定の話はほんの僅かしか出てきません*4。これは、そもそも「リスク分析」のツールボックスの中に「仮説検定」が基本的に含まれていないことを意味しています。

では、なぜリスク分析においては仮説検定が殆ど使われないのでしょうか？

なぜリスク分析においては仮説検定は使われないのか：本来的な理由と技術的な理由

リスク分析において仮説検定が殆ど使われない理由としては、本来的なものと技術的なものの２通りの理由があります。

まずは、本来的な理由の方から説明してみましょう。

本来的に、「リスク分析」とは、リスク（＝影響と不確実性の在りよう）を推定し、意思決定者の参考となる情報を提供することが基本的な目的です。つまり、「意思決定支援」こそがリスク分析の本懐であるわけです。

ここで、「仮説検定」がリスク分析と本来的に相性が悪いところは、仮説検定が「意思決定の自動代行機能」の側面を持っているところです。つまり、リスク分析の目的が「意思決定支援」であるにもかかわらず、仮説検定は仮説の採択に関して論点を先取りして肝心の「意思決定」をオートマティックに代行してしまう*5側面があり、この性質は「リスク分析のツール」としてはとても悪質なものであると言えます（→参考：具体的な弊害についての過去記事）。

もう一つの、技術的な理由の方は、リスク分析では推定の不確実性を見るときに基本的に「不確実性の分布全体を見る」ので、仮説検定がそもそも不要であることです。

さて。

本記事では以下、この「不確実性の分布全体を見る」ということについて、例を交えつつ説明していきたいと思います。

「不確実性の分布全体を見る」とはどういうことか：ブリの新エサを仮想例として

では、「不確実性の分布全体を見る」とはどういうことなのかを見て行きます。

説明のための仮想例として、「ブリの養殖業」における、ブリをより大きく育てるための「新エサの開発」についてのケースを考えてみましょう。

天然ブリ　1本-6kg前後

あなたはブリの養殖業を営んでおり、ブリを大きく育てるための新エサを開発しているとします。そして実際に、その新エサの効果を調べるために、新エサを用いて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。

ブートストラップ法というのは：

• (1) 現在得られているサンプル（今回の例の場合、30匹のブリの体長データ）から再び無作為に重複を許しながら値をサンプリングすることにより、新たな「サンプル」（今回の場合、新たな「30匹分のデータ」）を生成し、その「新たなサンプル」に基づき統計量（今回の場合、平均値）の推定値を計算する
• (2) (1)の計算をものすごくたくさん繰り返すことにより、ものすごくたくさんの数の統計量の推定値を算出し、その推定値の分布全体を得る

という方法です。（なんのこっちゃかよく分からない方は、基礎統計学講座＠ウィキã‚„フリーソフトによるデータ解析・マイニングã‚„奥村晴彦さんの記事のあたりの良記事をオススメいたします）

まあ案ずるより産むが易しということで、とりあえずやってみます*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)

という良書を拝読させていただき大変勉強になったのですが、この中に今回示したような「リスク分析」、つまり「データ解析」の結果を二次利用してさらなる解析を行い、意思決定により密接に関連するコンテクストにおける指標への影響を提示しうるような、そんな解析アプローチも載っていても良いのではなかろうか、という感想を抱きました。

スーパーで出来合いの「切り身」を買ってくるようなデータサイエンティストではなく、クライアントのニーズに沿って魚を一本まるごといかようにでも捌ける「寿司屋の大将」のようなデータサイエンティストを目指すのならば、そんな「リスク分析という包丁」も一つ持ってると身を助くこともあるのかも、とか思ったりする今日このごろです（←無責任）。

で、もし、本記事がそのような「寿司屋の大将」への入り口となりうればとても嬉しいと思っております。

# 「リスク分析」をもっと知りたい方にはこちらをオススメいたします↓

入門リスク分析―基礎から実践

（ちなみにわたくし、この本は「ベイズ統計入門」としてもベストと言っても良いかもしれないと密かに思っております）

Appendix：今回参考にさせていただいた記事のリスト（多謝）

仮説検定関連：

• そもそもビジネスの現場ではどういう「レベル」の統計学を使うべきなのか - 銀座で働くデータサイエンティストのブログ
• 統計学的検定に対するある拒絶反応: ニュースの社会科学的な裏側
• A/Bテストのガイドライン：仮説検定はいらない（Request for Comments｜ご意見求む） - 廿TT

ブートストラップ関連：

• 基礎統計学講座 @ ウィキ - ブートストラップ法
• Ｒとブートストラップ
• ブートストラップ
• 基礎統計学講座 @ ウィキ - ブートストラップ法

信頼区間関連：

• ブートストラップ法で信頼区間を求めるときの注意点 - ほくそ笑む
• ［統計］有意な差があるとは信頼区間が重ならないこと？: 馬車馬のように
• 研究者の多くはエラーバーの意味をろくに理解していない - Kamedoの音風景

Rによるベイズ解析（MCMC）の実装関連：

• Stanで統計モデリングを学ぶ(1): まずはStanの使い方のおさらいから - 銀座で働くデータサイエンティストのブログ
• Stanで統計モデリングを学ぶ(2): そもそもMCMCって何だったっけ？ - 銀座で働くデータサイエンティストのブログ

あんまり関係ないかもですがリスク分析における「変動性」と「不確実性」について自分で昔に書いた記事：

• http://takehiko-i-hayashi.hatenablog.com/entry/20100701/1277954027

*1:なので、もし読んでいて違和感がある場合には、「リスク分析」を「確率論的リスク分析」に適宜置き換えて読んでいただければと思います

*2:特に古典的な「統計学」に関するところが苦手です

*3:文献からの既往データとして仮説検定からの結果を（仕方なく）使うことはあります

*4:本書の全567ページ中で仮説検定の話が出てくるのは僅か10ページほどで、しかも内容は分布モデルの選択の文脈での適合度検定の話なので、本来なら仮説検定よりもAICを使うべき部分になります

*5:あるいは、代行しうると多くの人にみなされている

*6:ベイズ法についての例はたとえば→ tjoさんの記事

*7:まあ30匹程度でブートストラップやるのは精度が低いかもですが今回はその辺はあまり気にせずにやっていきます

*8:実装のコードはいろいろ方法はありますが、今回はsample関数を使いました。分布の形状が綺麗になるように多めに50000回の反復をさせています

*9:すみません。この式はテキトーに考えたもので特に意味はありません

*10:説明の単純化のため、ブリの体長の分散の影響と、この関係式自体に関する不確実性は無視できるものと仮定します

*11:今回の例は説明のためかなり単純化していますが、実際には多数のコンテクストからの多数のパラメータが絡み合う複雑な解析になることももちろんあります

*12:xパーセンタイル = 小さい方から数えて全体のx%目の順番に相当する値

*13:あんまりこのへん自信ない