brmspy - google colab¶
In [1]:
Copied!
import sys, os # for running from repo
sys.path.insert(0, os.path.abspath("../../"))
try: from brmspy import brms
except ImportError:
%pip install -q brmspy
from brmspy import brms
import sys, os # for running from repo
sys.path.insert(0, os.path.abspath("../../"))
try: from brmspy import brms
except ImportError:
%pip install -q brmspy
from brmspy import brms
[brmspy][_call_with_frames_removed] Running autoload!
In [2]:
Copied!
with brms.manage() as ctx:
ctx.install_brms(use_prebuilt=True)
with brms.manage() as ctx:
ctx.install_brms(use_prebuilt=True)
[brmspy][install_brms] Activating runtime /Users/sebastian/.brmspy/runtime/macos-arm64-r4.5-0.2.0 [brmspy][install_brms] lib paths are ['/Users/sebastian/.brmspy/environment/default/Rlib', '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library', '/Users/sebastian/.brmspy/runtime/macos-arm64-r4.5-0.2.0/Rlib'] [brmspy][install_brms] Setting cmdstan path to /Users/sebastian/.brmspy/runtime/macos-arm64-r4.5-0.2.0/cmdstan
In [3]:
Copied!
epilepsy = brms.get_brms_data("epilepsy")
formula = "count ~ zAge + zBase * Trt + (1|patient)"
family = "poisson"
brms.default_prior(formula, data=epilepsy, family=family)
epilepsy = brms.get_brms_data("epilepsy")
formula = "count ~ zAge + zBase * Trt + (1|patient)"
family = "poisson"
brms.default_prior(formula, data=epilepsy, family=family)
Out[3]:
| prior | class | coef | group | resp | dpar | nlpar | lb | ub | tag | source | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b | default | |||||||||
| 2 | b | Trt1 | default | ||||||||
| 3 | b | zAge | default | ||||||||
| 4 | b | zBase | default | ||||||||
| 5 | b | zBase:Trt1 | default | ||||||||
| 6 | student_t(3, 1.4, 2.5) | Intercept | default | ||||||||
| 7 | student_t(3, 0, 2.5) | sd | 0 | default | |||||||
| 8 | sd | patient | default | ||||||||
| 9 | sd | Intercept | patient | default |
In [4]:
Copied!
model = brms.brm(
formula=formula,
data=epilepsy,
family=family,
warmup=500,
iter=1000,
chains=4
)
idata = model.idata
model = brms.brm(
formula=formula,
data=epilepsy,
family=family,
warmup=500,
iter=1000,
chains=4
)
idata = model.idata
[brmspy][worker_main] Fitting model with brms (backend: cmdstanr)... Running MCMC with 4 chains, at most 2 in parallel... Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup) Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling) Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling) Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) Chain 1 finished in 1.8 seconds. Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup) Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) Chain 2 finished in 2.0 seconds. Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup) Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup) Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup) Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup) Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup) Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup) Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup) Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup) Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup) Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup) Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling) Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup) Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling) Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling) Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling) Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling) Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling) Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling) Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling) Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling) Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling) Chain 3 Iteration: 1000 / 1000 [100%] (Sampling) Chain 4 Iteration: 1000 / 1000 [100%] (Sampling) Chain 3 finished in 2.9 seconds. Chain 4 finished in 2.6 seconds. All 4 chains finished successfully. Mean chain execution time: 2.3 seconds. Total execution time: 5.0 seconds. [brmspy][worker_main] Fit done!
In [5]:
Copied!
brms.save_rds(model, "epilepsy_fixed_effects.rds")
# load later using: model = brms.load_rds_fit("epilepsy_fixed_effects.rds")
brms.save_rds(model, "epilepsy_fixed_effects.rds")
# load later using: model = brms.load_rds_fit("epilepsy_fixed_effects.rds")
In [6]:
Copied!
epilepsy
epilepsy
Out[6]:
| Age | Base | Trt | patient | visit | count | obs | zAge | zBase | |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 31.0 | 11.0 | 0 | 1 | 1 | 5.0 | 1 | 0.424995 | -0.757173 |
| 2 | 30.0 | 11.0 | 0 | 2 | 1 | 3.0 | 2 | 0.265284 | -0.757173 |
| 3 | 25.0 | 6.0 | 0 | 3 | 1 | 2.0 | 3 | -0.533274 | -0.944403 |
| 4 | 36.0 | 8.0 | 0 | 4 | 1 | 4.0 | 4 | 1.223553 | -0.869511 |
| 5 | 22.0 | 66.0 | 0 | 5 | 1 | 7.0 | 5 | -1.012408 | 1.302363 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 232 | 32.0 | 16.0 | 1 | 55 | 4 | 3.0 | 232 | 0.584707 | -0.569942 |
| 233 | 26.0 | 22.0 | 1 | 56 | 4 | 8.0 | 233 | -0.373562 | -0.345266 |
| 234 | 21.0 | 25.0 | 1 | 57 | 4 | 1.0 | 234 | -1.172120 | -0.232927 |
| 235 | 36.0 | 13.0 | 1 | 58 | 4 | 0.0 | 235 | 1.223553 | -0.682281 |
| 236 | 37.0 | 12.0 | 1 | 59 | 4 | 2.0 | 236 | 1.383264 | -0.719727 |
236 rows × 9 columns
In [7]:
Copied!
c = brms.call("as.data.frame", epilepsy)
c
c = brms.call("as.data.frame", epilepsy)
c
Out[7]:
| Age | Base | Trt | patient | visit | count | obs | zAge | zBase | |
|---|---|---|---|---|---|---|---|---|---|
| _obs_id_ | |||||||||
| 1 | 31.0 | 11.0 | 0 | 1 | 1 | 5.0 | 1 | 0.424995 | -0.757173 |
| 2 | 30.0 | 11.0 | 0 | 2 | 1 | 3.0 | 2 | 0.265284 | -0.757173 |
| 3 | 25.0 | 6.0 | 0 | 3 | 1 | 2.0 | 3 | -0.533274 | -0.944403 |
| 4 | 36.0 | 8.0 | 0 | 4 | 1 | 4.0 | 4 | 1.223553 | -0.869511 |
| 5 | 22.0 | 66.0 | 0 | 5 | 1 | 7.0 | 5 | -1.012408 | 1.302363 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 232 | 32.0 | 16.0 | 1 | 55 | 4 | 3.0 | 232 | 0.584707 | -0.569942 |
| 233 | 26.0 | 22.0 | 1 | 56 | 4 | 8.0 | 233 | -0.373562 | -0.345266 |
| 234 | 21.0 | 25.0 | 1 | 57 | 4 | 1.0 | 234 | -1.172120 | -0.232927 |
| 235 | 36.0 | 13.0 | 1 | 58 | 4 | 0.0 | 235 | 1.223553 | -0.682281 |
| 236 | 37.0 | 12.0 | 1 | 59 | 4 | 2.0 | 236 | 1.383264 | -0.719727 |
236 rows × 9 columns
In [ ]:
Copied!
from arviz_stats import summary
summary = summary(
idata,
kind="stats",
round_to=3
)
print("Posterior Summary")
print("="*60)
print(summary)
from arviz_stats import summary
summary = summary(
idata,
kind="stats",
round_to=3
)
print("Posterior Summary")
print("="*60)
print(summary)
Posterior Summary
============================================================
mean sd eti89_lb eti89_ub
b_Intercept 1.784 0.123 1.592 1.981
b_zAge 0.087 0.088 -0.050 0.229
b_zBase 0.695 0.115 0.516 0.884
b_Trt1 -0.296 0.167 -0.571 -0.033
b_zBase:Trt1 0.061 0.157 -0.191 0.312
... ... ... ... ...
r_patient[57,Intercept] -0.583 0.321 -1.115 -0.079
r_patient[58,Intercept] -1.217 0.423 -1.915 -0.556
r_patient[59,Intercept] -0.141 0.311 -0.655 0.328
lprior -3.183 0.009 -3.198 -3.170
lp__ -702.200 7.799 -714.586 -690.051
[68 rows x 4 columns]
In [15]:
Copied!
import matplotlib.pyplot as plt
from arviz_plots import plot_dist
import seaborn as sns
fig = plot_dist(
idata,
var_names=['b_Intercept', 'b_zAge', 'b_zBase', 'b_Trt1', 'b_zBase:Trt1']
)
plt.suptitle('Posterior Distributions - Fixed Effects', y=1.02, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
import matplotlib.pyplot as plt
from arviz_plots import plot_dist
import seaborn as sns
fig = plot_dist(
idata,
var_names=['b_Intercept', 'b_zAge', 'b_zBase', 'b_Trt1', 'b_zBase:Trt1']
)
plt.suptitle('Posterior Distributions - Fixed Effects', y=1.02, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
In [ ]:
Copied!