はじめに

AI研究開発室のS.Y.です。

タイトルの通り今回は、デルタ法を用いた変化率の区間推定を紹介します。
特に、ABテストを行った際のコントロール群 vs 処置群におけるメトリクスの変化率について扱います。
正しく区間推定することで、介入によってメトリクスが何%変動するかについて、正しい期待幅をもって意思決定を行えます。

また本稿では上記の紹介に加えて、合成データでの実験でデルタ法による推定の正しさを確認し、やってしまいがちないくつかのアンチパターンだと正しい推定ができないことも確認しようと思います。

結論は以下です。

• • 比率の分散は、デルタ法によって正しく推定できる。
  • ABテストにおける変化率の区間推定は、\(CI = (\frac{\bar{Y_t}}{\bar{Y_c}}-1) \pm z_{\alpha/2}*\sqrt{\frac{\frac{1}{\bar{Y_c}^2}*var(Y_t)+\frac{\bar{Y_t}^2}{\bar{Y_c}^4}*var(Y_c)}{n}}\)で計算する。
  • 変化率の区間推定になっていなかったり、推定の過程で分散を正しく推定できていなかったりすると、区間を過大・過小に推定してしまう。以下はやりがちなので注意。
    • アンチパターン① (変化率の区間推定になっていない): \(CI = \frac{lower(\bar{Y_t})}{\bar{Y_c}} \sim \frac{upper(\bar{Y_t})}{\bar{Y_c}}\)としてしまう。
    • アンチパターン② (変化率の分散推定が正しくない): \(var(\Delta\%) = \frac{var(\Delta)}{\bar{Y_c}^2}\)としてしまう。

デルタ法を用いた変化率の区間推定

前提

ABテストを実施し、とあるメトリクスについてのコントロール群\(Y_c\)と処置群\(Y_t\)間の変化について区間推定する場合を考えます。

一般的に、信頼度\(\alpha\)の区間推定は点推定値\(pest\)と分散推定値\(vest\)を用いて、以下の式で計算します。
$$CI = pest\pm z_{\alpha/2}*\sqrt{\frac{vest}{n}}$$

また、「メトリクスの平均値」、「メトリクスの標本分散」、「メトリクスの平均値の分散」の計算方法を確認しておきます。
・メトリクスの平均値の計算式: $$\bar{Y} = \frac{1}{n}\sum_{i=1}^{n}Y_{i}$$
・標本分散の計算式: $$var(Y) = \hat{\sigma}^{2} = \frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\bar{Y})$$
・メトリクスの平均値の分散の計算式: $$var(\bar{Y}) = var(\frac{1}{n}\sum_{i=1}^{n}Y_{i}) = \frac{1}{n^{2}}*n*var(Y) = \frac{\hat{\sigma}^{2}}{n}$$

変化量\(\Delta\)の区間推定

では、実際に区間推定について考えていきます。

変化量\(\Delta\)については、\(pest\), \(vest\)はそれぞれ下記のように計算できます。
$$pest = \bar{Y_t}-\bar{Y_c}$$
$$vest = var(\bar{Y_t}-\bar{Y_c}) = var(\bar{Y_t}) + var(\bar{Y_c}) = \frac{var(Y_t)}{n} + \frac{var(Y_c)}{n}$$
したがって変化量\(\Delta\)の推定区間\(CI\)は以下になります。
$$CI = (\bar{Y_t}-\bar{Y_c}) \pm z_{\alpha/2}*\frac{\sqrt{var(Y_t)+var(Y_c)}}{n}$$

変化率\(\Delta\%\)の区間推定

さて、ここからが本題です。変化率\(\Delta\%=\frac{\Delta}{\bar{Y_c}}\)の区間推定はどうやって計算しましょう？
やってしまいがちな(私が過去にやっていた)アンチパターンとして下記のようなものが挙げられるかと思います。
$$CI_{error} = \frac{lower(\bar{Y_t})}{\bar{Y_c}} \sim \frac{upper(\bar{Y_t})}{\bar{Y_c}}$$
ここでやっていることは
$$CI_{error} = \frac{\bar{Y_t}\pm z_{\alpha/2}*\sqrt{\frac{var(Y_t)}{n^2}}}{\bar{Y_c}}$$
で、これは\(Y_t\)についての点推定と分散推定を計算していることになります。

実験では処置群\(Y_t\)の結果により着目するため、\(Y_t\)の区間推定結果を中心に考えてしまう気持ちもわかりますが(わかってほしい)、変化率の区間推定としてはこれでは誤りです。

正しく区間推定するためには、変化率\(\Delta\%=\frac{\Delta}{\bar{Y_c}}\)についての点推定と分散推定をそれぞれ行う必要があります。そのやり方を見ていきましょう。

\(\Delta\%\)の点推定はシンプルで、以下のとおりです。
$$pest = \frac{\bar{Y_t}-\bar{Y_c}}{\bar{Y_c}} = \frac{\bar{Y_t}}{\bar{Y_c}}-1$$

