2022-11-23

統計的検定でよく使われる手法をまとめる

思い出すのに時間が掛かってしまうことがあるので、よく行われる検定においてどの手法がどのように使われるか、なぜ使うことができるかをまとめようと思います。

母比率の検定

母比率の検定においては、 $n$ が十分大きいときに二項分布が正規分布に近似できることを用いて $z$ 検定を行います。詳細な流れは以下のようになります。

まず、成功確率が $p$ のベルヌーイ試行を考え、これを $n$ 回独立に繰り返すとします。このとき、 $n$ 回の中で出てきた成功回数を $X$ という確率変数で表すと、 $X$ は二項分布に従います。すなわち、

$\displaystyle X \sim B(n, p)$

と表されます。また、このとき $n$ が大きいと二項分布は正規分布に近似できるため、 $X$ は平均 $np$ 、分散 $np(1-p)$ の正規分布に従います *1 *2 。すなわち、

$\displaystyle X \sim N(np, np(1-p))$

となります。更にこの $X$ を $n$ で割った確率変数、すなわち標本の比率（これを $\hat{p}$ とします）は正規分布の再生性により次のように表されます *3 。

$\displaystyle \hat{p} \sim N(p, \frac{p(1-p)}{n})$

そしてこの確率変数を標準化することで検定が行われます。

このように、母比率の検定では $n$ が十分大きいときに二項分布が正規分布に近似する性質を利用し $z$ 検定を行います。

補足

少し混乱しやすいですが、母比率の検定では上記のように $n$ が大きいときの近似の性質を使っているため、「得られたデータが正規分布に従っているか」ということを考える必要はありません。すなわち、正規性の検定を行う必要はありません。

母比率の差の検定

母比率の差の検定においては、各群が正規分布に従うことと正規分布の再生性を利用し、 $z$ 検定が使われます。詳細な説明としては以下の通りです。

まず、1群目のサンプルサイズを $n_1$ 、母比率を $p_1$ 、標本比率を $\hat{p_1}$ とし、2群目も同様に $n_2$ 、 $p_2$ 、 $\hat{p_2}$ とします。このとき、 $n_1$ と $n_2$ が十分に大きければ、標本比率 $\hat{p_1}$ と $\hat{p_2}$ はそれぞれ正規分布に従います。

また、2群の標本比率の差 $\hat{p_1} - \hat{p_2}$ も正規分布の再生性の性質により平均 $p_1 - p_2$ 、分散 $\frac{p_1(1-p_1)}{n_1} + \frac{p_2(1-p_2)}{n_2}$ の正規分布に従います。すなわち、以下のように表すことができます。

$\displaystyle \hat{p_1} - \hat{p_2} \sim N(p_1 - p_2, \frac{p_1(1-p_1)}{n_1} + \frac{p_2(1-p_2)}{n_2})$

そして、この確率変数を標準化することで検定が行われます。

補足

母比率の差の検定においても事前に正規性の検定を行う必要はなく、また等分散性の検定も不要です *4 。

母平均の検定

母平均の検定では、母分散が既知で $n$ が大きいときには中心極限定理により標本平均が正規分布に従うため、 $z$ 検定を行います *5 。

一方、母分散が未知または $n$ が小さいときには、母集団が正規分布に従っていると考えられる場合、 $t$ 検定を使います。母集団が正規分布に従っていると考えられるかは正規性の検定によって調べ *6 、仮定ができる場合は $t$ 検定を行う、という流れになります *7 。

母平均の差の検定

ここでは2標本に対応がない場合を考えます。このとき、母平均の差の検定では両標本の母集団が正規分布に従っていると考えられる場合、ウェルチの $t$ 検定を使います。

そのため、流れとしてはまず両標本について正規性の検定を行い、その結果正規分布の仮定を満たせる場合にウェルチの $t$ 検定を行う、となります。

なお、等分散性の検定を行ったのち、スチューデントの $t$ 検定 or ウェルチの $t$ 検定のどちらを使うかを判断する方法もありますが、こちらは検定の多重性の問題が生じたり *8 、等分散性の検定結果に関する問題が生じたりする *9 ため、等分散性の検定を行わずウェルチの $t$ 検定を行うことが推奨されています。

正規性の検定で正規分布が仮定できなかった場合は、ノンパラメトリックな手法であるマン・ホイットニーのU検定などを使います *10 。

