二元分割表のベイズ推定 - 行和と列和のいずれかが与えられる場合

こんにちは、リブセンスでデータサイエンティストをしている北原です。今回は、二元分割表のベイズ推定を扱います。実務のデータ分析ではクロス集計をよく使うと思いますが、その集計で得られた表を詳細に分析するときに使います。コードはRとStanです。なお、二元分割表のベイズ推定に関する別記事は下記をご参照ください。

analytics.livesense.co.jp

二元分割表を対象にした主な統計モデル
行和と列和のいずれかが与えられる場合のモデル
実装
まとめ
参考文献

二元分割表を対象にした主な統計モデル

　二元分割表がどのようにしてできあがったかによって統計モデルが異なります。まず(宮川, 青木, 2018)を参考に二元分割表の主な統計モデルを紹介します。

サンプリング	統計モデル	ユースケース
総度数のみが与えられる場合	多項分布	ある学生の2つの試験成績
行和と列和のいずれかが与えられる場合	多項分布(行または列ごと)	クラスごとの試験成績
行和と列和のいずれもが与えられる場合	多項超幾何分布	官能検査
総度数も与えられない場合	ポアソン分布	2因子実験

　総度数のみが与えられる場合というのは、行と列の双方が反応、つまり結果的な変量、であるケースです。ある学生を選定して、その学生たちに2つの試験を行う場合について考えてみましょう。学生数が固定なので試験成績の総数は一定で、どのような成績が多いかあるいは少ないかだけが異なります。また、どちらの試験の成績も結果的な変量です。そのため、総度数のみが与えられる場合に該当します。

　行和と列和のいずれかが与えられる場合というのは、行あるいは列はあらかじめ定められており、もう一方が反応であるケースです。クラスごとに試験を行う場合について考えてみましょう。クラスの人数は一定で、クラス内でどの成績が多いかあるいは少ないかだけが異なります。また、試験成績のみが結果的な変量です。そのため、行和と列和のいずれかが与えられる場合に該当します。

　行和と列和のいずれもが与えられる場合は、各行もしくは各列の度数が試験計画で決まり、もう一方は度数のみが知られている判定のようなケースです。例えば、コーヒーやワインを試飲して銘柄を当てるケースです。正解の銘柄数は決まっているので度数は固定です。一方で、各銘柄の数がわかっており判定はその数に合わせるため度数も固定です。正解の銘柄数も、判定される銘柄数も決まっており、間違え方のみが異なっているので、このケースに該当します。実務上、あまり適用例が思いつかないですが、フィッシャーの正確検定はこのケースに基づいて提唱されたようなので、理論上は重要性の高い問題設定です。

　総度数も与えられない場合は、行と列の双方が反応で、総度数も結果的に決まるケースです。何らかの不具合件数に関する2因子実験について考えてみましょう。各因子に関わる不具合は結果的に決まるので行和ともに反応です。不具合の総数も結果的に決まるので総度数は与えられません。そのため、総度数も与えられない場合に該当します。

　本記事では、これらのうち実務上、最も利用頻度が高いと思われる行和と列和のいずれかが与えられる場合について扱います。

行和と列和のいずれかが与えられる場合のモデル

　行和と列和のいずれかが与えられる場合というのは、一方が属性などで分類されていて行和もしくは列和が決まっていて、もう一方はその属性ごとの結果であるケースです。前節で紹介した4つのケースの中では一番利用頻度が高いのではないかと思います。属性が異なると、結果の比率がどうなるかを調べるときに使います。なお、行と列には対称性があるので以下では行和が与えられる場合のみを考えます。

　例えば、カスタマジャーニーにおける各フェーズでの離脱率がユーザー属性ごとに違うかなどを調べるのに使えます。求人サービスで年齢層によってどこで離脱するかを分析する場合を考えましょう。会員登録したうち応募にすら至らず離脱、応募後に採用に至らず離脱、内定後に採用に至らず離脱、無事採用それぞれの数を年齢層ごとに集計します。年齢層ごとにユーザー数は決まっているので行和は与えられていると考えられます。ただし、離脱率は年齢層によって異なるかもしれないので、その離脱率を推定します。

