# How to express power prior in stan

Let D0 denote the historical data. I would like to calculate the posterior distribution in the picture.

But I am not sure how stan to express it. If there is no power delta, it is very straightforward.

[Please include Stan program and accompanying data if possible]

You can just take the log of your expression and increment the log probability. The power turns into a factor to multiply the log-likelihood with. For example, if your likelihood function is Gaussian, you can write

``````target += delta*normal_lpdf(D_0|mu, sigma);
``````
2 Likes

You usually want to build a better model rather than hacking the likelihood in this way. You have to make sure you include non-constant normalizing terms. So if your `D_0` is a parameter, you may have some extra work to do.

Iâ€™m very late to this discussion, but in case someone else stumbles upon this thread with the same question (like I just did), hereâ€™s a simple example solution that I cooked up. This uses the â€śeight_schoolsâ€ť data from Gelman et alâ€™s (2003) Bayesian Data Analysis book. In these data, we have have means (y) and standard errors (sigma) summarizing the effect of coaching at eight high schools. In the below, I fit a very simple partial pooling model but there are other (more efficient?) approaches.

To set up for the power prior, I create a grouping variable and then use the power prior to give more weight to group #2. This mimics a situation where perhaps group #1 reflects a historical observation period or a pre-test that, although useful, warrants less weight than the observations in group #2.

``````library(rstan)
library(shinystan)

# the original schools data from Gelman et al (2003):
J <- 8
y <- c(28,  8, -3,  7, -1, 1, 18, 12)
sigma <- c(15, 10, 16, 11, 9, 11, 10, 18)
schools.dat <- list(J=J, y=y, sigma=sigma)

# I now add a hypothetical grouping variable to demo a power prior.
# This variable is NOT in Gelman et al (2003)
# Note that group #1 is associated with schools with lower coaching effects
schools.dat\$period <- c(2,   1,  2,  1,  1, 1,  2, 2)

schools.code5 <- '
data {
int<lower=0> J;          // number of schools
real y[J];               // estimated treatment effect (school j)
real<lower=0> sigma[J];  // std err of effect estimate (school j)
real<lower=0> period[J];  // group number
}
parameters {
real mu;
real theta[J];
real<lower=0> tau;        // variance between schools
}
model {
for (n in 1:J) {
theta[n] ~ normal(mu, tau);
if (period[n] == 1)
target += normal_lpdf(y[n] | theta[n], sigma[n]);
if (period[n] == 2)
target += 2 * normal_lpdf(y[n] | theta[n], sigma[n]); // Note the power prior of 2
}
}
'
schools.mod5 <- stan(model_code=schools.code5, data=schools.dat, chains=4, iter=2000)
launch_shinystan(schools.mod5)
summary(schools.mod5, probs = c(0.1, 0.9))\$summary
``````

â€¦And thatâ€™s it.

4 Likes