*1: $n$ が大きいと言えるための条件は、文献によって書かれていることが異なるようです。例えば、以下の Wikipedia のページでは「期待値 $np$ および $np(1-p)$ が 5 よりも大きいとき」と書かれていますが、もう一つ貼った別の解説記事では「 $np > 5$ かつ $n(1-p) > 5$ のとき」と書かれています。 ja.wikipedia.org techtipshoge.blogspot.com

*2:証明は以下のページに載っています。physnotes.jp

*3:正規分布の再生性の内容・その証明は以下の記事に載っています。mathlandscape.com

*4:等分散性の検定が不要なのは、母比率の差の検定では正規分布の再生性を用いて検定しており、 $t$ 検定のように等分散という仮定を置いていないためです（そのため、その仮定が合っているかという検証もする必要がありません）。実際、以下の記事でも等分散性の検定が行われていません。bellcurve.jp

*5: $n$ が大きいと言えるための条件には厳密な定義はなく、サンプルサイズが26以上または31以上と慣習的に言われています。ただ、 $t$ 分布は極限として正規分布を含むので、 $n$ が大きいときでも $z$ 検定ではなく $t$ 検定を使うという考え方もあります。書籍「Rで学ぶ確率統計学　一変量統計編」の P124 より。

*6:正規性の検定をする際には、ダゴスティーノのK二乗検定などを使います。詳細については以下のページが参考になります。bellcurve.jp docs.scipy.org

*7:正規分布が仮定できなかった場合の検定方法についてはうろ覚えですが次のようなものだったはずです。まず、そもそもの $n$ を大きくし、そうすると

不偏分散 $s^{2}$ が母分散に近づいていく（不偏分散の一致性による）
$\frac{\sigma}{\sqrt{n}}$ と $\frac{s}{\sqrt{n}}$ の差が分母の $\sqrt{n}$ の作用により小さくなっていく（ $n$ が大きくなると分母も大きくなり差が小さくなる、という意味です）

となります。そのため、中心極限定理が適用できる（すなわち標本平均が $N(\mu_0, \frac{s^{2}}{n})$ に従う）と考え、 $z$ 検定を行います。出典元を忘れてしまったため、後日見つけたときに内容の正誤を確認し修正も行います。

*8:詳細は以下のページの通りです。ja.wikipedia.org

*9:等分散性の検定の結果、帰無仮説（等分散である）が棄却されなかった場合、これは「等分散である」ということが支持される訳ではなく、等分散かが分からないという情報が得られたにすぎない、という問題です。書籍「Rで学ぶ確率統計学　一変量統計編」の P136 より。

*10:詳細は以下のページに載っています。ja.wikipedia.org

2022-11-12

データ分析をするときに意識していること

データ分析者の果たすべき役割として良く言われることに、「価値ある意思決定をできるようにサポートすること」と「意思決定を効率化できるようにサポートすること」というものがあるが、このサポートの塩梅が非常に難しいことがある。

例えば、「データ上大きな効率化が期待できそう」と思われる案があったとしても、現場の人にとっては非現実的なこともある。また、分析の結果非常に価値があると思った施策案があったとしても、ビジネス的な観点で「今やるべきではない」と判断されることもある。

前者は実際の実現性を理解できていないがために無茶な効率化になってしまい、後者は事業に対する理解が足りていないために施策案が実際のアクションに結実しない。この二つに共通しているのは、データ分析の結果「価値がある」「効率化できる」と思っていることが、他者の視点から見て同じように思われるわけではない、ということだ。

なので、データ分析をするときには相手の視点も入れながら提案等を行っていく必要がある。

具体的な方法としては、例えば分析に入る前に事前に「根本的な課題は何か」「どこまでの解決ができると良いか」についてヒアリングをして、自分が分析依頼者と同じ視点で物事を見られるようにする、という方法がある。また、提案を作る過程においても、ときおり依頼者からフィードバックをもらったり、案を複数用意してその中から相手がベストなものを選べるようにする、という方法もある。

一方で、上記の話とは別に、分析をしても "意思決定の質を高めるような" or "意思決定を効率化できるような" 案が見つからないこともある。そのような場合は無理して案を出すのではなく、むしろ分析依頼者とともにどうPDCAを回していくかを決めた方が良い。

例えば、手元に十分なデータがない場合などは分析をしたとしても断片的な状況しか分からない事が多い。そのようなときには、データから分かること・分からないことを依頼者に共有し、その上で何をしていくか（まずはデータを収集できる仕組みを作るのか・断片的な情報から施策案を考えていくのか、その後はどのように改善を進めていくのかなど）を決めた方が、分析チームと依頼者のいるチーム両者にとってスムーズな動きができる。

