Stan and brms output comparison

Interfaces brms
,
1

Hi, I’m trying to recreate the brms model found here: Bayesian meta-analysis in brms | A. Solomon Kurz in Stan.
Specifically, this one:

b14.5 <- 
  brm(data = spank, family = gaussian,
      d | se(se) ~ 1 + (1 | study),
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(cauchy(0, 1), class = sd)),
      iter = 2000, warmup = 1000, cores = 4, chains = 4,
      seed = 14)

I know that there’s stancode() function in brms, but the output is too difficult for me to parse. I just want the bare-bones structure to start with.

Here’s my Stan code:

data {
  int <lower=0> N; // total number of observations
  vector[N] d; // dependent variable
  vector <lower=0>[N] se; // standard errors
  // data for group-level effects
  int <lower=1> K; // number of groups
  array[N] int<lower=1, upper=K> J; // grouping indicator per observation
}
parameters {
  real<lower=0> tau;
  real Intercept;
  vector[K] theta;
}
model {
  for (n in 1:N) {
      target += normal_lpdf(d[n] | theta[J[n]], se[n]);
      target += normal_lpdf(theta[J[n]] | Intercept, tau);
  }
  target += normal_lpdf(Intercept | 0, 1);
  target += cauchy_lpdf(tau | 0, 1);
}

However, I get slightly different output:
brms output:

Intercept = 0.38 [.30, .45]
sd(Intercept) = 0.26 [.21, .33]

And the output from my model is:
Intercept = 0.33 [.30, .38]
tau = 0.21 [.18, .24]

I think the differences are a bit too large for mcmc simulation variation, so I must be misspecifying the model, but can’t figure what exactly am I doing wrong? Thanks.

1 Like
2

Can you also post the output from the make_stancode call?

3

Here’s the brms generated code:

// generated with brms 2.16.1
functions {
}
data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  vector<lower=0>[N] se;  // known sampling error
  // 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
  int<lower=1> J_1[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_1;
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  vector<lower=0>[N] se2 = square(se);
}
parameters {
  real Intercept;  // temporary intercept for centered predictors
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  vector[N_1] z_1[M_1];  // standardized group-level effects
}
transformed parameters {
  real<lower=0> sigma = 0;  // dispersion parameter
  vector[N_1] r_1_1;  // actual group-level effects
  r_1_1 = (sd_1[1] * (z_1[1]));
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = Intercept + rep_vector(0.0, N);
    for (n in 1:N) {
      // add more terms to the linear predictor
      mu[n] += r_1_1[J_1[n]] * Z_1_1[n];
    }
    target += normal_lpdf(Y | mu, se);
  }
  // priors including constants
  target += normal_lpdf(Intercept | 0, 1);
  target += cauchy_lpdf(sd_1 | 0, 1)
    - 1 * cauchy_lccdf(0 | 0, 1);
  target += std_normal_lpdf(z_1[1]);
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
}