　行和が与えられる場合では、行ごとに多項分布でモデル化します。$N$行$M$列の分割表の度数$x_{i,j}$を推定する場合、$N_i$を$i$行の総度数、$i$行$j$列の確率を$\theta_{i,j}$として

$$ \begin{align} P(\mathbf{x}_i | \mathbf{\theta}_i) &= \frac{N_i!}{\prod_{j=1}^M x_{i,j}!} \prod_{j=1}^M \theta_{i,j}^{x_{i,j}} \\ \sum_{j=1}^M \theta_{i,j} &= 1 \end{align} $$

となります。

実装

　では、Stanで推定してみましょう。　

　Stanコードは次の通りです。ここでは多項分布の事前分布にディリクレ分布を使っていますが、無情報事前分布を使うのであれば必須ではありません。多項分布のパラメータをsimplex型にするところがポイントです。行和は確率$\theta_{i,j}$の推定だけであれば指定しなくともよいですが、度数を推定する場合は必要です。

// ct1.stan
data {
  int<lower=1> N;             // 行数
  int<lower=1> M;             // 列数
  array[N, M] int<lower=0> X; // 分割表データ
}

transformed data {
  array[N] int<lower=0> X_r; // Xの行ごとの総和
  for (i in 1:N) {
    X_r[i] = sum(X[i, ]); 
  }
}

parameters {
  simplex[M] theta[N]; // 多項分布のパラメータ
}

model {
  // 事前分布
  for (i in 1:N) {
    theta[i] ~ dirichlet(rep_vector(1.0, M));
  }

  // 多項分布
  for (i in 1:N) {
    X[i] ~ multinomial(theta[i]);
  }
}

generated quantities {
  array[N, M] int X_est;

  for (n in 1:N) {
    X_est[n] = multinomial_rng(theta[n], X_r[n]);
  }
}

　Rの呼び出しコードは次のとおりです。サンプルデータとしてfactoextraパッケージに入っているhousetasksデータを使っています。

library(tidyverse)
library(cmdstanr)
library(knitr)
library(factoextra)

data(housetasks)
dat <- housetasks %>% as.matrix()

data_list <- 
  list(
    N = nrow(dat),
    M = ncol(dat),
    X = dat)
mod <- cmdstan_model("ct1.stan")
fit <- mod$sample(
  data = data_list, 
  seed = 123
)
fit$summary() %>% head(10) %>% kable() %>% print()

以下のような出力が得られます。実際には分割表の度数を推定するより、着目する行同士の確率の差を推定して本当に差があるかを調べることが多いと思います。

variable	mean	median	sd	mad	q5	q95	rhat	ess_bulk	ess_tail
lp__	-1515.51	-1515.22	4.49	4.40	-1523.48	-1508.85	1	1940.66	2606.41
theta[1,1]	0.87	0.87	0.02	0.02	0.83	0.91	1	5340.01	2257.17
theta[2,1]	0.80	0.80	0.03	0.03	0.74	0.85	1	5012.71	2581.44
theta[3,1]	0.70	0.70	0.04	0.04	0.63	0.76	1	4923.69	3348.61
theta[4,1]	0.58	0.58	0.04	0.04	0.51	0.64	1	5409.83	3127.79
theta[5,1]	0.43	0.43	0.04	0.04	0.36	0.50	1	5976.73	2897.42
theta[6,1]	0.28	0.28	0.04	0.04	0.22	0.35	1	5272.16	3058.75
theta[7,1]	0.27	0.27	0.04	0.04	0.21	0.34	1	5001.32	2855.54
theta[8,1]	0.13	0.13	0.03	0.03	0.08	0.19	1	5270.87	3245.14
theta[9,1]	0.08	0.08	0.02	0.02	0.04	0.12	1	5810.50	2816.92