こんにちは、リブセンスでデータサイエンティストをしている北原です。今回は前回の続きで二元分割表のベイズ推定を扱います。前回は行和と列和のいずれかが与えられる場合でしたが、今回は総度数のみが与えられる場合です。コードはRとStanです。なお、二元分割表のベイズ推定に関する別記事は下記をご参照ください。
総度数のみが与えられる場合のモデル
総度数のみが与えられる場合というのは、行も列も結果的に得られる変量のケースです。変量の分布に特徴が見られるかや、変量間に相関があるかを調べる時に使います。
例えば、弊社のサービスknewでビデオチャットを行った後のお互いの評価の関係を調べるのに使えます。knewではビデオチャット後に相手への満足度を11段階、評価を7段階で取得しています。ビデオチャットごとに男性から女性へ評価と女性から男性への評価を集計すると二元分割表ができます。いずれの評価も結果であり、ある時点ではビデオチャットの総数は決まっているので総度数のみが与えられる場合に該当します。この分割表を調べると、お互いの評価の関係性に特徴があるか、相関があるかなどを調べることができます。
総度数のみが与えられる場合では、個々のセルの度数を多項分布でモデル化します。$N$行$M$列の分割表の度数$x_{i,j}$を推定する場合、$N_a$を総度数、$i$行$j$列の確率を$\theta_{i,j}$として
となります。
行と列が独立しているならば周辺確率$q_i$ と$r_j$を
としたとき
となります。この関係性($\ref{relation}$)を利用して特徴のあるセルを見つけることもできます。
実装
では、Stanで推定してみましょう。モデル($\ref{model}$)を使ってセルの確率や度数を推定するだけでなく、($\ref{relation}$)を使ってセルの確率が周辺確率同士の積から乖離の大きい特徴的なセルを見つけます。そのため、周辺確率の推定も行います。
Stanコードは次の通りです。ここでは多項分布の事前分布にディリクレ分布を使っていますが、無情報事前分布を使うのであれば必須ではありません。入力データを2次元の分割表としているので一次元の配列として扱えるようにto_array_1d()関数を使っています。
// ct2_1.stan
data {
int<lower=1> N; // 行数
int<lower=1> M; // 列数
array[N, M] int<lower=0> X; // 分割表データ
}
transformed data {
int<lower=1> T; // 総度数
array[N] int<lower=0> X_r; // 行和
array[M] int<lower=0> X_c; // 列和
T = sum(to_array_1d(X));
for (i in 1:N) {
X_r[i] = sum(X[i,]);
}
for (j in 1:M) {
X_c[j] = sum(X[,j]);
}
}
parameters {
simplex[N * M] theta; // 多項分布のパラメータ(全セル)
simplex[N] theta_r; // 多項分布のパラメータ(行)
simplex[M] theta_c; // 多項分布のパラメータ(列)
}
model {
// パラメータの事前分布
theta ~ dirichlet(rep_vector(1.0, N * M));
theta_r ~ dirichlet(rep_vector(1.0, N));
theta_c ~ dirichlet(rep_vector(1.0, M));
to_array_1d(X) ~ multinomial(theta);
X_r ~ multinomial(theta_r);
X_c ~ multinomial(theta_c);
}
generated quantities {
matrix[N, M] diff;
array[N * M] int X_est;
// 度数の推定
X_est = multinomial_rng(theta, T);
// theta_{i,j} - q_i * r_j
for (i in 1:N) {
for (j in 1:M) {
diff[i, j] = theta[(i - 1) * M + j] - theta_r[i] * theta_c[j];
}
}
}
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) mod2_1 <- cmdstan_model("ct2_1.stan") fit2_1 <- mod2_1$sample( data = data_list, seed = 123 ) # セルの確率と周辺確率同士の積との差分(theta_{i,j} - q_i * r_j)が90%区間に0を含まないセル fit2_1$summary("diff") %>% filter(q5 * q95 > 0) %>% head(10) %>% mutate_at(vars(-variable), ~ round(., 2)) %>% kable() %>% print()
以下はセルの確率と周辺確率同士の積の差の90%信用区間が0を含んでいないセルを抽出した結果の一部です。適合度検定と異なりセル単位で調べられるのが利点です。
| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| diff[1,1] | 0.05 | 0.05 | 0.01 | 0.01 | 0.04 | 0.06 | 1 | 4765.11 | 2763.49 |
| diff[2,1] | 0.04 | 0.04 | 0.01 | 0.01 | 0.03 | 0.05 | 1 | 5111.12 | 2996.99 |
| diff[3,1] | 0.02 | 0.02 | 0.01 | 0.01 | 0.01 | 0.03 | 1 | 5738.20 | 3193.49 |
| diff[4,1] | 0.02 | 0.02 | 0.01 | 0.01 | 0.01 | 0.03 | 1 | 4792.35 | 2953.83 |
| diff[8,1] | -0.01 | -0.01 | 0.00 | 0.00 | -0.02 | -0.01 | 1 | 5105.75 | 3075.50 |
| diff[9,1] | -0.02 | -0.02 | 0.00 | 0.00 | -0.03 | -0.02 | 1 | 4714.39 | 3164.94 |
| diff[10,1] | -0.01 | -0.01 | 0.00 | 0.00 | -0.02 | -0.01 | 1 | 4938.30 | 2806.89 |
| diff[11,1] | -0.02 | -0.02 | 0.00 | 0.00 | -0.03 | -0.02 | 1 | 4973.46 | 3116.60 |
| diff[12,1] | -0.03 | -0.03 | 0.00 | 0.00 | -0.04 | -0.03 | 1 | 4924.92 | 3077.45 |
| diff[13,1] | -0.03 | -0.03 | 0.00 | 0.00 | -0.04 | -0.03 | 1 | 4519.26 | 2904.93 |
行和と列和のいずれかが与えられる場合との違い
前回の記事で扱った、行和と列和のいずれかが与えられる場合との違いは、分析者の問題設定次第であるところもあります。2変量が結果的な情報であっても、一方の変量について特定の結果が得られたということを属性とみなして、行和と列和のいずれかが与えられる場合のモデルを考えることも可能なように思います。問題設定を変えるとどのくらい推定結果が異なるか調べてみましょう。
総度数のみが与えられる場合と行和と列和のいずれかが与えられる場合とを比較するStanとRのコードは次の通りです。
// ct2_2.stan
data {
int<lower=1> N; // 行数
int<lower=1> M; // 列数
array[N, M] int<lower=0> X; // 分割表データ
}
transformed data {
int<lower=1> T; // 総度数
array[N] int<lower=0> X_r; // 行和
T = sum(to_array_1d(X));
for (i in 1:N) {
X_r[i] = sum(X[i,]);
}
}
parameters {
simplex[N * M] theta; // 多項分布のパラメータ(総度数のみが与えられる場合)
simplex[M] theta_r[N]; // 多項分布のパラメータ(行和と列和のいずれかが与えられる場合)
}
model {
// パラメータの事前分布
theta ~ dirichlet(rep_vector(1.0, N * M));
for (n in 1:N) {
theta_r[n] ~ dirichlet(rep_vector(1.0, M));
}
to_array_1d(X) ~ multinomial(theta);
for (n in 1:N) {
X[n] ~ multinomial(theta_r[n]);
}
}
generated quantities {
matrix[N, M] diff; // 推定度数の差(総度数のみが与えられる場合 - 行和と列和のいずれかが与えられる場合)
array[N * M] int X_est; // 推定度数(総度数のみが与えられる場合)
array[N, M] int X_r_est; // 推定度数(行和と列和のいずれかが与えられる場合)
for (n in 1:N) {
X_r_est[n] = multinomial_rng(theta_r[n], X_r[n]);
}
X_est = multinomial_rng(theta, T);
for (i in 1:N) {
for (j in 1:M) {
diff[i, j] = X_est[(i - 1) * M + j] - X_r_est[i, j];
}
}
}
mod2_2 <- cmdstan_model("ct2_2.stan") fit2_2 <- mod2_2$sample( data = data_list2, seed = 123, chains = 4, parallel_chains = 4 ) fit2_2$summary(c("theta[1]", "theta[5]", "theta[9]")) %>% rbind(fit2_2$summary("theta_r") %>% head(3)) %>% rbind(fit2_2$summary(c("X_est[1]", "X_est[5]", "X_est[9]"))) %>% rbind(fit2_2$summary("X_r_est") %>% head(3)) %>% rbind(fit2_2$summary("diff") %>% head(3)) %>% mutate_at(vars(-variable), ~ round(., 2)) %>% kable() %>% print()
総度数のみが与えられる場合と行和と列和のいずれかが与えられる場合とで、セルの確率、度数の推定値、度数の差を、分割表の1列目1~3行について出力した結果は次の通りです。セルの確率の推定値は両者で推定しているものが違うので異なりますが、今回のサンプルデータでは度数の推定値に大きな乖離はありません。ただし、総度数のみが与えられる場合のほうが標準偏差が大きく信用区間の幅も広くなっていることがわかります。総度数のみより行ごとの行和のほうが情報が多いので妥当な結果とも言えます。総度数のみが与えられる場合であっても、一方の結果変量の傾向がほぼ明らかな場合は、行和と列和のいずれかが与えられる場合とみなしたり、事前分布を指定して調整してもよいかもしれません。
| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| theta[1] | 0.09 | 0.09 | 0.01 | 0.01 | 0.08 | 0.10 | 1 | 8412.47 | 2747.98 |
| theta[5] | 0.07 | 0.07 | 0.01 | 0.01 | 0.06 | 0.08 | 1 | 6760.74 | 2655.30 |
| theta[9] | 0.04 | 0.04 | 0.00 | 0.00 | 0.04 | 0.05 | 1 | 6893.18 | 3242.64 |
| theta_r[1,1] | 0.87 | 0.87 | 0.02 | 0.03 | 0.83 | 0.91 | 1 | 6639.61 | 3022.40 |
| theta_r[2,1] | 0.80 | 0.80 | 0.03 | 0.03 | 0.74 | 0.85 | 1 | 8242.93 | 2254.09 |
| theta_r[3,1] | 0.70 | 0.70 | 0.04 | 0.04 | 0.62 | 0.77 | 1 | 8060.61 | 3037.86 |
| X_est[1] | 152.23 | 152.00 | 16.16 | 16.31 | 126.00 | 179.00 | 1 | 5183.52 | 3667.87 |
| X_est[5] | 121.36 | 121.00 | 14.95 | 14.83 | 98.00 | 146.05 | 1 | 5256.62 | 3841.17 |
| X_est[9] | 75.52 | 75.00 | 11.88 | 11.86 | 57.00 | 96.00 | 1 | 4956.25 | 3791.91 |
| X_r_est[1,1] | 153.76 | 154.00 | 6.17 | 5.93 | 143.00 | 163.00 | 1 | 4700.32 | 3817.66 |
| X_r_est[2,1] | 121.89 | 122.00 | 6.85 | 7.41 | 110.00 | 133.00 | 1 | 5287.65 | 3684.45 |
| X_r_est[3,1] | 75.22 | 75.00 | 6.55 | 5.93 | 64.00 | 86.00 | 1 | 5823.34 | 3755.41 |
| diff[1,1] | -1.54 | -2.00 | 17.38 | 16.31 | -30.00 | 28.00 | 1 | 4943.48 | 3748.22 |
| diff[2,1] | -0.53 | -1.00 | 16.48 | 16.31 | -27.00 | 27.00 | 1 | 5426.36 | 4103.59 |
| diff[3,1] | 0.31 | 0.00 | 13.56 | 13.34 | -21.00 | 23.00 | 1 | 4894.10 | 3660.11 |
まとめ
今回は、総度数のみが与えられる場合の二元分割表のベイズ推定を紹介しました。今回も初歩的な推定ですが、結果変量間の関係を調べたいことはよくあるので、すぐに使えるようにしておくと便利です。
