ケーススタディ
ケーススタディとして、アフリカの牛の CBPP
病の発生を記述したlme4パッケージの cbpp
データを使用することにする。このデータセットには、period
（期間）、herd
（牛群を特定する因子）、incidence
（特定の牛群と期間における新しい病気の症例数）、および size
（特定の期間の開始時の牛群サイズ）の4つの変数が含まれている。
data("cbpp", package = "lme4")
head(cbpp)
herd incidence size period
1 1 2 14 1
2 1 3 12 2
3 1 4 9 3
4 1 0 5 4
5 2 3 22 1
6 2 1 18 2
まず最初に、ベースラインモデルとして単純な二項モデルを用いて
incidence を予測することにする。観測されたイベント数 \(y\) ( 私たちの場合は incidence
) と試行の総数 \(T\) (
size )
に対して、二項分布の確率質量関数は次のように定義される。
\[
P(y | T, p) = \binom{T}{y} p^{y} (1 - p)^{N-y}
\]
ここで、 \(p\)
はイベント確率である。古典的な2項モデルでは、 \(p\)
をlogitスケールで直接予測する。これは、各オブザベーション \(i\) に対して、成功確率 \(p_i\)
を次のように計算することを意味する。
\[
p_i = \frac{\exp(\eta_i)}{1 + \exp(\eta_i)}
\]
ここで、 \(\eta_i\)
は、オブザベーション \(i\)
の線形予測項である（brmsの線形予測項の詳細については、vignette("brms_overview")
を参照）。 period と変化する切片 herd によって
incidence を予測することは、brms
では直感的である。
fit1 <- brm(incidence | trials(size) ~ period + (1|herd),
data = cbpp, family = binomial())
要約出力では、発生確率は群によって大きく異なるが、period
の負の係数が示すように、時間の経過とともに減少することがわかる。
Family: binomial
Links: mu = logit
Formula: incidence | trials(size) ~ period + (1 | herd)
Data: cbpp (Number of observations: 56)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~herd (Number of levels: 15)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.75 0.22 0.40 1.27 1.00 1426 2101
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.40 0.27 -1.94 -0.90 1.00 1607 1615
period2 -0.99 0.31 -1.61 -0.40 1.00 4444 2863
period3 -1.14 0.33 -1.79 -0.51 1.00 4313 2706
period4 -1.61 0.43 -2.52 -0.83 1.00 4116 2396
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
二項モデルの欠点は、線形予測変数が考慮された後、その分散が \(\text{Var}(y_i) = T_i p_i (1 - p_i)\)
に固定されることである。この値を超えるすべての分散は、モデルによって考慮されないことがある。
このいわゆる「過剰分散」を扱う方法は複数あり、以下に説明する解決策は、brms
でカスタム・ファミリーを定義する方法の例として役立つ。
ベータ二項分布
ベータ二項モデルは、二項モデルに過分散を考慮するためのパラメータを追加した一般化モデルである。ベータ二項モデルでは、二項確率
\(p_i\)
を直接予測するのではなく、ハイパーパラメータ \(\alpha > 0\) と \(\beta > 0\) でベータ分布と仮定する。
\[
p_i \sim \text{Beta}(\alpha_i, \beta_i)
\]
\(\alpha\) と \(\beta\)
のパラメータはいずれも解釈が難しく、一般に回帰モデルでの使用は推奨されない。したがって、私たちはパラメータ
$$ および \(\phi > 0\)、これを \(\text{Beta2}\) と呼ぶことにする。
\[
\text{Beta2}(\mu, \phi) = \text{Beta}(\mu \phi, (1-mu) \phi)
\]
パラメータ \(\mu\) と \(\phi\)
は、それぞれ平均値と精度パラメータを指定する。を定義することによって
\[
\mu_i = \frac{path(\eta_i)}{1 + \exp(\eta_i)}.
\]
は、変換された線形予測変数によって期待される確率を（元の2項モデルのように）予測するが、パラメータ
\(\phi\)
によって、潜在的な過剰分散を考慮する。
カスタムファミリーモデルのあてはめ
ベータ二項分布は、現在では brms
でネイティブにサポートされているが、ここでは例として、custom_family
関数を用いて自分自身で定義することにする。この関数では、分布の名前、パラメータの名前（ここでは
mu と phi
）、対応するリンク関数（パラメータが予測される場合のみ適用）、理論下限と上限（パラメータが予測されない場合のみ適用）、分布が離散か連続かの情報、最後にパラメータ以外の追加の変数が分布に渡される必要があるかどうかを指定する。ベータ二項分布の例では、この結果、次のカスタムファミリーが作成される。
beta_binomial2 <- custom_family(
"beta_binomial2", dpars = c("mu", "phi"),
links = c("logit", "log"),
lb = c(0, 0), ub = c(1, NA),
type = "int", vars = "vint1[n]"
)
試行回数を含む変数の名前 vint1
は、後述するように任意に選ばれたものではない。次に、分布がStan自体で定義されていない場合、関連するStan関数を提供する必要がある。beta_binomial2
分布については、順序分布である beta_binomial
分布がすでに実装されているので、これは簡単なことである。
stan_funs <- "
real beta_binomial2_lpmf(int y, real mu, real phi, int T) {
return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);
}
int beta_binomial2_rng(real mu, real phi, int T) {
return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);
}
"
モデルフィッティングには、beta_binomial2_lpmf
だけが必要であるが、beta_binomial2_rng
は、後処理をするときに便利だろう。以下のように定義する。
stanvars <- stanvar(scode = stan_funs, block = "functions")
試行回数(整数変数)の情報を提供するために、カスタムファミリでのみ使用できる追加引数
vint()
を使用するつもりである。同様に、実データの追加ベクタを含める必要がある場合は、vreal()
を使用する。実は、この特別な例では、基本的な二項モデルのように、vint()
の代わりに、よりエレガントに追加引数 trials()
を適用することができる。しかし、この vignette
は、このトピックの一般的な概観を与えることを目的としているので、より一般的な方法を採用することにする。
これで、カスタムベータ二項モデルに適合させるためのすべてのコンポーネントが揃いた。
fit2 <- brm(
incidence | vint(size) ~ period + (1|herd), data = cbpp,
family = beta_binomial2, stanvars = stanvars
)
要約出力を見ると、period
の係数の不確実性が基本二項モデルよりいくらか大きいことがわかる。これは、モデルに過分散パラメータ
phi
を含めた結果である。それを除けば、結果はかなり似ている。
Family: beta_binomial2
Links: mu = logit; phi = identity
Formula: incidence | vint(size) ~ period + (1 | herd)
Data: cbpp (Number of observations: 56)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~herd (Number of levels: 15)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.38 0.25 0.02 0.95 1.00 930 1751
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.34 0.25 -1.85 -0.85 1.00 3729 2512
period2 -0.99 0.41 -1.83 -0.18 1.00 4421 3141
period3 -1.27 0.46 -2.23 -0.42 1.00 3848 2889
period4 -1.54 0.52 -2.60 -0.57 1.00 3525 2576
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 17.27 14.79 5.45 58.40 1.00 1782 1318
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
ポスト処理 カスタムファミリーモデル
summary や plot
などのいくつかの後処理方法は、カスタムファミリーモデルにそのまま適用できる。しかし、特に重要な3つのメソッドがあり、これらはユーザーによる追加入力が必要である。これらは、posterior_epred、posterior_predict、log_lik
で、それぞれ予測平均値、予測応答値、対数尤度値を計算する。これらはそれ自体に関係があるだけでなく、他の多くの後処理方法の基礎となるものである。例えば、私たちは、loo
の方法で実装された近似的なleave-one-out cross-validation
によって、二項モデルの適合性とベータ二項モデルの適合性を比較することに興味があるだろうが、そのためには
log_lik が動作していることが必要である。
あるファミリーの log_lik
関数は次のような名前でなければならない。
log_lik_<family-name> と名付け、2つの引数
i (観測値を示す) と prep
を持つ。prep
がどのように作られるかはあまり気にする必要はない
(興味があれば、prepare_predictions
関数をチェックしてみてみよう)。代わりに知っておくべきことは、パラメータはスロット
dpars に、データはスロット data
に格納されるということである。
一般的に、パラメータは、予測される場合（私たちの例では mu
）には \(S \times N\)
行列（事後ドローの数 \(S =\)
とオブザベーションの数 \(N =\)
）、予測されない場合（ phi ）にはサイズ \(N\) のベクトルの形式を取る。
R で完全な対数尤度関数を直接定義することもでくし、自己定義された
Stan
関数を公開し、それを適用することも可能である。後者の方が便利であるが、前者の方が安定しており、brms
をベースにした他のRパッケージでカスタムファミリーを実装する場合の唯一の選択肢となる。
今回の vignette では、後者の方法を採用する。
expose_functions(fit2, vectorize = TRUE)
を作成し、必要な log_lik
関数を数行のコードで定義する。
log_lik_beta_binomial2 <- function(i, prep) {
mu <- brms::get_dpar(prep, "mu", i = i)
phi <- brms::get_dpar(prep, "phi", i = i)
trials <- prep$data$vint1[i]
y <- prep$data$Y[i]
beta_binomial2_lpmf(y, mu, phi, trials)
}
get_dpar
関数は、分布パラメータが各行で別々に予測される場合と、フィット全体について同じである場合の両方を扱うために必要な変換を行う。
これによって、log_lik
を必要とするすべての後処理方法が同様に機能するようになる。例えば、モデルの比較は、単純に
Output of model 'fit1':
Computed from 4000 by 56 log-likelihood matrix
Estimate SE
elpd_loo -99.6 10.0
p_loo 21.8 4.2
looic 199.2 20.1
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 43 76.8% 1299
(0.5, 0.7] (ok) 10 17.9% 215
(0.7, 1] (bad) 3 5.4% 81
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.
Output of model 'fit2':
Computed from 4000 by 56 log-likelihood matrix
Estimate SE
elpd_loo -94.7 8.2
p_loo 10.5 1.9
looic 189.4 16.4
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 49 87.5% 398
(0.5, 0.7] (ok) 6 10.7% 784
(0.7, 1] (bad) 1 1.8% 184
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.
Model comparisons:
elpd_diff se_diff
fit2 0.0 0.0
fit1 -4.9 4.4
ELPD
の値が大きいほど適合度が高いので、β二項モデルはやや適合度が高いことがわかるが、対応する標準誤差から、その差はそれほど大きくはないことが明らかになった。
次に、posterior_predict 方式に必要な関数を定義する。
posterior_predict_beta_binomial2 <- function(i, prep, ...) {
mu <- brms::get_dpar(prep, "mu", i = i)
phi <- brms::get_dpar(prep, "phi", i = i)
trials <- prep$data$vint1[i]
beta_binomial2_rng(mu, phi, trials)
}
posterior_predict
関数は、対数尤度値の代わりに応答のランダムドローを作成していることを除けば、対応する
log_lik
関数とほとんど同じように見える。ここでも、公開された
便宜上、Stan
関数を使用する。posterior_predict
関数を使用しない場合でも、...
の引数を追加するようにしてみよう。いくつかのファミリーでは追加の引数が必要だことがある。posterior_predict
が動作することで、例えば事後予測的なチェックに従事することができる。
posterior_epred 関数を定義するとき、それは
prep
引数だけを持っていて、一度にすべてのオブザベーションの平均応答値を計算する必要があることを心に留めておく必要がある。ベータ二項分布の平均は
\(\text{E}(y) = \mu T\)
ですので、対応する posterior_epred
関数の定義はそれほど複雑ではないが、パラメータとデータの次元を揃える必要がある。
posterior_epred_beta_binomial2 <- function(prep) {
mu <- brms::get_dpar(prep, "mu")
trials <- prep$data$vint1
trials <- matrix(trials, nrow = nrow(mu), ncol = ncol(mu), byrow = TRUE)
mu * trials
}
posterior_epred に直接依存する後処理法は
conditional_effects
であり、予測因子の効果を視覚化することができる。
conditional_effects(fit2, conditions = data.frame(size = 1))
解釈を容易にするため、size を 1
に設定し、上記のプロットの Y 軸が確率を示すようにした。