二元分割表のベイズ推定 - 対数線形モデル

　こんにちは、リブセンスでデータサイエンティストをしている北原です。今回も二元分割表のベイズ推定ですが、今回は一般化線形モデルの一つである対数線形モデルを扱います。コードはRとStanです。なお、二元分割表のベイズ推定に関する別記事は下記をご参照ください。

analytics.livesense.co.jp

対数線形モデル

　対数線形モデルは、期待度数の対数を各効果の線形和で表現したモデルです。総度数や行和、列和を与えず、ポアソン分布を仮定するのが一般的です。そのため、総度数も与えられない場合に対応したモデルとして使われます。二元分割表だけでなく、多元分割表の分析にも使えますが、本記事では二元分割表のみを扱います。

　行と列を独立と仮定した場合、定数項と主効果項のみで交互作用項を含まない加法モデルを使います。加法モデルの$i$行$j$列の度数$y_{i,j}$は、定数項$\mu$と行の主効果$\alpha_i$、列の主効果$\beta_j$を使って

$$ \begin{align} \log(y_{i,j}) = \mu + \alpha_i + \beta_j \label{add_model} \end{align} $$

となります。

　行と列を独立と仮定した場合、交互作用項を含むモデルを使います。全ての交互作用項を含むモデルは飽和モデルと呼ばれます。一部の交互作用がないことが明らかな場合は、対応する交互作用項をモデルから適宜除外します。交互作用項を含むモデルの$i$行$j$列の度数$y_{i,j}$は、交互作用項$(\alpha \beta)_{i,j}$を使って

$$ \begin{align} \log(y_{i,j}) = \mu + \alpha_i + \beta_j + (\alpha \beta)_{i, j} \label{sat_model} \end{align} $$

となります。

　対数線形モデルを使うことで、行と列が独立しているか、独立していない場合はどのような交互作用があるかを調べることができます。加法モデルと飽和モデルを比較し、両者に差がないならば、行と列は独立していると判断できます。独立していない場合は、飽和モデルの交互作用項を調べることで、どのような交互作用があるかがわかります。

実装

　では、Stanで推定してみましょう。ここでは加法モデルと飽和モデルの二つを推定します。

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

　ここでは、加法モデル、飽和モデルともに、パラメータの事前分布に0を平均とした正規分布を設定しています。これは、いずれの効果も存在しない確率が高いことを表しています。この事前分布を使っても0にならないパラメータに対応する効果は無視できないものであると考えられます。なお、飽和モデルは推定パラメータが多いため、無情報事前分布を使うと推定が不安定になりやすいので注意が必要です。

　実装上の注意点は、推定結果が不定にならないように、推定対象のパラメータ数を行数や列数より1小さくしているところです。ここでは行、列ともに添字番号が一番大きい効果は、主効果や交互作用のパラメータを0としたときに対応するようにしています。

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

parameters {
  // 飽和モデルのパラメータ
  real m_s;
  vector[N - 1] a_s;
  vector[M - 1] b_s;
  matrix[N - 1, M - 1] ab_s;
  real<lower=0> a_s_sig;
  real<lower=0> b_s_sig;
  real<lower=0> ab_s_sig;
  
  // 加法モデルのパラメータ
  real m_a;
  vector[N - 1] a_a;
  vector[M - 1] b_a;
  real<lower=0> a_a_sig;
  real<lower=0> b_a_sig;
}

model {
  // 飽和モデルの事前分布
  a_s ~ normal(0, a_s_sig);
  b_s ~ normal(0, b_s_sig);
  to_vector(ab_s) ~ normal(0, ab_s_sig);
  
  // 加法モデルの事前分布
  a_a ~ normal(0, a_a_sig);
  b_a ~ normal(0, b_a_sig);
  
  for (i in 1:(N - 1)) {
    for (j in 1:(M - 1)) {
      // 飽和モデル
      X[i, j] ~ poisson_log(m_s + a_s[i] + b_s[j] + ab_s[i, j]);
      // 加法モデル
      X[i, j] ~ poisson_log(m_a + a_a[i] + b_a[j]);
    }
  }
  for (i in 1:(N - 1)) {
    X[i, M] ~ poisson_log(m_s + a_s[i]);
    X[i, M] ~ poisson_log(m_a + a_a[i]);
  }
  for (j in 1:(M - 1)) {
    X[N, j] ~ poisson_log(m_s + b_s[j]);
    X[N, j] ~ poisson_log(m_a + b_a[j]);
  }
  X[N, M] ~ poisson_log(m_s);
  X[N, M] ~ poisson_log(m_a);
}

