p値と仲良くなる - エムスリーテックブログ

＜エムスリー Advent Calendar 2020 まで残り5日となりました。Advent Calendar本編に先んじて新卒1〜2年目メンバーが執筆します。＞

エムスリーエンジニアリンググループ AI・機械学習チームの金山 (@tkanayama_)です。

エムスリーでは、主に施策の継続判断などを行う際に、統計的仮説検定などを利用してデータに基づいた総合的な判断をしています。エムスリーの場合はエンジニアリンググループとは別に、データ分析を専門に扱うデータ分析グループが存在するので、エンジニア全員が必ずしも仮説検定の数理的な背景まで理解している必要はないかもしれません。しかし個人的には、意思決定のツールとして利用している以上、誤用しないだけの最低限の素養はエンジニアにもあると望ましいと思います。

今回は、データに基づいた意思決定を行う中でよく使われる用語の1つである「p値」に焦点を当てたいと思います。

f:id:tepppei:20201125154719j:plain:w600

用語の簡単な説明

統計にあまり馴染みがない方向けに説明すると、統計的仮説検定では、否定したい仮説（例えば、「この施策は効果がない」など）を「帰無仮説」、示したい仮説（例えば、「この施策には効果がある」など）を「対立仮説」と呼びます。そして、p値は「帰無仮説を仮定した時に、今手元にあるデータよりも極端な値が得られる確率」と定義されます。もしp値が十分小さければ、「帰無仮説が正しいと仮定すると、手元にあるデータが得られる確率はとても低い→帰無仮説は間違っているだろう」というロジックで帰無仮説を否定できます。

p値は意思決定のための定量的で明確な基準を与えてくれるためとても便利ですが、盲目的に数字だけを追いかけると誤った意思決定に繋がる可能性があります。p値や統計的仮説検定に関する理論上重要な話は素晴らしい書籍*1がたくさん存在するのでそちらを参照していただくとして、今回は少し変わった観点からp値を眺めてみます。直接実務に役立つことはないと思いますが、p値をより多面的に知ることで、誤用や思い込みによる誤った意思決定を防ぐことに繋がれば幸いです。

この記事でやりたいこと

統計的仮説検定において、検定対象の真のパラメータが

帰無仮説で仮定したパラメータと大きく乖離している場合
帰無仮説で仮定したパラメータとあまり乖離していない場合

それぞれにおいて、p値の実現値がどのような形の分布従うのかを計算してみます。

準備

以下では、広告や記事などの推薦システムを想定して、あるアイテムに対して発生するクリックモデル化する場合を考えます。必要に応じて、関心のあるドメインに読み替えてください。

平均 $p_{actual}$ のベルヌーイ分布に従う確率変数 $X$ を導入し、あるアイテムがクリックされる事象を $X=1$ 、クリックされない事象を $X=0$ としてモデル化します。このとき、アイテムの真のCTRは $p_{actual}$ であり、独立な $n$ 個の観測値から推定されるCTRは $\bar{X} = \frac{1}{n} \Sigma_{i=1}^n X_{i}$ と表現できます。

ここで、nが十分大きい場合*2、 $\bar{X}$ は以下のような正規分布に従うと近似できます。

$\bar{X} \sim N(p_{actual}, \frac{1}{n} p_{actual}(1-p_{actual}))$

仮説検定

ここで、定数 $p_{hypo}$ を指定して、

帰無仮説： $p_{actual} = p_{hypo}$
対立仮説： $p_{actual} > p_{hypo}$

という仮説検定を行います*3。

以下では、 $p_{actual}$ と $p_{hypo}$ にそれぞれ具体的な値を設定して、 $\bar{X}$ の実現値から計算されるp値がどのような形の分布に従うか計算してみたいと思います。

仮説との乖離が小さい場合

まず、 $p_{hypo}$ と $p_{actual}$ の乖離が小さい例として、下記のような状況を考えます。

$p_{hypo} = 0.010$

$p_{actual} = 0.011$

$n=10000$

このとき、「 $X$ が帰無仮説に従う場合」と「 $X$ が真の分布に従う場合」それぞれについて、 $\bar{X}$ の分布を描くと下記のようになります。

f:id:tepppei:20201124121948p:plain — 正規分布に近似した $\bar{X}$ の確率密度分布

例として、 $p値 \leq 0.025$ となる確率を図中に記入しました。ここで、「 $p値 \leq 0.025$ となる $\bar{X}$ の範囲は『帰無仮説が成立する場合の分布』により決まり、 $\bar{X}$ がその範囲に収まる確率は『真の分布』によって決まる」という点に注意してください。

これを踏まえて、p値の分布を下記のようなコードで描いてみました。

import scipy.stats
import numpy as np
import matplotlib.pyplot as plt

# 実験条件
n = 10000
p_hypo = 0.01
p_actual = 0.011
s_hypo = np.sqrt(p_hypo * (1 - p_hypo) / n)
s_actual = np.sqrt(p_actual * (1 - p_actual) / n)

# p値の確率密度関数を計算する
def p_value_pdf(p_value, p_hypo, s_hypo, p_actual, s_actual):
    x_bar = scipy.stats.norm.ppf(1 - p_value, loc=p_hypo, scale=s_hypo)
    return scipy.stats.norm.pdf(x_bar, loc=p_actual, scale=s_actual) / scipy.stats.norm.pdf(x_bar, loc=p_hypo, scale=s_hypo)

