Power Prior for a Hierarchical Model in Stan

I am attempting to implement a Power prior for mixed effects models, where both the current and historical data are modelled using hierarchical models.The method is outlined in section 4 of https://projecteuclid.org/euclid.ss/1009212673, - it involves specifying a power prior which is the integral of the historical data likelihood (multiplied by the density of the random effects) with respect to the random effects parameters. This part I can’t figure out how to implement in Stan. The code below implements the power prior but does not correctly account for the integration of the random effect from the historical likelihood.

  int<lower=1> N; //Number of pupils - current
  int<lower=1> J; //Number of schools (current)
  int<lower=1> N0; // Number of pupils - historic
  int<lower = 1> J0; //number of schools (historic)
  real y[N]; //outcome - current
  real y0[N0]; //outcome - historic
  int<lower=0,upper=J> SchoolCode[N];
  int<lower=0,upper=J0> SchoolCode0[N0];
  vector[J] school_data; // current school-level predictors
  vector[J0] school_data0; //historic school level predictors
  real<lower=0,upper=1> a0; //discounting parameter

  real<lower=0> sigma;
  real alpha;
  //varying intercepts - current
  vector[J] eta_raw; //N(0,1) params of school-specific intercepts
  real<lower = 0> sigma_eta;
  real zeta; // school-level coefficients
  //varying intercepts - historical
  vector[J0] eta0_raw;

transformed parameters{
  vector[J] eta = alpha + school_data * zeta + sigma_eta * eta_raw;
  vector[J0] eta0 = alpha + school_data0 * zeta + sigma_eta * eta0_raw;

  vector[N] mu = eta[SchoolCode];
  vector[N0] mu0 = eta0[SchoolCode0];
  target += a0 * normal_lpdf(y0|mu0, sigma);
  eta_raw ~ normal(0,1);
  eta0_raw ~ normal(0,1);
  alpha ~ normal(0,5);
  zeta ~ normal(0,15);
  sigma ~ normal(0, 10);
  sigma_eta ~ exponential(1);
  y ~ normal(mu, sigma);

Sorry it took so long to get to your question. I’m afraid I don’t follow what

means. If a_0 is fixed - as opposed to being assigned a prior measure and allowed to vary - there’s no need for integration. Your program actually looks fine to me. As a general approach, implementing a fixed a_0 power prior entails only finding the line corresponding to the likelihood (of the historical data) and multiply that by a0.

Please let me know what/if I misunderstand.


Thanks a lot for taking a look at this.

I think I agree now that what I did initially was ok. What had confused me was the integral in equation (4.3) in the paper I have linked above, as the random effects parameters in the power prior are different to those in the current likelihood and have therefore been integrated out.