公式VignetteのHow do I install brms?に記載されているように，brmsはCRANから容易にインストール可能です． ただしその前に，使用予定のRのバージョンに応じて，適切なRtoolsをインストールしておいてください． URLも，上記のVignetteに書いてある通りです．
※例えば，Rバージョン4.4.2を使用するつもりなら，Rtools4.4をインストールする，という具合です．
fit <-brm(mpg ~ wt * am,
data = mtcars,
iter =2000,
warmup =1000,
seed =1234,
chain =4,
cores =4# 使用するPCの性能に応じて調整)
サンプリング結果を格納したオブジェクトをsummary()に入れると要約された結果が返されます．
summary(fit)# 以下，出力 ---------------------------------
Family: gaussian
Links: mu = identity
Formula: mpg ~ wt * am
Data:mtcars (Number of observations:32)
Draws:4 chains, each with iter =2000; warmup =1000; thin =1;
total post-warmup draws =4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 31.553.1925.2337.831.0018852339
wt -3.820.83-5.45-2.191.0019232322
am1 14.674.495.9723.541.0016191935
wt:am1 -5.231.50-8.24-2.281.0017072229
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 2.680.362.083.511.0025682498
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).
brms::waic(fit)
Computed from 4000 by 32 log-likelihood matrix.
Estimate SE
elpd_waic -78.84.0
p_waic 4.11.0
waic 157.67.94 (12.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.
brms::loo(fit)
Computed from 4000 by 32 log-likelihood matrix.
Estimate SE
elpd_loo -79.04.0
p_loo 4.31.1
looic 158.08.0------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in[0.5,0.9]).
All Pareto k estimates are good (k <0.7).
See help('pareto-k-diagnostic')for details.
data {
int<lower=0> N;// number of data items
int<lower=0> K;// number of predictors
matrix[N, K] x;// predictor matrixvector[N] y;// outcome vector}
parameters {
real alpha;// intercept
vector[K] beta;// coefficients for predictors
real<lower=0> sigma;// error scale
}
model {for(n in1:N){
y[n]~normal(x[n]* beta, sigma);}}
brms::hypothesis(fit,"wt = 0")# 以下，出力 -----------------------
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (wt)=0-3.820.83-5.45-2.19NANA*---'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.
# 無情報でない事前分布 -----------------
priors <-c(set_prior("normal(0,5)", class ="b"),set_prior("normal(0,5)", class ="Intercept"),set_prior("normal(0,5)", class ="sigma"))
fit <-brm(mpg ~ wt * am,
data = mtcars,
iter =20000,# サンプリング回数を大きくした
warmup =2000,
seed =1234,
chain =4,
cores =4,
save_pars =save_pars(all =TRUE),
sample_prior ="yes",# ここが重要！！
prior = priors # 無情報でない事前分布の指定)
brms::hypothesis(fit,"wt = 0")# 以下，出力 -----------------------
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (wt)=0-4.590.75-6.09-3.1400*---'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.
brms::hypothesis(fit,"wt + wt:am1 = 0")# 以下，出力 -----------------------
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (wt+wt:am1)=0-7.591.16-9.8-5.2200*---'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.
stats::glm()と同様に，目的変数が従う確率分布を，引数family =に指定しましょう．今回は，オートマ車かマニュアル車かamは2値の変数なので，family = "bernoulli"とベルヌーイ分布を指定します（余談ですが，stats::glm()では2値の目的変数に対しては二項分布（family = "binomial"）を指定しますが，brmsでは2値の目的変数に対して二項分布を指定すると，Only 2 levels detected so that family 'bernoulli' might be a more efficient choice.とサジェストされます）．
#amをfactor型に変換していたのを，いったんリセットdata(mtcars)
fit <-glm(cbind(am,1- am)~ mpg + wt,
family ="binomial",
data = mtcars)#目的変数が2値の場合は、#fit <- glm(am ~ mpg + wt, family = "binomial", data = mtcars) でOKsummary(fit)
#amをfactor型に変換していたのを，いったんリセットdata(mtcars)
fit <-brm(am |trials(1)~ mpg + wt,
family ="binomial",
data = mtcars)#目的変数が2値の場合は、#fit <- brm(am ~ mpg + wt, family = "bernoulli", data = mtcars) でOKsummary(fit)
fit_zip <-brm(num ~ minute,
data =data.frame(num, minute),
family ="zero_inflated_poisson",
seed =1234)summary(fit_zip)# 以下，出力 ----------------------------
Family: zero_inflated_poisson
Links: mu = log
Formula: num ~ minute
Data:data.frame(num, minute)(Number of observations:300)
Draws:4 chains, each with iter =2000; warmup =1000; thin =1;
total post-warmup draws =4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.010.193.664.361.0031102872
minute -0.080.01-0.09-0.071.0025762570
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi 0.300.030.230.361.0023612280
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).
fit_zip2 <-brm(brms::bf(num ~ minute, zi ~ minute),
data =data.frame(num, minute),
family ="zero_inflated_poisson",
seed =1234)summary(fit_zip2)# 以下，出力 -----------------------
Family: zero_inflated_poisson
Links: mu = log; zi = logit
Formula: num ~ minute
zi ~ minute
Data:data.frame(num, minute)(Number of observations:300)
Draws:4 chains, each with iter =2000; warmup =1000; thin =1;
total post-warmup draws =4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.050.193.674.411.0026003025
zi_Intercept -0.410.81-1.971.221.0023972665
minute -0.080.01-0.09-0.071.0022192575
zi_minute -0.010.02-0.060.031.0020862205
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).
lme4::glmer(y ~ x + (1 | group), family = "poisson", data = dat)
brmsでは，関数をbrm()に変えるだけで同じように階層線形モデルを推定できます．
ただし，ベイズモデリングの強みは柔軟なモデルを表現・推定できることですから，最尤法では推定できなかったような複雑なモデルも表現できます．やや高度な話になるので詳述はしませんが，詳しくは以下のリンクからAdvanced Multilevel Modeling with brmsを参照してください．このページにはそのほかにも高度な書き方を要する情報が掲載されています．
話を簡単にするために，ここでは目的変数に「病気の再発までの時間time」，説明変数に「患者の年齢age」を投入した，単回帰分析を考えます．また，本当は右側打ち切り（right-censored）が起こっているため，それを考慮したモデリングをするべきですが（brmsなら容易に実現できます），打ち切りが起こっていないデータだけを抽出しています．
なおbrmsでは，brm(y | cens(censored) ~ x, data = dat)のように書くことで，打ち切りデータのモデリングも容易に実現できます（cens()のなかに，打ち切りの有無を判別する変数を投入します）．
#データセットの作成 ---------------------------------------------#1回目の再発までの時間 & 非打ち切りデータのみを使用
dat <- brms::kidney %>%
dplyr::filter(recur =="1"& censored ==0)#データセットの行列数 -------------------------------------------dim(dat)#38人分のデータが抽出された[1]387#冒頭6行分を表示 ------------------------------------------------head(dat)
time censored patient recur age sex disease
1801128 male other
22302148 female GN
32203132 male other
444704131 female other
53005110 male other
62406116 female other
fit_q <-brm(bf(mpg ~1, quantile =0.2),
data = mtcars,
family =asym_laplace(),
seed =1234)summary(fit_q)# 以下，出力 ----------------------------
Family: asym_laplace
Links: mu = identity
Formula: mpg ~1
quantile =0.2
Data:mtcars (Number of observations:32)
Draws:4 chains, each with iter =2000; warmup =1000; thin =1;
total post-warmup draws =4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 15.130.6613.7916.451.0025731785
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.510.281.052.151.0028702695
quantile 0.200.000.200.20NANANA
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).