そして分散推定の計算式は以下のようになります。
$$vest = var(\Delta\%) = var(\frac{\bar{Y_t}-\bar{Y_c}}{\bar{Y_c}}) = var(\frac{\bar{Y_t}}{\bar{Y_c}})$$
このような比率についての分散推定は直接的な計算は難しく、近似的な手法を使って求める必要があることが知られています。

ちなみに、こちらの書籍[1]では上記の分散推定のよくあるアンチパターンとして、\(var(\Delta\%) = \frac{var(\Delta)}{\bar{Y_c}^2}\)としてしまう、というものが挙げられています。\(\bar{Y_c}\)自体も観測された値であり、母集団の真のパラメータではないため、これでは正しく推定することはできません。

話を戻して比率の分散推定についてですが、今回はタイトルにもあるデルタ法を用いて行います。
デルタ法はデータがある程度(n>1000)あれば比較的簡単に推定分散を計算でき、データが集まりやすいwebサービスでの実験などのケースにマッチしていると思います。

デルタ法による比率の分散推定は以下のようになります[2]。

$$var(\Delta\%) = \frac{1}{\bar{Y_c}^2} * var(Y_t) + \frac{\bar{Y_t}^2}{\bar{Y_c}^4} * var(Y_c) – 2*\frac{\bar{Y_t}}{\bar{Y_c}^3}*conv(\bar{X},\bar{Y})$$

コントロール群と処置群は独立なので共分散の項は0となり、上記を踏まえて最終的に、Δ%の区間推定は以下のような計算になります。
$$CI = (\frac{\bar{Y_t}}{\bar{Y_c}}-1) \pm z_{\alpha/2}*\sqrt{\frac{\frac{1}{\bar{Y_c}^2}*var(Y_t)+\frac{\bar{Y_t}^2}{\bar{Y_c}^4}*var(Y_c)}{n}}$$

コントロール群、処置群それぞれの標本平均と標本分散のみを使い、変化率の区間推定が行えました。

推定の正しさの確認

下記のそれぞれのパターンで変化率の区間推定を行った場合について、推定の正しさを確認してみます。

• • デルタ法
  • アンチパターン①: 分散推定が誤っている
    • \(var(\Delta\%) = \frac{var(\Delta)}{\bar{Y_c}^2}\)
  • アンチパターン②: 変化率の区間推定になっていない
    • \(CI = \frac{lower(\bar{Y_t})}{\bar{Y_c}} \sim \frac{upper(\bar{Y_t})}{\bar{Y_c}}\)

各推定法の実装

デルタ法

import numpy as np
import scipy.stats as stats

def delta_method_ci(data_c, data_t, alpha=0.05):
    mean_c = data_c.mean()
    mean_t = data_t.mean()
    s_c = data_c.std()
    s_t = data_t.std()
    n_c = len(data_c)
    n_t = len(data_t)
    s_ct = 0
    
    z_alpha_2 = stats.norm.ppf(1 - alpha / 2)  # z_alpha/2の値

    # 各項を分けて定義
    var_t_term = s_t**2 / (n_t * mean_c**2)  # 処置群の分散に基づく項
    cov_term = 2 * (mean_t / (n_c * mean_c**3)) * s_ct  # 共分散に基づく項
    var_c_term = (mean_t**2 / (n_c * mean_c**4)) * s_c**2  # コントロール群の分散に基づく項

    # 分散推定の最終式
    variance_estimate = var_t_term + var_c_term - cov_term

    # 信頼区間の計算
    point_estimate = (mean_t / mean_c) - 1  # 点推定値
    ci_lower = point_estimate - z_alpha_2 * np.sqrt(variance_estimate)  # 信頼区間下限
    ci_upper = point_estimate + z_alpha_2 * np.sqrt(variance_estimate)  # 信頼区間上限

    return point_estimate, ci_lower, ci_upper

アンチパターン①

import numpy as np
import scipy.stats as stats

def error_variance_estimate_ci(data_c, data_t, alpha=0.05):
  mean_c = data_c.mean()
  mean_t = data_t.mean()
  var_c = data_c.var()
  var_t = data_t.var()

  z_alpha_2 = stats.norm.ppf(1 - alpha / 2)  # z_alpha/2の値

  # 信頼区間の計算
  point_estimate = (mean_t / mean_c) - 1  # 点推定値
  error_variance_estimate = (var_c + var_t) / (mean_c**2) # 分散推定値(誤り)
  ci_lower = point_estimate - z_alpha_2 * np.sqrt(error_variance_estimate)  # 信頼区間下限
  ci_upper = point_estimate + z_alpha_2 * np.sqrt(error_variance_estimate)  # 信頼区間上限

  return point_estimate, ci_lower, ci_upper

アンチパターン②

import numpy as np
from scipy.stats import t

