二元分割表のベイズ推定 - 総度数も与えられない場合

　こんにちは、リブセンスでデータサイエンティストをしている北原です。今回は前々回、前回の続きで二元分割表のベイズ推定を扱います。前々回は行和と列和のいずれかが与えられる場合、前回は総度数のみが与えられる場合でしたが、今回は総度数も与えられない場合です。コードはRとStanです。なお、二元分割表のベイズ推定に関する別記事は下記をご参照ください。

analytics.livesense.co.jp

総度数も与えられない場合のモデル

　総度数も与えられない場合というのは、行も列のいずれかが結果的に得られる変量であり度数も定まっていないケースです。総度数のみが与えられる場合と同様に、変量の分布に特徴が見られるかや、変量間に相関があるかを調べる時に使います。

　前回の記事で、総度数のみが与えられる場合の例として弊社のサービスknewでビデオチャットを行った後のお互いの評価を紹介しましたが、総度数も与えられない場合として扱うこともできます。ビデオチャット件数が既知ならば総度数のみが与えられる場合に該当し、件数を既知としないならば総度数も与えられない場合に該当します。やりたいことや問題設定に応じて適宜使い分けます。すでに得られている分割表の分析をするのであれば、総度数のみが与えられる場合として扱うほうが結果の解釈が行いやすいと思います。一方で、今後得られる結果を予測する場合、件数すらわからないとするならば総度数も与えられない場合として扱います。一定の件数ビデオチャットが行われた時にどれくらいの度数が得られるかを予測するならば総度数のみが与えられる場合として扱います。

　総度数も与えられない場合では、個々のセルの度数をポアソン分布でモデル化します。$N$行$M$列の分割表の度数$x_{i,j}$を推定する場合、$i$行$j$列のパラメータを$\lambda_{i,j}$として

$$ \begin{align} P(\mathbf{x} | \mathbf{\lambda}) = \prod_{i=1}^N \prod_{j=1}^M \frac{\lambda^{x_{i,j}}_{i,j}}{x_{i,j}} e^{-\lambda_{i,j}} \label{poisson_model} \end{align} $$

となります。

　独立にポアソン分布にしたがう場合に総度数を条件付けると多項分布になります。総度数のみが与えられる場合と総度数も与えられない場合との違いは、統計モデルでも対応がとれているということです。導出はQiitaの記事などをご参照ください。

実装

　では、Stanで推定してみましょう。モデル($\ref{poisson_model}$)を使ってセルの度数を推定するだけでなく、総度数のみが与えられる場合との比較も行います。

　StanとRのコードは次の通りです。ここでは無情報事前分布としています。以前の記事と同様に、サンプルデータとしてfactoextraパッケージに入っているhousetasksデータを使っています。

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

transformed data {
  int<lower=1> T;             // 総度数
  T = sum(to_array_1d(X));
}

parameters {
  vector<lower=0>[N * M] lam; // 総度数も与えられない場合
  simplex[N * M] theta;       // 総度数のみが与えられる場合
}

model {
  to_array_1d(X) ~ poisson(lam);       // 総度数も与えられない場合
  to_array_1d(X) ~ multinomial(theta); // 総度数のみが与えられる場合
}

generated quantities {
  array[N * M] int<lower=0> X_p_est; // 総度数も与えられない場合
  array[N * M] int<lower=0> X_m_est; // 総度数のみが与えられる場合
  matrix[N, M] diff;                 // 推定度数の差(総度数のみが与えられる場合 - 総度数も与えられない場合)
  
  X_p_est = poisson_rng(lam);
  X_m_est = multinomial_rng(theta, T);
  for (i in 1:N) {
    for (j in 1:M) {
      diff[i, j] = X_m_est[(i - 1) * M + j] - X_p_est[(i - 1) * M + j];
    }
  }
}

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)
mod3 <- cmdstan_model("ct3.stan")
fit3 <- mod3$sample(
  data = data_list2, 
  seed = 123
)
fit3$summary("X_p_est") %>% head(3) %>% 
  rbind(fit3$summary("X_m_est") %>% head(3)) %>% 
  rbind(fit2_2$summary("diff") %>% head(3)) %>% 
  mutate_at(vars(-variable), ~ round(., 2)) %>% 
  kable() %>% 
  print()

総度数も与えられない場合と総度数のみが与えられる場合とで、度数の推定値と度数の差を、分割表の1行目1~3列について出力した結果は次の通りです。今回のサンプルデータでは度数の推定値に大きな乖離はありません。総度数も与えられない場合のほうが、標準偏差がわずかに大きく信用区間の幅もわずかに広くなっていることがわかります。総度数が与えられるかの違いなので、それほど大きな差はないようです。度数0が多いなど、度数をポアソン分布でモデル化するのが不適切な場合は別のモデルを使うことになります。そのような場合は差が大きくなることが予想されます。

variable	mean	median	sd	mad	q5	q95	rhat	ess_bulk	ess_tail
X_p_est[1]	156.47	156	17.89	17.79	128	187	1	5553.14	3689.25
X_p_est[2]	14.92	14	5.50	5.93	7	25	1	6215.01	3437.59
X_p_est[3]	2.96	2	2.45	1.48	0	8	1	5473.71	3560.23
X_m_est[1]	152.73	152	16.61	16.31	126	180	1	5866.96	3576.88
X_m_est[2]	14.49	14	5.34	4.45	7	24	1	4729.36	3824.86
X_m_est[3]	2.97	2	2.44	1.48	0	8	1	5560.17	3566.79
diff[1,1]	-1.54	-2	17.38	16.31	-30	28	1	4943.48	3748.22
diff[2,1]	-0.53	-1	16.48	16.31	-27	27	1	5426.36	4103.59
diff[3,1]	0.31	0	13.56	13.34	-21	23	1	4894.10	3660.11