JMM’s notes on

BRMS (Bayesian Regression Models using Stan)

I typically just use plain Stan (see /stan for notes there), but I’ll try to document some BRMS stuff since other people in the lab use it.

Nix setup

Here’s the default.nix I use to set up BRMS. 

{ pkgs ? import <nixpkgs> {}
, ... }:

with pkgs;
{
  myenv = let
    my-R-packages = with rPackages; [
      R
      tidyverse
      jsonlite
      data_table
      rstan
      ggplot2
      brms
    ];
    newR = rWrapper.override {
      packages = my-R-packages;
    };
  in
    buildEnv {
      name = "myenv";
      paths = [
        newR
        gcc                     # For rstan
        gnumake                 # Ditto
      ];
    };
}

Then you can have a direnv (see also my direnv notes) like so: 

#!/usr/bin/env bash
# .envrc files are run in bash,
# They just set any variables that are exported.

if [ ! -d "myenv" ]; then
    echo "You should build the myenv directory"
    echo "Probably using: nix build -vL -f . myenv -o myenv"
else
    :
fi

PATH_add myenv/bin

# Local Variables:
# mode: sh
# sh-shell: bash
# End:

Intro example

I think this is from some intro to BRMS. Can’t remember where I got it from, but I used this back in 2022. 

> library(brms)
> library(rstan)
Loading required package: StanHeaders
Loading required package: ggplot2
rstan (Version 2.21.3, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
> rstan_options(auto_write = TRUE)
> options(mc.cores = parallel::detectCores())
> 
+ fitmodel <- brm(count ~ zAge + zBase * Trt + (1|patient),
+                 family = poisson())
+ 
+ Error: Data must be specified using the 'data' argument.
> 
+ fit1 <- brm(count ~ zAge + zBase * Trt + (1|patient),
+             data = epilepsy, family = poisson())
+ 
+ Compiling Stan program...
Start sampling

SAMPLING FOR MODEL '04f7ce655d7a14cc1ac82ece51dcda88' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00026 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.6 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '04f7ce655d7a14cc1ac82ece51dcda88' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000164 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.64 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '04f7ce655d7a14cc1ac82ece51dcda88' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.00035 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.5 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '04f7ce655d7a14cc1ac82ece51dcda88' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000152 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.52 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.99756 seconds (Warm-up)
Chain 1:                4.24329 seconds (Sampling)
Chain 1:                10.2408 seconds (Total)
Chain 1: 
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 5.47962 seconds (Warm-up)
Chain 4:                4.57296 seconds (Sampling)
Chain 4:                10.0526 seconds (Total)
Chain 4: 
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 6.13173 seconds (Warm-up)
Chain 3:                4.47731 seconds (Sampling)
Chain 3:                10.609 seconds (Total)
Chain 3: 
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 6.14743 seconds (Warm-up)
Chain 2:                5.10291 seconds (Sampling)
Chain 2:                11.2503 seconds (Total)
Chain 2: 
> 
+ plot(fit1, variable = c("b_Trt1", "b_zBase"))
+ 
> >

I guess I couldn’t just compile a model without data?

Generating Stan code

This is probably also from the same intro. Just found in an old file on my hard drive: 

scode <- make_stancode(count ~ zAge + zBase * Trt + (1|patient),
                       data = epilepsy, family = poisson())

sdata <- make_standata(count ~ zAge + zBase * Trt + (1|patient),
                       data = epilepsy, family = poisson())

stanmodel <- stan_model(model_code = scode)

stanfit <- sampling(stanmodel, data = sdata)