 # 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.

``````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);
}
``````
3 Likes

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 (log) likelihood (of the historical data) and multiply that by `a0`.

Please let me know what/if I misunderstand.

2 Likes

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.