イントロダクション
この vignette は、brms
を使用して非線形マルチレベル・モデルを適合させる方法について紹介する。非線形モデルは非常に柔軟で強力であるが、一般的な線形モデルよりもモデル仕様とプライオに注意が必要である。グループレベルの効果をとりあえず無視し、観測値
\(n\) の一般化線形モデルの予測項 \(\eta_n\)
は以下のように書くことができる。
\[\eta_n = \sum_{i = 1}^K b_i
x_{ni}\]
ここで、 \(b_i\) は予測変数 \(i\) の回帰係数、 \(x_{ni}\) はオブザベーション \(n\) に対する予測変数 \(i\)
のデータである。これは交互作用項や他のさまざまなデータ変換も含んでいる。しかし、
\(\eta_n\) の構造は、回帰係数 \(b_i\)
にいくつかの予測変数の値を掛けて合計するという意味で、常に線形である。このことは、仮想的な予測変数項
\[\eta_n = b_1 \exp(b_2 x_n)\]
は、もはや線形予測変数ではなく、一般化線形モデルの古典的な手法を使って適合させることはできないだろう。したがって、私たちはより一般的なモデル・クラスを必要とし、それを
非線形
モデルと呼ぶことにする。「非線形」という用語は、応答変数の仮定された分布について何も言っていないことに注意。特に、非線形予測項をあらゆる種類の応答分布に適用することができるので、「正規分布ではない」という意味ではない（brms
で利用可能な応答分布の詳細については、vignette("brms_families")
を参照）。
簡単な非線型モデル
まず、シミュレーションデータを使った簡単な例から始める。
b <- c(2, 0.75)
x <- rnorm(100)
y <- rnorm(100, mean = b[1] * exp(b[2] * x))
dat1 <- data.frame(x, y)
前述のように、一般化線形モデルを用いて \(b\)
を推定することはできないので、あえて非線形モデルを指定する。
prior1 <- prior(normal(1, 2), nlpar = "b1") +
prior(normal(0, 2), nlpar = "b2")
fit1 <- brm(bf(y ~ b1 * exp(b2 * x), b1 + b2 ~ 1, nl = TRUE),
data = dat1, prior = prior1)
上記のコードを見ると、まず明らかになるのは、formula
の構文を変更して、予測子（すなわち x
）とパラメータ（すなわち b1 と b2
）を含む非線形式を、bf
の呼び出しでラップして表示すようにしたことである。これは、予測変数のみが与えられ、パラメータが暗示的である古典的な
R 式と対照的である。引数 b1 + b2 ~ 1
には2つの目的がある。まず、formula
のどの変数がパラメータであるかという情報を提供し、次に、各パラメータの線形予測子の項を指定する。実際、非線形パラメータは、パラメータそのものというよりも、線形予測変数のためのプレースホルダと考えるべきである（以下の例も参照）。今回のケースでは、b1
と b2 を予測する他の変数がないので、上記のモデル式に \(b_1\) と \(b_2\)
の私たちの推定値を表す切片を当てはめるだけである。b1 + b2 ~ 1
という式は b1 ~ 1, b2 ~ 1
の短縮形であり、複数の非線形パラメータが同じ式を共有している場合に使用できる。nl = TRUE
を設定すると、brms
にその式が非線形として扱われるようになる。
一般化線形モデルとは対照的に、非線形モデルを同定するためには、母集団レベルのパラメータ（すなわち「固定効果」）に対する事前分布がしばしば必須となる。
そのため、brms
では、これらの事前分布を明示的に指定する必要がある。この例では、（母集団レベルの切片）
b1 に対して normal(1, 2)
の事前分布を用い、（母集団レベルの切片） b2 に対しては
normal(0, 2)
の事前分布を用いた。事前分布を設定することは、あらゆる種類のモデル、特に非線形モデルにおいて、自明な作業ではないので、常に適切な事前分布を考えるために時間をかけるべきである。非線形モデルを初めてフィットさせた後、異なるMCMCチェーンが異なる事後領域に収束するのを観察すると、かなり頻繁にプライアを変更することを余儀なくされることがある。これは同定問題の明確な兆候であり、一つの解決策はより強い（すなわち、より狭い）プライヤーを設定することである。
適合したモデルの要約を得るために、次のように適用する。
Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ b1 * exp(b2 * x)
b1 ~ 1
b2 ~ 1
Data: dat1 (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b1_Intercept 1.90 0.11 1.67 2.13 1.00 1605 1824
b2_Intercept 0.81 0.04 0.73 0.90 1.00 1521 1801
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.95 0.07 0.83 1.10 1.00 2329 1785
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).
plot(conditional_effects(fit1), points = TRUE)
summary
メソッドによると、私たちはかなりうまく真のパラメータ値を回復することができたことがわかる。plot
メソッドによると、私たちの MCMC チェーンはうまく収束し、同じ posterior
になった。conditional_effects
メソッドはモデルが示唆する（model
implied）（非線形）回帰直線を可視化する。
また、この非線形モデルを古典的な線形モデルと比較することにも興味がある。
fit2 <- brm(y ~ x, data = dat1)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ x
Data: dat1 (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.63 0.13 2.37 2.88 1.00 4173 3084
x 1.93 0.12 1.69 2.18 1.00 4046 2900
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.26 0.09 1.10 1.45 1.00 4060 3240
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).
モデルの適合度を調査・比較するために、バックエンドの
bayesplot
パッケージを利用したグラフィカルな事後予測チェックを適用することができる。
また、リーブワンアウトクロスバリデーションを用いて、モデルの適合度を簡単に比較することができる。
Output of model 'fit1':
Computed from 4000 by 100 log-likelihood matrix
Estimate SE
elpd_loo -137.9 7.5
p_loo 2.6 0.6
looic 275.9 15.0
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
Output of model 'fit2':
Computed from 4000 by 100 log-likelihood matrix
Estimate SE
elpd_loo -166.8 8.7
p_loo 4.1 1.2
looic 333.7 17.4
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
Model comparisons:
elpd_diff se_diff
fit1 0.0 0.0
fit2 -28.9 8.4
LOOIC
の値が小さいほどモデルの適合度が高いことを示すので、非線形モデルがよりデータに適合していることはすぐにわかる。もちろん、まさにそのモデルからデータをシミュレートしたので、これはあまり驚くことではないのである。
実世界の非線型モデル
Markus
Gesmannは自身のブログで、異なる起源年に由来する累積保険金支払額の経時的な伸びを予測している（https://www.magesblog.com/post/2015-11-03-loss-developments-via-growth-curves-and/
参照）。
ここでは、彼のモデルを少し簡略化したものを使って説明する。それは次のようなものである。
\[cum_{AY, dev} \sim N(\mu_{AY, dev},
\sigma)\] \[\mu_{AY, dev} = ult_{AY}
\left(1 - \exp\left(- \left( \frac{dev}{\theta} \right)^\omega \right)
\right)\]
累積保険金 \(cum\)
は時間とともに増加するが、この依存性を変数 \(dev\) を使ってモデル化する。さらに、 \(ult_{AY}\)
は、毎年の事故の最終損失額（推定値）である。これは、累積損失の成長を担うパラメータ
\(\theta\) および \(\omega\)
とともに、私たちのフレームワークにおける非線形パラメータを構成し、年度をまたいで同じであると仮定している。このデータは
brms と一緒に出力されている。
AY dev cum premium
1 1991 6 357.848 10000
2 1991 18 1124.788 10000
3 1991 30 1735.330 10000
4 1991 42 2182.708 10000
5 1991 54 2745.596 10000
6 1991 66 3319.994 10000
を提案し、提案したモデルを非線形brmsモデルに変換する。
fit_loss <- brm(
bf(cum ~ ult * (1 - exp(-(dev/theta)^omega)),
ult ~ 1 + (1|AY), omega ~ 1, theta ~ 1,
nl = TRUE),
data = loss, family = gaussian(),
prior = c(
prior(normal(5000, 1000), nlpar = "ult"),
prior(normal(1, 2), nlpar = "omega"),
prior(normal(45, 10), nlpar = "theta")
),
control = list(adapt_delta = 0.9)
)
最終損失額 ult について、事故年（変数 AY
）のグループレベルの効果を推定した。これはまた、非線形パラメータが実際には線形予測変数のプレースホルダーであることをうまく示している。ult
の場合は、年ごとに変化する切片だけを含んでいる。ここでも、母集団レベルの効果に関する事前分布が必要であり、このモデルでは、識別可能性を確保するために、実際には必須である。私たちは、よく知られた方法を用いて、モデルを要約する。
Family: gaussian
Links: mu = identity; sigma = identity
Formula: cum ~ ult * (1 - exp(-(dev/theta)^omega))
ult ~ 1 + (1 | AY)
omega ~ 1
theta ~ 1
Data: loss (Number of observations: 55)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~AY (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(ult_Intercept) 740.19 228.51 420.08 1287.53 1.00 1182 1854
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
ult_Intercept 5298.17 290.96 4757.71 5893.84 1.00 1168 1395
omega_Intercept 1.34 0.05 1.24 1.43 1.00 2465 2624
theta_Intercept 46.18 2.13 42.40 50.90 1.00 2075 2079
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 140.54 15.86 112.98 176.93 1.00 2683 2299
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).
plot(fit_loss, N = 3, ask = FALSE)
conditional_effects(fit_loss)
次に、限界効果を各年度に分けて示す。
conditions <- data.frame(AY = unique(loss$AY))
rownames(conditions) <- unique(loss$AY)
me_loss <- conditional_effects(
fit_loss, conditions = conditions,
re_formula = NULL, method = "predict"
)
plot(me_loss, ncol = 5, points = TRUE)
例えば、自然災害が特定の年にしか発生しないなどの理由で、事故年によって累積損失にばらつきがあることがわかる。
さらに、予測される累積損失の不確実性は、利用可能なデータポイントが少ない後年ほど大きいことが分かる。このデータセットに関するより詳細な議論については、Gesmann
& Morris (2020)の4.5節を参照。
項目反応モデルの応用
3つ目の例として、brms
の非線形モデルの枠組みを使って、より高度な項目反応モデルをモデル化する方法を示したいと思われる。簡単のために、1つの強制選択項目があり、3つの選択肢のうち1つだけが正しいとする。私たちの応答変数は、人がその項目に正しく答えるか(1)、答えないか(0)である。人は、その項目に正しく答える能力に差があると仮定する。しかし、すべての人が推測だけでその項目を正解する確率は33%である。そこで、このような状況を反映したデータをシミュレートする。
inv_logit <- function(x) 1 / (1 + exp(-x))
ability <- rnorm(300)
p <- 0.33 + 0.67 * inv_logit(ability)
answer <- ifelse(runif(300, 0, 1) < p, 1, 0)
dat_ir <- data.frame(ability, answer)
最も基本的な項目反応モデルは、単純なロジスティック回帰モデルと同等である。
fit_ir1 <- brm(answer ~ ability, data = dat_ir, family = bernoulli())
しかし、このモデルでは推測確率を完全に無視しているため、偏った推定値や予測値になる可能性が高い。
Family: bernoulli
Links: mu = logit
Formula: answer ~ ability
Data: dat_ir (Number of observations: 300)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.54 0.13 0.28 0.80 1.00 2693 2405
ability 0.79 0.14 0.53 1.06 1.00 3560 2720
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).
plot(conditional_effects(fit_ir1), points = TRUE)
推測確率を取り入れたより洗練されたアプローチは次のようになる。
fit_ir2 <- brm(
bf(answer ~ 0.33 + 0.67 * inv_logit(eta),
eta ~ ability, nl = TRUE),
data = dat_ir, family = bernoulli("identity"),
prior = prior(normal(0, 5), nlpar = "eta")
)
bernoulli ファミリーのリンク関数を identity
に設定することは非常に重要で、そうしないと2つのリンク関数が適用されることになる。これは、非線形予測項にはすでに望ましいリンク関数(
0.33 + 0.67 * inv_logit
)が含まれているのに、bernoulli
ファミリーはその上にデフォルトの logit
リンクを適用してしまうことがある。これはもちろん、奇妙で解釈しがたい結果につながる。
したがって、非線形予測項がすでに希望のリンク関数を含んでいる場合は、必ずリンク関数を
identity に設定してみよう。
Family: bernoulli
Links: mu = identity
Formula: answer ~ 0.33 + 0.67 * inv_logit(eta)
eta ~ ability
Data: dat_ir (Number of observations: 300)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
eta_Intercept -0.43 0.24 -0.94 0.01 1.00 1974 1759
eta_ability 1.37 0.29 0.86 1.99 1.00 1814 1838
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).
plot(conditional_effects(fit_ir2), points = TRUE)
リーブワンアウトクロスバリデーション（leave-one-out
cross-validation）によるモデル適合度の比較
Output of model 'fit_ir1':
Computed from 4000 by 300 log-likelihood matrix
Estimate SE
elpd_loo -184.9 6.6
p_loo 2.1 0.2
looic 369.9 13.2
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
Output of model 'fit_ir2':
Computed from 4000 by 300 log-likelihood matrix
Estimate SE
elpd_loo -184.6 6.8
p_loo 2.3 0.4
looic 369.2 13.7
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
Model comparisons:
elpd_diff se_diff
fit_ir2 0.0 0.0
fit_ir1 -0.3 1.9
しかし、最初のモデルの予測は、低い能力値に対する推測確率を下回っている可能性があり、まだ誤解を招く可能性があることを覚えておくと良い。さて、推測確率がわからず、データから推定したいとする。これは、前のモデルを少し変えるだけで簡単にできる。
fit_ir3 <- brm(
bf(answer ~ guess + (1 - guess) * inv_logit(eta),
eta ~ 0 + ability, guess ~ 1, nl = TRUE),
data = dat_ir, family = bernoulli("identity"),
prior = c(
prior(normal(0, 5), nlpar = "eta"),
prior(beta(1, 1), nlpar = "guess", lb = 0, ub = 1)
)
)
ここで、推測確率を非線形パラメータとしてモデル化し、それが区間 \([0, 1]\) . eta
の切片を推定しなかったのは、推定された推測パラメータにバイアスがかかるからである（試してみてほしい。これは、非線形モデルにおいていかに注意しなければならないかを示す良い例である）。
Family: bernoulli
Links: mu = identity
Formula: answer ~ guess + (1 - guess) * inv_logit(eta)
eta ~ 0 + ability
guess ~ 1
Data: dat_ir (Number of observations: 300)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
eta_ability 1.10 0.22 0.70 1.55 1.00 2634 2470
guess_Intercept 0.24 0.05 0.14 0.33 1.00 2308 1680
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).
plot(conditional_effects(fit_ir3), points = TRUE)
その結果、この非線形モデルでシミュレーションしたモデルのパラメータを回復することができた。もちろん、実際の項目応答データには複数の項目があり、項目や人ごとに複数の観測値があるため、項目や人の変動性を考慮する（例えば、切片を変化させたマルチレベルモデルを使用する）必要が出てくるのである。幸いなことに、これはすべて
brms の非線形フレームワークの中で行うことができ、この
vignette が良い出発点になることを望む。