# 描画
delta = 0.0001
p_values = np.arange(delta, 1, delta)
p_value_densities = np.array([p_value_pdf(p_value,  p_hypo, s_hypo, p_actual, s_actual) for p_value in p_values])
plt.figure(figsize=(10, 5), dpi=100)
plt.xlabel('p-value')
plt.ylim(0, 40)
plt.xlim(0, 1)
plt.plot(p_values, p_value_densities, linewidth = 3.0)

f:id:tepppei:20201124113906p:plain — p値の確率密度関数

このようなグラフになりました。次に、具体的な積分値を下記のような関数で計算しました*4。

# x_minとx_maxで挟まれた領域の面積を計算する
def calcualte_probability(x_min, x_max, p_hypo, s_hypo, p_actual, s_actual, delta=0.0001):
    p_values = np.arange(x_min + delta, x_max, delta)
    p_value_densities = np.array([p_value_pdf(p_value,  p_hypo, s_hypo, p_actual, s_actual) for p_value in p_values])
    return sum(p_value_densities) * delta

これを用いて、例えばp値が0.01以下になる確率は10%程度であることが計算できます。

仮説との乖離が大きい場合

次に、 $p_{hypo}$ と $p_{actual}$ の乖離が大きい例として、下記のような状況を考えます。

$p_{hypo} = 0.010$

$p_{actual} = 0.015$

$n=10000$

このとき、先ほどと同様に「 $X$ が帰無仮説に従う場合」と「 $X$ が真の分布に従う場合」それぞれについて、 $\bar{X}$ の分布を描くと下記のようになります。

f:id:tepppei:20201124114524p:plain — 正規分布に近似した $\bar{X}$ の確率密度分布

さらに、p値の分布を図示すると下記のようになります。

f:id:tepppei:20201124113941p:plain — p値の確率密度関数

そして、先ほどと同様にp値が0.01以下になる確率を計算すると99%となります。仮説との乖離が小さい場合と比べると、p値の確率密度がほとんど0付近に集中していることがわかります。

まとめと教訓

今回は、真のパラメータが帰無仮説と近い場合と遠い場合の2パターンに対して、p値の分布がどのように変わるかを具体的に見てみました。ここからどのような教訓が得られるでしょうか。

例えば、p値が0.05付近でギリギリ帰無仮説を棄却できなかった時、「偶然この値になってしまったのかもしれないので、もう少しサンプルを増やせば結果が変わるかも…」などと粘りたくなるのが人情です（人は俗にこれをp-hackingと呼びます）*5。しかし、上記のグラフを見ると、十分効果が期待できる施策において、"偶然"p値が0.05付近になってしまうことはほぼありえないということがわかります*6。そのため、無駄に粘るのではなく、「もし効果があったとしてもごくわずかだと予想されるので、潔くこの施策はあきらめて次に行こう」という発想の方が合理的かもしれません（もちろん、改善幅がどれだけビジネスインパクトを持つかどうかにも左右されるので、一概には言えませんが）。

参考文献

ベルヌーイ分布を仮説検定の形に落とし込む際に、「サンプルサイズの決め方（永田靖著）」の3章「1つの母平均の検定」および11章「その他の手法」を参考にしました。「仮説の設定→有意水準の設定→検出力の設定→サンプルサイズの計算」の流れがわかりやすく整理されています。

www.asakura.co.jp

$\bar{X}$ の分布をp値の分布に変換する際に、「プログラミングのための確率統計（平岡和幸著・堀玄著）」の4章「連続値の確率分布」を参考にしました。解析学の教科書のような厳密さを求めるのではなく、確率に対して直感的イメージを持ちたい場合はおすすめの書籍です。

www.ohmsha.co.jp

We are hiring!

エムスリーは、ソフトウェアエンジニアやビジネスサイドの方々であっても、データを根拠にしながら客観的な意思決定をしようという風土が根付いているように感じます。

特に私のように若い社員にとって、このようにデータを根拠に物事を進めていく風土はとても居心地が良いです。なぜなら、「まだ若いから」「まだ経験が浅いから」などの理由で意見が通らないことが起きにくいからです。一緒にデータに基づいた意思決定で事業を推進していきましょう！

エンジニアリンググループ募集一覧

jobs.m3.com

*1:参考文献の章で何冊か紹介しています。

*2:参考文献「サンプルサイズの決め方」には、「 $np \geq 5$ かつ $n(1-p) \geq 5$ のとき、 $n$ が十分大きいとみなす」という基準が記されています。

*3:今回は簡単のため $p_{hypo}$ を外から指定する問題設定にしていますが、実務上はコントロール群・テスト群のパラメータをそれぞれ $p_{control}, p_{test}$ として、真の値はどちらも未知であるという問題設定（2つの母平均の検定）が多いと思います。

*4:端の扱いが雑なので、p値が0付近の確率を計算する場合は、先に反対側の領域を計算して1から引くのが良いと思います。

*5:理想的には、必要な検出力を事前に定義した上でサンプルサイズを決めればこのようなことは起きないのですが…

*6:厳密に議論するなら、p値がgivenの場合の真のパラメータの条件付き分布を見る必要がありそうです。