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.

```
data{
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
}
parameters{
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;
}
model{
vector[N] mu = eta[SchoolCode];
vector[N0] mu0 = eta0[SchoolCode0];
//priors
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);
//likelihood
y ~ normal(mu, sigma);
}
```