このように、分析を用いた問題の解決がすぐにできそうにないときには、PDCAの進め方から共に決めたり、それがスムーズにできる環境を整えたりするのがデータ分析者の役割になってくる。

上記二点（「依頼者の視点も入れながら提案を行うこと」と「状況に応じてPDCAの進め方から決めること」）が、私がデータ分析をするときに意識していることである。今回はあとで振り返ることができる形で残しておきたいと思い、記事に起こした。

2022-11-03

XGBoostについての備忘録

XGBoost の公式ドキュメントを見るたびに、毎回「どのような式変形でこの形になるのだっけ」となり、思い出すのに時間がかかる箇所があったので、備忘録としてまとめようと思います。

xgboost.readthedocs.io

式変形の箇所

「Tree Boosting」という見出しの中の「Additive Training」の部分です。

Additive Training

... We write the prediction value at step $t$ as $\hat{y}^{(t)}$ . Then we have

$\displaystyle \hat{y}^{(t)} = \sum_{k=1}^{t} f_k(x_i) = \hat{y_i}^{(t-1)} + f_t(x_i)$

It remains to ask: which tree do we want at each step? A natural thing is to add the one that optimizes our objective.

$\displaystyle obj^{(t)} = \sum_{i=1}^{n} l(y_i, \hat{y_i}^{(t)}) + \sum_{i=1}^{t} \omega(f_t) \\ = \sum_{i=1}^{n} l(y_i, \hat{y_i}^{(t-1)} + f_t(x_i)) + \omega(f_t) + constant \qquad (1)$

... So in the general case, we take the Taylor expansion of the loss function up to the second order:

$\displaystyle obj^{(t)} = \sum_{i=1}^{n} [ l(y_i, \hat{y_i}^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i f_t^{2}(x_i) ] + \omega(f_t) + constant \quad (2)$

where the $g_i$ and $h_i$ are defined as

$\displaystyle g_i = \partial_{\hat{y_i}^{(t-1)}} l(y_i, \hat{y_i}^{(t-1)}) \\ h_i = \partial^{2}_{\hat{y_i}^{(t-1)}} l(y_i, \hat{y_i}^{(t-1)})$

この最後のテーラー展開の仕方を思い出すときに時間が掛かってしまいます。

式変形の詳細

では、上記の変形がどのように行われているか詳細を見ていきましょう。

まず、テーラー展開そのものについて復習です。テーラー展開は関数 $f$ を以下のような無限和で表現したものを言うのでした。

$\displaystyle f(x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x-a)^{n} \\\ = f(a) + f^{'}(a)(x-a) + \frac{f^{''}(a)}{2!}(x-a)^{2} + \frac{f^{'''}(a)}{3!}(x-a)^{3} + \cdots$

なお、このとき点 $a$ での情報を用いているため、「 $a$ まわりでテーラー展開する」と言います。

また、この $f(x)$ の $x$ を $x+\Delta$ とし、 $x$ まわりで展開することで、以下のように別の表現をすることができます。

$\displaystyle f(x + \Delta x) \\\ = \sum_{n=0}^{\infty} \frac{f^{(n)}(x)}{n!} (x+\Delta x - x)^{n} \\\ = \sum_{n=0}^{\infty} \frac{f^{(n)}(x)}{n!} (\Delta x)^{n} \\\ = f(x) + f'(x) \Delta x + f''(x) (\Delta x)^{2} + f'''(x) (\Delta x)^{3} + \cdots$

XGBoost のドキュメントに載っている式変形では、こちらの表現が使われています。具体的には、前述の引用文にある (1) 式

$\displaystyle obj^{(t)} = \sum_{i=1}^{n} l(y_i, \hat{y_i}^{(t-1)} + f_t(x_i)) + \omega(f_t) + constant$

に対して、テーラー展開の $f$ ・ $x$ ・ $\Delta x$ を以下に置き換えて展開しています。

$f$ を $l$
$x$ を $\hat{y_i}^{(t-1)}$
$\Delta x$ を $f_t(x)$

すなわち、以下のように関数 $l$ を表しています。

$\displaystyle l(y_i, \hat{y_i}^{(t-1)} + f_t(x_i)) \\\ = \sum_{n=0}^{\infty} \frac{l^{(n)}(y_i, \hat{y_i}^{(t-1)})}{n!} (f_t(x_i))^{n} \\\ = l(y_i, \hat{y_i}^{(t-1)}) + l'(y_i, \hat{y_i}^{(t-1)}) f_t(x_i) + \frac{1}{2} l''(y_i, \hat{y_i}^{(t-1)}) (f_t(x_i))^{2} + \cdots$