def naive_ratio_ci(data_c, data_t, alpha=0.05):
  n_t = len(data_t)
  mean_t = data_t.mean()
  std_t = sem(data_t)
  treat_intervals = t.interval(
      confidence=1 - alpha,
      df=n_t - 1,
      loc=mean_t,
      scale=std_t,
  )
  treat_ci = np.max(treat_intervals) - np.mean(treat_intervals)
  mean_c = data_c.mean()
  point_estimate = (mean_t / mean_c) - 1
  ci_lower = (mean_t - treat_ci)/mean_c - 1 # 信頼区間下限
  ci_upper = (mean_t + treat_ci)/mean_c - 1 # 信頼区間上限

  return point_estimate, ci_lower, ci_upper

 

実験

推定の正しさの指標として、\(coverage = \frac{推定区間が真の変化率を含む回数m}{試行回数n}\)を用います。

区間推定の際に設定する信頼率\(1-\alpha\)と実際の\(coverage\)が等しいほど、その推定法は正しい信頼率で区間推定ができていると言えます。

これ以降いくつかの分布のケースについて、試行回数\(n=10000\)、信頼率95%(=\(\alpha=0.05\))でそれぞれの区間推定法を試した結果の\(coverage\)を見てみます。
サンプル数は、コントロール群・処置群ともに10000件のデータを生成します。

ポワソン分布における変化率の区間推定

コントロール群: lam=50, 処置群: lam=52として、ポワソン分布からそれぞれ10000件のデータを生成します。サイト上のあるボタンについて、ユーザーあたりのクリック数がコントロール群で50回、処置群で52回、といったようなケースです。

推定区間が真の変化率(\(\frac{52}{50}-1=0.04\))を含んでいる割合は以下のようになりました。

• デルタ法: 0.9468
• アンチパターン①: 1.0
• アンチパターン②: 0.8289

デルタ方は0.95に近いcoverageで推定できています。

一方でアンチパターン①では区間を過大(過広?)に、アンチパターン②では過小(過狭?)に推定してしまっています。

ベルヌーイ分布における変化率の区間推定

コントロール群: p=0.10, 処置群: p=0.104として、ベルヌーイ分布からそれぞれ10000件のデータを生成します。サイト上のある機能について、接触したユーザー数がコントロール群で1000人、処置群で1040人、といったようなケースです。

推定区間が真の変化率(\(\frac{0.104}{0.10}-1=0.04\))を含んでいる割合は以下のようになりました。

• デルタ法: 0.9491
• アンチパターン①: 1.0
• アンチパターン②: 0.8171

推定の正しさの傾向としてはポワソン分布と同じですね。

正規分布における変化率の区間推定

正規分布についても試してみましょう。

コントロール群: mean=100, std=5, 処置群: mean=104, std=4でやります。

推定区間が真の変化率(\(\frac{104}{100}-1=0.04\))を含んでいる割合は以下のようになりました。

• デルタ法: 0.9487
• アンチパターン①: 1.0
• アンチパターン②: 0.8137

はい、こちらも推定の正しさの傾向は変わりませんでした。

結果

デルタ法による区間推定は、期待している信頼度で正しく推定できているといえます。

一方アンチパターンな推定は、区間を過大・過小に推定してしまいました。

まとめ

デルタ法を用いて、メトリクスの変化率を正しく区間推定する方法を紹介しました。

また合成データによる実験を通じて、デルタ法を用いた推定法と、アンチパターンとされる方法で推定した場合とで、それぞれ推定がどのようになるか確認しました。

以下が結論になります。

• • 比率の分散は、デルタ法によって正しく推定できる。
  • ABテストにおける変化率の区間推定は、\(CI = (\frac{\bar{Y_t}}{\bar{Y_c}}-1) \pm z_{\alpha/2}*\sqrt{\frac{\frac{1}{\bar{Y_c}^2}*var(Y_t)+\frac{\bar{Y_t}^2}{\bar{Y_c}^4}*var(Y_c)}{n}}\)で計算する。
  • 変化率の区間推定になっていなかったり、推定の過程で分散を正しく推定できていなかったりすると、区間を過大・過小に推定してしまう。以下はやりがちなので注意。
    • アンチパターン① (変化率の区間推定になっていない): \(CI = \frac{lower(\bar{Y_t})}{\bar{Y_c}} \sim \frac{upper(\bar{Y_t})}{\bar{Y_c}}\)としてしまう。
    • アンチパターン② (変化率の分散推定が正しくない): \(var(\Delta\%) = \frac{var(\Delta)}{\bar{Y_c}^2}\)としてしまう。

さいごに

グループ研究開発本部 AI研究開発室では、データサイエンティスト/機械学習エンジニアを募集しています。ビッグデータの解析業務などAI研究開発室にご興味を持って頂ける方がいらっしゃいましたら、ぜひ 募集職種一覧 からご応募をお願いします。皆さんのご応募をお待ちしています。

参考文献

[1] 『A/Bテスト実践ガイド 真のデータドリブンへ至る信用できる実験とは』

[2]『Applying the Delta Method in Metric Analytics: A Practical Guide with Novel Ideas』