Bayesian school effects

This example fits the school effects model, score ~ study_hours + (1 | school), as a Bayesian model in Stan and compares the posterior with the Phit\fit_lmer() estimates. CmdStanPHP builds the Stan data from the same formula through phit, then runs CmdStan.

Compiling and sampling a Stan model takes far longer than a page request, so the page does not run CmdStan. It reads draws from a seeded fit that is stored in the repository. The command php artisan stan:fit-school-posterior regenerates them, and only that command needs CmdStan.

Code

$df = DataFrame::fromQuery(StudentScore::query());
$data = StanData::fromFormula('score ~ study_hours + (1 | school)', $df);
// N, K, X (with the intercept column), y, J, and school_id

$model = cmdstan_model('school_intercepts.stan');
$fit = $model->sample(data: $data, seed: 20261016, chains: 4, thin: 4);
$draws = $fit->drawsColumns(['beta', 'sigma', 'tau', 'school_intercept']);

// score ~ study_hours + (1 | school), with a non-centered random intercept.
// Data come from CmdStan\StanData::fromFormula(); X includes the intercept column.
data {
  int<lower=0> N;
  int<lower=1> K;
  matrix[N, K] X;
  vector[N] y;
  int<lower=1> J;
  array[N] int<lower=1, upper=J> school_id;
}
parameters {
  vector[K] beta;
  real<lower=0> sigma;
  real<lower=0> tau;
  vector[J] z;
}
transformed parameters {
  vector[J] u = tau * z;
}
model {
  beta[1] ~ normal(50, 25);
  beta[2:K] ~ normal(0, 10);
  sigma ~ student_t(3, 0, 10);
  tau ~ student_t(3, 0, 10);
  z ~ std_normal();
  y ~ normal(X * beta + u[school_id], sigma);
}
generated quantities {
  vector[J] school_intercept = beta[1] + u;
}

The priors are weakly informative: normal(50, 25) on the intercept, normal(0, 10) on the slope, and half-Student-t(3, 0, 10) on σ and τ. The random intercepts are non-centered.