なお、

$\displaystyle l^{(n)}(y_i, \hat{y_i}^{(t-1)}) = \partial_{\hat{y_i}^{(t-1)}}^{(n)} l(y_i, \hat{y_i}^{(t-1)})$

です。

そして、テーラー展開後の2次の項までで表される式を $l$ の近似式とし、 $l'(y_i, \hat{y_i}^{(t-1)})$ を $g_i$ 、 $l''(y_i, \hat{y_i}^{(t-1)})$ を $h_i$ と置くと、関数 $l$ は次のように表されます。

$\displaystyle l(y_i, \hat{y_i}^{(t-1)} + f_t(x_i)) \approx l(y_i, \hat{y_i}^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i (f_t(x_i))^{2}$

この近似式を (1) 式に入れると、

$\displaystyle obj^{(t)} = \sum_{i=1}^{n} [ l(y_i, \hat{y_i}^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i (f_t(x_i))^{2} ] + \omega(f_t) + constant \qquad$

となり、(2) 式と同じになります。このようにして展開されています。

補足

上記では (1) 式

$\displaystyle obj^{(t)} = \sum_{i=1}^{n} l(y_i, \hat{y_i}^{(t-1)} + f_t(x_i)) + \omega(f_t) + constant$

に対して、テーラー展開 $f(x + \Delta x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(x)}{n!} (\Delta x)^{n}$ の $f$ ・ $x$ ・ $\Delta x$ を以下に置き換えて展開しました。

$f$ を $l$
$x$ を $\hat{y_i}^{(t-1)}$
$\Delta x$ を $f_t(x)$

しかし、XGBoost のドキュメントを読むとステップ t においては既に $\hat{y_i}^{(t-1)}$ の値が定まっており、定数として扱われています。

そのため、変数ではないのに $x$ の置き換え先として $\hat{y_i}^{(t-1)}$ を使って良いのかと思うかもしれませんが、論理としては、一旦 $\hat{y_i}^{(t-1)}$ を変数として扱いテーラー展開しても、入出力の関係を表す関数としては変わりがないため、このようにしても問題ないのだと思われます。

また、これと同じように、ステップ t において $f_t(x)$ は値が定まっていない変数ですが、一旦定数として見て $\Delta x$ の置き換え先にしたとしても、入出力の関係を表す関数として変わりがないため問題ないのだと思われます。

参考記事

テーラー展開の説明：

computationjp.com

テーラー展開を用いた式変形の大まかな流れ：

qiita.com

2022-10-12

BigQueryで統計的検定をできるようにする

BigQueryで検定をしたい

分析データの格納先として BigQuery（BQ）が多くの会社で使われていますが、実際に分析をするときに統計的検定を行いたいことがあると思います。

しかし、BQ で検定をするためには面倒なことも多く、「もっと手軽にできるようにしたい」と感じることがあります。

そんな場面で役立つようなツールを今回作りました。現時点で対応している検定の種類は少ないですが、このツールを使えば BigQuery 内で簡単に検定を行えるようになります。

BigQueryでの検定はなぜ面倒か？

そもそも BigQuery 内での検定がなぜ面倒かというと、検定が BQ 内部で完結せず、何かしら外部のツールを使わなければいけないからです。

BQ ではデフォルトで備わっている統計関連の処理として、算術平均を計算する AVG 関数や不偏標本分散を計算する STDDEV_SAMP などの関数があります。ただ、残念ながら現状これらの結果を用いて検定を行う機能はありません。

また、BQ には User Defined Functions（UDF）という機能もあり、有志の方が統計処理を行うための関数を作っていらっしゃいますが、例えばt検定を行う t_test 関数では算出されるのがt値のみでp値は自分で確認しなければなりません（以下の GitHub ページから関数の詳細を確認できます）*1。

github.com

そのため、現状 BQ で統計的検定をしようとすると、おおよそ以下のような方法で行うことになります。

前述のUDFを実行し、t分布表などを元にp値を出す。それを元に検定結果を出す。
BQ のデータを Python で扱い、scipy や statsmodels パッケージを使って検定する。

ここで問題なのが、どちらの方法においても何かしら外部のツールが必要になることです（一点目の方法なら分布表、二点目の方法なら Jupyter Notebook などの Python の実行環境）。

検定をする毎にこのように何かしら使わなければならないのは非常に骨が折れます。

どのような解決を行ったか？

そこで今回利用したのが、現在プレビューで提供されているリモート関数という機能です。こちらを使うことで検定を BQ 内で完結させ、前述の手間を無くすことができます。

cloud.google.com

この機能について簡単に説明すると、リモート関数とは BigQuery から Cloud Functions / Cloud Run を実行することができる機能で、これにより BQ のデータに対して何らか加工を行い、その結果を BQ に返すことができるようになっています。

今回私が作成したツールでは、この機能を使って検定に必要なデータを Cloud Functions に送り、それを元に検定を行い結果を返す、という処理を行うようにしました。図示すると以下のようになります。

BQのクエリ内で分析データから検定に必要な情報を取得し、Cloud Functionsで検定を行う

実装の詳細

以下は検定をする上では必須の情報ではないため、必要に応じて「## 実際に検定をしてみる」まで読み飛ばすなどして頂ければと思います。

今回実装したリモート関数の詳細ですが、次のようになっています。なお、ソースコード・Cloud Functions へのデプロイの方法に関しては次の GitHub のレポジトリに載せています。

github.com

bq-remote-functions-stats/
├── .env.template
├── .gcloudignore
├── .gitignore
├── main.py
├── Makefile
├── README.md
├── requirements.txt
├── statstest.py
└── statstests/
    ├── __init__.py
    └── prop.py

この中で主要なファイルの役割は以下のようになっています。

main.py
- BQ から受け取ったデータをパースし、実施したい検定の種類（母比率の差の検定/母平均の差の検定など）と、検定に必要な標本に関する情報を取得する。
- これらの情報を元に statstest.py で定義されたクラスをインスタンス化し、検定を行い、その結果を BQ に返す。
statstest.py
- 検定を行うためのクラス（StaticalTest クラス）を定義したファイル。
- このクラスでは、行いたい検定の種類に応じて statstests ディレクトリ配下で定義されたモジュールを呼び出し、検定を行う。
statstests ディレクトリ
- 各種検定を定義したディレクトリ。現在は母比率の差の検定（prop.py）にのみ対応（今後対応できる検定を増やして行こうと考えています）。
statstests/prop.py
- 母比率の差の検定を定義したモジュール。現在は両側検定にのみ対応。
- 検定自体は statsmodels パッケージで行い、p値などを出す。

実際に検定をしてみる

では、上記のツールを実際に使って BQ で検定を行ってみましょう。

例えば、ある施策を AB テストで実施し、A を従来のパターン、B を新パターンとして CVR に対して有意水準10%の両側検定をしたいとします。すなわち、帰無仮説と対立仮説はそれぞれ以下になります。

帰無仮説： AパターンとBパターンの母比率は等しい
対立仮説： AパターンとBパターンの母比率は等しくない（差がある）

そして、データの集計が終わり以下のような結果が出たとします。

Aパターン：対象者100人、20人がCV
Bパターン：対象者100人、10人がCV

この結果を元に BQ で以下のようなクエリを実行します。

WITH
-- 施策結果のデータを作成
prep AS (
  SELECT
    10 AS successes,   -- AパターンのCV人数
    100 AS trials  -- Aパターンの対象者数
  UNION ALL
  SELECT
    20,  -- BパターンのCV人数
    100  -- Bパターンの対象者数
)

SELECT
  `analysis-test-337607.stats.prop_test`(successes, trials)  -- リモート関数の実行
FROM
  prep

すると、以下のような結果になります。

画像より Pval: 0.047670 となっておりp値が0.1よりも小さいことから、帰無仮説は棄却され、AパターンとBパターンは等しくないと考えられます。このようにして検定結果を導くことができました。

まとめ

今回は BigQuery で検定を行うのに手間が掛かるという問題に対して、リモート関数を使って BQ 内で検定を完結できるようにしました。

まだ対応しているものは少ないですが、Python で出来る検定はリモート関数でも出来るはずなので、これから増やしていこうと考えています。

今回の記事がもし興味深いと思われましたら、スターを付けて頂けると嬉しいです。（対応可能な検定を増やすモチベーションになり、また、次回以降の記事を書く励みにもなります！）

*1:JavaScript で定義した UDF を使えばp値も一度に出すことが出来ますが、両側/片側検定、対応のある/なし、t検定以外の検定をしたくなったときに、BQ 内が非常に複雑になってしまうというデメリットがあります。