Estimating Distributional Models with
brms
Paul Bürkner
2022-12-14
イントロダクション
この vignette は、brms
を使って分布回帰モデルを適合させる方法について紹介するものである。分布モデルという用語は、想定される応答分布のすべてのパラメータに対して予測項を指定することができるモデルを指すために使用する。回帰モデルの実装の大部分では、応答分布の位置パラメータ（通常は平均）だけが予測変数と対応する回帰パラメータに依存する。他のパラメータ（たとえば、スケールまたは形状パラメータ）は、オブザベーション間で一定であると仮定して、補助パラメータとして推定される。この仮定は、回帰モデルを適用するほとんどの研究者が、しばしば（私の経験では）それを緩和する可能性に気づかないほど一般的である。この仮定を緩和すると、モデルの複雑さが大幅に増加し、モデルの適合が困難になるため、これは理解できることである。幸い、brms
はバックエンドに Stan
を使用しており、ベイズモデルを推定するための非常に柔軟で強力なツールなので、モデルの複雑さはそれほど問題ではない。
正規分布の応答変数があるとする。そして、基本的な線形回帰では、正規分布の平均パラメータ
\(\mu\) に対して予測項 \(\eta_{\mu}\)
を指定する。正規分布の2番目のパラメータ–残差標準偏差 \(\sigma\)
–は、オブザベーション間で一定であると仮定される。私たちは \(\sigma\)
を推定するが、それを予測しようとはしない。しかし、分布モデルでは、予測項
\(\eta_{\mu}\) に加えて、 \(\sigma\) の予測項 \(\eta_{\sigma}\)
を指定することによって、まさにこれを実行する。グループ-レベル効果を無視し、オブザベーション
\(n\) のパラメータ \(\theta\)
の線形予測項は、次のような形式をとる。
\[\eta_{\theta n} = \sum_{i =
1}^{K_{\theta}} b_{\theta i} x_{\theta i n}\]
ここで \(x_{\theta i n}\)
は、オブザベーション \(n\)
に対するパラメータ \(\theta\) の \(i\) 番目の予測変数の値を示し、 \(b_{\theta i}\) は、パラメータ \(\theta\) の \(i\) 番目の回帰係数を示す。応答変数 \(y\)
を持つ分布正規モデルは、次のように書くことができる。
\[y_n \sim \mathcal{N}\left(\eta_{\mu n},
\, \exp(\eta_{\sigma n}) \right)\]
\(\sigma\)
は標準偏差を構成するため正の値のみを取り、線形予測変数は任意の実数であることを反映し、
\(\eta_{\sigma}\)
の周りに指数関数を使用した。
簡単な分布モデル
不等分散モデルは、分布モデルの中で最も単純であるが、それにもかかわらず非常に重要な応用例である。2つの患者グループがあるとする。一方のグループは治療（例えば、抗うつ薬）を受け、もう一方のグループはプラセボを受ける。治療がすべての患者に等しく効くとは限らないので、治療の数週間後には、治療群の症状分散がプラセボ群の症状分散よりも大きくなる可能性がある。簡単のために、治療後の値だけを調査するとする。
group <- rep(c("treat", "placebo"), each = 30)
symptom_post <- c(rnorm(30, mean = 1, sd = 2), rnorm(30, mean = 0, sd = 1))
dat1 <- data.frame(group, symptom_post)
head(dat1)
group symptom_post
1 treat -1.223302
2 treat 5.297133
3 treat 2.117939
4 treat 4.628694
5 treat 2.527661
6 treat -2.478216
次のモデルは、正規の応答分布の平均と残留標準偏差の両方について、group
の影響を推定する。
fit1 <- brm(bf(symptom_post ~ group, sigma ~ group),
data = dat1, family = gaussian())
有用な要約統計量やプロットは、以下の方法で取得できる。
summary(fit1)
plot(fit1, N = 2, ask = FALSE)
plot(conditional_effects(fit1), points = TRUE)
対数スケールでの2つの残留標準偏差の対比である母集団レベル効果
sigma_grouptreat
は、両グループの分散が確かに異なることを明らかにする。この印象は、group
の conditional_effects
を見ると確認できる。さらに一歩進んで、hypothesis
の方法を用いて、元のスケールでの残留標準偏差を計算することができる。
hyp <- c("exp(sigma_Intercept) = 0",
"exp(sigma_Intercept + sigma_grouptreat) = 0")
hypothesis(fit1, hyp)
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (exp(sigma_Interc... = 0 1.17 0.16 0.91 1.53 NA NA *
2 (exp(sigma_Interc... = 0 2.55 0.35 1.98 3.36 NA NA *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
また、直接比較し、その差の事後分布をプロットすることもできる。
hyp <- "exp(sigma_Intercept + sigma_grouptreat) > exp(sigma_Intercept)"
(hyp <- hypothesis(fit1, hyp))
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (exp(sigma_Interc... > 0 1.39 0.38 0.81 2.05 Inf 1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
確かに、治療群の残留標準偏差は、プラセボ群より大きいようだ。しかも、この差の大きさは、私たちがデータシミュレーションに入れた値から予想されるものとほぼ同じである。
ゼロインフレートモデル
分布回帰の枠組みのもう1つの重要なアプリケーションは、いわゆるゼロインフレート・モデルである。これらのモデルは、応答変数に当然予想されるよりも多くのゼロがあるときに役立つ。たとえば、人々が1日に吸うタバコの数を予測しようとして、非喫煙者も含めると、適切にモデル化されないと、パラメータの推定をひどく歪める可能性のある膨大な量のゼロが存在することになる。ここでは、様々な集団による魚の捕獲数を扱った例を考えてみる。UCLAのホームページ(˶‾᷄
-̫
‾᷅˵)では、次のように説明されている。「州の野生生物学者は、ある州立公園で釣り人がどれだけの魚を釣っているかをモデル化したいと考えている。訪問者には、滞在時間、グループの人数、子供の有無、釣った魚の数などを聞いている。釣りをしない人もいるが、釣りをしたかしないかのデータはない。釣りをした人の中には、釣りをしなかった人もいるので、データにはゼロが余っている。”
zinb <- read.csv("https://paul-buerkner.github.io/data/fish.csv")
head(zinb)
nofish livebait camper persons child xb zg count
1 1 0 0 1 0 -0.8963146 3.0504048 0
2 0 1 1 1 0 -0.5583450 1.7461489 0
3 0 1 0 1 0 -0.4017310 0.2799389 0
4 0 1 1 2 1 -0.9562981 -0.6015257 0
5 0 1 0 1 0 0.4368910 0.5277091 1
6 0 1 1 4 2 1.3944855 -0.7075348 0
予測変数として、グループあたりの人数、子供の数、グループがキャンパーで構成されているかどうかを選んだ。多くのグループは、まったく魚を釣ろうとしないかもしれないので（したがって、多くのゼロ回答をもたらす）、ゼロ膨張ポアソン・モデルをデータに適合させる。今のところ、オブザベーション全体にわたって一定のゼロ膨張確率を仮定する。
fit_zinb1 <- brm(count ~ persons + child + camper,
data = zinb, family = zero_inflated_poisson())
ここでも、通常の方法で結果をまとめている。
Family: zero_inflated_poisson
Links: mu = log; zi = identity
Formula: count ~ persons + child + camper
Data: zinb (Number of observations: 250)
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 -1.01 0.18 -1.36 -0.67 1.00 2524 2653
persons 0.87 0.05 0.79 0.97 1.00 2569 2155
child -1.36 0.09 -1.55 -1.19 1.00 2552 2743
camper 0.80 0.09 0.62 0.99 1.00 3200 2715
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi 0.41 0.04 0.32 0.49 1.00 3016 2432
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_zinb1), ask = FALSE)
パラメータ推定値によると、大きなグループはより多くの魚を捕り、キャンプをする人はしない人に比べてより多くの魚を捕り、子供の多いグループはより少ない魚を捕るということである。ゼロインフレ確率
zi
は、平均が41%とかなり大きい。なお、魚を釣らない確率は実際には41%より高いのであるが、この確率の一部はすでにポアソン分布自体でモデル化されている（だからゼロインフレという名前）。もし、すべてのゼロを別のプロセスに由来するものとして扱いたい場合は、代わりにハードルモデル（ここでは示していない）を使用することができる。
ここで、子供の数によってゼロインフレ確率を追加予測しようとする。その理由は、子供の数が多いグループは、魚を捕ろうともしないと予想されるからである。ほとんどの子どもは、何かが起こるまで何時間も待つのがひどく苦手なのである。純粋に統計学的な観点からは、ゼロインフレート（およびハードル）分布は2つのプロセスの混合であり、モデルの両方の部分を予測することは自然であり、データを十分に活用するためにしばしば非常に合理的なことなのである。
fit_zinb2 <- brm(bf(count ~ persons + child + camper, zi ~ child),
data = zinb, family = zero_inflated_poisson())
Family: zero_inflated_poisson
Links: mu = log; zi = logit
Formula: count ~ persons + child + camper
zi ~ child
Data: zinb (Number of observations: 250)
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 -1.08 0.19 -1.46 -0.71 1.00 2731 2450
zi_Intercept -0.96 0.26 -1.48 -0.47 1.00 3832 2991
persons 0.89 0.05 0.80 0.99 1.00 2882 2463
child -1.17 0.10 -1.37 -1.00 1.00 2915 3010
camper 0.78 0.09 0.60 0.96 1.00 3305 2680
zi_child 1.22 0.28 0.68 1.77 1.00 4054 2905
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_zinb2), ask = FALSE)
zi
の線形予測器を確率に変換するために、brms は logit-link
を適用する。
\[logit(zi) =
\log\left(\frac{zi}{1-zi}\right) = \eta_{zi}\]
ロジットリンクは $ [0, 1] $
内の値を取り、実線上の値を返す。したがって、それは確率と線形予測変数の間の推移を可能にする。
このモデルによると、子供と一緒に釣りをしようとすると、全体の釣果が減少するだけでなく（モデルのポアソン部分が示唆するように）、全く釣れない確率も大幅に増加する（ゼロ・インフレーション部分が示唆するように）。
Additive Distributional Models (加法的分布モデル)
これまでの例では、マルチレベルのデータがなかったため、brmsの分布回帰のフレームワークの機能を十分に活用することができなかった。以下に紹介する例では、分布モデルでマルチレベルデータを扱う方法だけでなく、モデルに平滑項（スプラインなど）を組み込む方法も紹介する。多くのアプリケーションにおいて、予測因子と応答の関係がどのようなものか、全く分からない、あるいは非常に曖昧なものでしかない。この問題に取り組むための非常に柔軟なアプローチは、スプラインを使用して、関係の形を把握させることである。説明のために、mgcvパッケージでいくつかのデータをシミュレートする。これはbrmsでも使用され、平滑項を準備する。
dat_smooth <- mgcv::gamSim(eg = 6, n = 200, scale = 2, verbose = FALSE)
Gu & Wahba 4 term additive model
y x0 x1 x2 x3 f
1 10.014562 0.36690008 0.1709383 0.03360030 0.53224372 6.504797
2 7.497204 0.09132215 0.5228089 0.02731544 0.73259184 9.565616
3 21.980127 0.29860651 0.5517069 0.30395678 0.73082890 21.169294
4 18.579703 0.55628934 0.3029421 0.80063964 0.10799495 16.890296
5 8.114506 0.80286967 0.2876062 0.48606145 0.72254170 8.733341
6 11.753057 0.12556843 0.3270885 0.94736718 0.05777094 8.694570
このデータは、予測変数 x0 から x3
と、データの入れ子構造を示すグループ化因子 fac
を含んでいる。私たちは、応答変数 y を、x1 と
x2 の平滑項と fac
の変化する切片を使って予測する。さらに、残差標準偏差 sigma
は、x0 の平滑化項と fac
の変化する切片によって変化すると仮定する。
fit_smooth1 <- brm(
bf(y ~ s(x1) + s(x2) + (1|fac), sigma ~ s(x0) + (1|fac)),
data = dat_smooth, family = gaussian(),
chains = 2, control = list(adapt_delta = 0.95)
)
Family: gaussian
Links: mu = identity; sigma = log
Formula: y ~ s(x1) + s(x2) + (1 | fac)
sigma ~ s(x0) + (1 | fac)
Data: dat_smooth (Number of observations: 200)
Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 2000
Smooth Terms:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sds(sx1_1) 2.68 2.38 0.17 8.99 1.00 711 1091
sds(sx2_1) 18.73 5.49 11.26 32.56 1.00 777 1252
sds(sigma_sx0_1) 1.17 0.95 0.05 3.53 1.00 654 839
Group-Level Effects:
~fac (Number of levels: 4)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 4.73 1.99 2.25 9.62 1.00 1123 1316
sd(sigma_Intercept) 0.17 0.20 0.00 0.68 1.00 599 843
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 15.33 2.25 10.62 19.77 1.00 1075 1218
sigma_Intercept 0.70 0.13 0.39 0.92 1.00 1260 552
sx1_1 16.63 6.71 7.39 34.19 1.01 1122 761
sx2_1 25.66 15.68 -6.59 54.80 1.00 1522 1395
sigma_sx0_1 -1.49 2.34 -7.35 2.17 1.00 1260 1061
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_smooth1), points = TRUE, ask = FALSE)
このモデルは、手元のデータに対して過剰なものである可能性があるが、brmsで複雑なモデルを指定し、バックエンドでStanを使用して適合させることが簡単にできることをうまく示している。