generated quantities {
  array[N, M] int<lower=0> X_s_est;   // 飽和モデル
  array[N, M] int<lower=0> X_a_est;   // 加法モデル
  
  for (i in 1:(N - 1)) {
    for (j in 1:(M - 1)) {
      X_s_est[i, j] = poisson_log_rng(m_s + a_s[i] + b_s[j] + ab_s[i, j]);
      X_a_est[i, j] = poisson_log_rng(m_a + a_a[i] + b_a[j]);
    }
  }
  for (i in 1:(N - 1)) {
    X_s_est[i, M] = poisson_log_rng(m_s + a_s[i]);
    X_a_est[i, M] = poisson_log_rng(m_a + a_a[i]);
  }
  for (j in 1:(M - 1)) {
    X_s_est[N, j] = poisson_log_rng(m_s + b_s[j]);
    X_a_est[N, j] = poisson_log_rng(m_a + b_a[j]);
  }
  X_s_est[N, M] = poisson_log_rng(m_s);
  X_a_est[N, M] = poisson_log_rng(m_a);
}

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

// このモデルは推定に時間がかかるので並列処理する
library(doParallel)
registerDoParallel(detectCores())
options(mc.cores = parallel::detectCores())

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

data_list <- 
  list(
    N = nrow(dat),
    M = ncol(dat),
    X = dat)
mod4 <- cmdstan_model("ct4.stan")
fit4 <- mod4$sample(
  data = data_list, 
  seed = 123
)
fit4$summary("a_s") %>% head(3) %>% 
  rbind(
    fit4$summary("b_s") %>% head(3) 
  ) %>% 
  rbind(
    fit4$summary("ab_s") %>% head(3) 
  ) %>% 
  mutate_at(vars(-variable), ~ round(., 2)) %>% 
  kable() %>% 
  print()
fit4$summary("X_s_est") %>% head(5) %>% 
  rbind(
    fit4$summary("X_a_est") %>% head(5) 
  ) %>% 
  mutate_at(vars(-variable), ~ round(., 2)) %>% 
  kable() %>% 
  print()

飽和モデルの行と列の主効果、交互作用を3つずつ出力すると以下のようになります。交互作用が強すぎる結果となっており、行と列が独立していないことがわかります。

variable	mean	median	sd	mad	q5	q95	rhat	ess_bulk	ess_tail
a_s[1]	-3.39	-3.36	0.47	0.46	-4.22	-2.71	1.00	2580.42	2298.99
a_s[2]	-3.36	-3.34	0.44	0.43	-4.13	-2.69	1.00	2453.33	2496.50
a_s[3]	-2.35	-2.34	0.28	0.28	-2.81	-1.90	1.00	2872.59	3274.16
b_s[1]	-4.08	-4.01	0.64	0.63	-5.24	-3.15	1.01	665.77	706.41
b_s[2]	-3.79	-3.74	0.55	0.54	-4.78	-2.97	1.01	741.40	1329.18
b_s[3]	-2.95	-2.92	0.38	0.37	-3.64	-2.38	1.00	977.60	1395.83
ab_s[1,1]	7.56	7.51	0.81	0.80	6.34	8.98	1.00	785.13	842.87
ab_s[2,1]	7.30	7.24	0.80	0.78	6.09	8.75	1.00	788.99	1101.15
ab_s[3,1]	5.80	5.77	0.71	0.69	4.71	7.05	1.00	749.83	879.95

度数の予測結果を飽和モデルと加法モデルから5つずつ出力すると以下のようになります。サンプルデータの度数は156、124、77、82、53なので、行と列を独立と仮定している加法モデルでは、今回のサンプルデータの度数を予測するのは難しいことがわかります。

variable	mean	median	sd	mad	q5	q95	rhat	ess_bulk	ess_tail
X_s_est[1,1]	155.43	155	17.89	17.79	128	186	1	4050.36	3754.34
X_s_est[2,1]	123.55	123	15.80	16.31	98	149	1	4188.74	3796.00
X_s_est[3,1]	76.38	76	12.09	11.86	57	97	1	4090.50	3673.82
X_s_est[4,1]	81.73	81	12.63	11.86	62	103	1	3745.98	3791.11
X_s_est[5,1]	52.33	52	10.26	10.38	36	70	1	4061.25	3683.46
X_a_est[1,1]	58.33	58	9.08	8.90	44	74	1	4129.89	3800.70
X_a_est[2,1]	51.52	51	8.20	8.90	38	66	1	4497.67	3850.48
X_a_est[3,1]	39.59	39	7.21	7.41	28	52	1	4570.38	3785.05
X_a_est[4,1]	48.12	48	7.94	7.41	36	61	1	4762.14	4007.21
X_a_est[5,1]	42.89	43	7.55	7.41	31	56	1	4295.31	3902.87

まとめ

　今回は、対数線形モデルを使った二元分割表のベイズ推定を紹介しました。対数線形モデルは、総度数や行和、列和が与えられない時に、交互作用の強さを直接的に推定するのに有用なモデルだと思います。

参考文献

[1] 宮川雅巳, 青木敏, 分割表の統計解析 : 二元表から多元表まで. 朝倉書店, 2018

[2] A. J. Dobson, An introduction to generalized linear models, second edition. in Chapman & Hall/CRC texts in statistical science. Taylor & Francis, 2010. (田中豊, 森川敏彦, 山中竹春, 冨田誠(訳). 一般化線形モデル入門. 共立出版, 2008)

[3] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, Bayesian Data Analysis, Third Edition. CRC Press, 2013.