BRMS to STAN Compatibility

Modeling
,
1

Greetings, I was using the following code to model data and decided to use the make_stan code:

ClassicGP_ROI0 <- brm(RT ~ AMBUAMB * (SZM1 + SZM2) +
                        (1 + AMBUAMB * (SZM1 + SZM2) || MD5) +
                        (1 + AMBUAMB * (SZM1 + SZM2) || item),
                      data=ClassicGP_test[ClassicGP_test$ROI==0,],
                      prior = brms_parms$prior,
                      cores = brms_parms$ncores,
                      iter = brms_parms$niters,
                      seed = brms_parms$seed,
                      warmup = brms_parms$warmup,
                      control = list(adapt_delta = brms_parms$adapt_delta))
stancode(ClassicGP_ROI0,save_model="brms_to_stan")

And it gave me the following code:

functions {
}
data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int<lower=1> Kc;  // number of population-level effects after centering
  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  array[N] int<lower=1> J_1;  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_1;
  vector[N] Z_1_2;
  vector[N] Z_1_3;
  vector[N] Z_1_4;
  vector[N] Z_1_5;
  vector[N] Z_1_6;
  // data for group-level effects of ID 2
  int<lower=1> N_2;  // number of grouping levels
  int<lower=1> M_2;  // number of coefficients per level
  array[N] int<lower=1> J_2;  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_2_1;
  vector[N] Z_2_2;
  vector[N] Z_2_3;
  vector[N] Z_2_4;
  vector[N] Z_2_5;
  vector[N] Z_2_6;
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  }
}
parameters {
  vector[Kc] b;  // regression coefficients
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> sigma;  // dispersion parameter
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  array[M_1] vector[N_1] z_1;  // standardized group-level effects
  vector<lower=0>[M_2] sd_2;  // group-level standard deviations
  array[M_2] vector[N_2] z_2;  // standardized group-level effects
}
transformed parameters {
  vector[N_1] r_1_1;  // actual group-level effects
  vector[N_1] r_1_2;  // actual group-level effects
  vector[N_1] r_1_3;  // actual group-level effects
  vector[N_1] r_1_4;  // actual group-level effects
  vector[N_1] r_1_5;  // actual group-level effects
  vector[N_1] r_1_6;  // actual group-level effects
  vector[N_2] r_2_1;  // actual group-level effects
  vector[N_2] r_2_2;  // actual group-level effects
  vector[N_2] r_2_3;  // actual group-level effects
  vector[N_2] r_2_4;  // actual group-level effects
  vector[N_2] r_2_5;  // actual group-level effects
  vector[N_2] r_2_6;  // actual group-level effects
  real lprior = 0;  // prior contributions to the log posterior
  r_1_1 = (sd_1[1] * (z_1[1]));
  r_1_2 = (sd_1[2] * (z_1[2]));
  r_1_3 = (sd_1[3] * (z_1[3]));
  r_1_4 = (sd_1[4] * (z_1[4]));
  r_1_5 = (sd_1[5] * (z_1[5]));
  r_1_6 = (sd_1[6] * (z_1[6]));
  r_2_1 = (sd_2[1] * (z_2[1]));
  r_2_2 = (sd_2[2] * (z_2[2]));
  r_2_3 = (sd_2[3] * (z_2[3]));
  r_2_4 = (sd_2[4] * (z_2[4]));
  r_2_5 = (sd_2[5] * (z_2[5]));
  r_2_6 = (sd_2[6] * (z_2[6]));
  lprior += normal_lpdf(b | 0,150);
  lprior += normal_lpdf(Intercept | 300,1000);
  lprior += normal_lpdf(sigma | 0,500)
    - 1 * normal_lccdf(0 | 0,500);
  lprior += normal_lpdf(sd_1 | 0,200)
    - 6 * normal_lccdf(0 | 0,200);
  lprior += normal_lpdf(sd_2 | 0,200)
    - 6 * normal_lccdf(0 | 0,200);
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = rep_vector(0.0, N);
    mu += Intercept;
    for (n in 1:N) {
      // add more terms to the linear predictor
      mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_1_2[J_1[n]] * Z_1_2[n] + r_1_3[J_1[n]] * Z_1_3[n] + r_1_4[J_1[n]] * Z_1_4[n] + r_1_5[J_1[n]] * Z_1_5[n] + r_1_6[J_1[n]] * Z_1_6[n] + r_2_1[J_2[n]] * Z_2_1[n] + r_2_2[J_2[n]] * Z_2_2[n] + r_2_3[J_2[n]] * Z_2_3[n] + r_2_4[J_2[n]] * Z_2_4[n] + r_2_5[J_2[n]] * Z_2_5[n] + r_2_6[J_2[n]] * Z_2_6[n];
    }
    target += normal_id_glm_lpdf(Y | Xc, mu, b, sigma);
  }
  // priors including constants
  target += lprior;
  target += std_normal_lpdf(z_1[1]);
  target += std_normal_lpdf(z_1[2]);
  target += std_normal_lpdf(z_1[3]);
  target += std_normal_lpdf(z_1[4]);
  target += std_normal_lpdf(z_1[5]);
  target += std_normal_lpdf(z_1[6]);
  target += std_normal_lpdf(z_2[1]);
  target += std_normal_lpdf(z_2[2]);
  target += std_normal_lpdf(z_2[3]);
  target += std_normal_lpdf(z_2[4]);
  target += std_normal_lpdf(z_2[5]);
  target += std_normal_lpdf(z_2[6]);
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
}

This is no where near compatible with my data formats in brms, and I wonder if there is a way for stancode to align with the data in brms?

Thank you for any help provided!

2

Under the hood, brms formats your data for use with this Stan program using brms::make_standata. You can do the same thing explicitly if you’d like to run this Stan program yourself, e.g. using rstan or cmdstanr.

2 Likes
3

Thank you so much! So does this imply that if I also were to run make_standata on my own data, it would produce data compatible to the stan code generated?

4

Yes :)

Note that make_standata formats the data for a particular model under consideration, when you say “run make_standata on my own data” what you really mean is to run it on your own data and model specification.

2 Likes