How to express power prior in stan

Modeling
1

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

Screen Shot 2017-06-15 at 10.50.54 AM.png1016Ă—106 5.59 KB

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]

2

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
3

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.

4

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