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);
}