Integrating previous model's parameters as priors for Bayesian modeling of new data

Question

I'm working on a project for scientific publication (in the social sciences) in which I am structuring my manuscript with a Study 1, Study 2, and so on format. Basically, I have a compelling finding with one dataset which itself in my field is probably sufficient for a publication. However, I know that my colleagues would see this as more interesting if it were replicated in a separate dataset, so that is what I am doing.

So let's make this a more concrete example (this is generated data, but poses a statistically similar problem to my real one). I am predicting that the average number of donuts a person consumes on a weekly basis is associated with greater body weight. For the purposes of the question, let's just assume this effect is independent of all other demographics that we might normally control for in a more complex model.

Here is the data generating process:

set.seed(3)
donuts1 <- runif(100, 0, 15)
weight1 <- 5 * donuts1 + rnorm(100, 0, 50) + 185

As expected, running lm(weight1 ~ donuts1) reveals a statistically significant effect of donut consumption on weight.

Great! Now I will go to another population and collect more data. While the effect is constant in this new population, there is more random error.

set.seed(8)
donuts2 <- runif(100, 0, 15)
weight2 <- 5 * donuts2 + rnorm(100, 0, 100) + 185

Now if I report my Study 2 as lm(weight2 ~ donuts2), we will see a p-value of about .22 and therefore a failure to replicate. Now, a reasonable person could eyeball the Study 2 data and say that despite statistical insignificance, Study 2's data is consistent with the data in Study 1. Still, eyeballing it and saying "close enough" isn't usually what we're going for when it comes to statistical inference.

The way I see it, there are a few options:

• Combine the two samples into one and run analyses on the combined data. This would be seen as odd in my field and many others.
• Treat them like I would a meta-analysis, analyzing the effect sizes and weighting by precision.
• Incorporate Study 1's information as an informative prior for a Bayesian approach to analyzing Study 2.

I am particularly interested in that last option even if for purposes of demonstration only—I am learning Bayesian analysis and want to see if I understand things correctly.

So here's the question: How do I incorporate the results of Study 1 as an informative prior for the analysis of Study 2?

My inclination would be to do this, demonstrated below with code to run an analysis in JAGS.

library(rjags)
library(runjags)

model.inform <- '
model {
  for (i in 1:N) {
    weight2[i] ~ dnorm(y.hat[i], tau)
    y.hat[i] <- a + b * donuts2[i]
  }

  a ~ dnorm(183.41, 9.11)
  b ~ dnorm(5.07, 1.08)

  tau <- pow(sigma, -2)
  sigma ~ dunif(0, 100)
}
'

I've set the priors on the a and b terms (assuming a linear model: $Y = a + bx$) to the parameter estimates from Study 1. The prior's standard deviation is set to the OLS-derived standard error. If you run the analysis (code at end of the question), you'll see that this better reflects what we know from the two datasets in that the 95% credibility interval excludes zero with a median estimate close to the true effect. I'll note that using a flat prior results in a lower estimated effect with a much wider confidence interval.

With that said, I've never seen anything like this in published research, though I confess I rarely see any Bayesian analysis in my field. This is similar to Empirical Bayes from what I understand, but in this case I'm not using the same data for estimating both prior and posterior. It strikes me that perhaps what I've done does not express enough uncertainty about the priors, but I'm not sure.

I understand that many call these subjective priors for a reason, but I think this question is itself not so subjective that it can't be answered on Cross-Validated even if there may be multiple ways to approach this problem and to model the same general approach.

Code to fit the model above in JAGS:

fit.inform <- run.jags(model.inform, c('a','b'),
                data = list('donuts2' = donuts2,
                            'weight2' = weight2,
                            'N' = 100),
                inits=list(list(.RNG.name="lecuyer::RngStream"),
                           list(.RNG.name="lecuyer::RngStream"),
                           list(.RNG.name="lecuyer::RngStream")),
                burnin = 1000, sample = 5000, thin = 15, summarise = TRUE,
                n.chains = 3)

Answer 1

In general, informing a prior requires a lot of judgment calls (and justification in the write up). There are several steps:

1. Collect the relevant previous studies that could inform the present one. This step is much like collecting previous studies for meta-analysis. You want to be sure it doesn't suffer from the file-drawer problem, whereby you'd only use published research that happened to reach significance and not the unpublished studies that happened not reach significance. In your case, this step might be easy because you've already restricted yourself to considering one specific previous study. (But then be careful not to dismiss "inconvenient" results of previous studies.)

2. Decide how to use the results of the previous studies to inform the present study. In general, previous studies might use different designs, different variables, different analyses. It's an artful translation from the parameters of different previous models to priors on the parameters in a new model. In your case, this step might be easy if your second study is the same structure as the previous study --- the mapping of parameters from previous analysis to new analysis is direct.

3. If you have an MCMC-sampled posterior distribution of the previous study, then figure out how to inform the mathematically-specified prior of the new study. In general, this step is not trivial because the posterior distribution of the previous study is, presumably, a complex multi-dimensional distribution almost certainly with correlations of parameters, and probably with skewed and/or kurtotic tails. You need to decide what sort of mathematical distribution can reasonably approximate the MCMC sample. Then you need to put that mathematical expression in the prior of JAGS. In your case, since you're doing linear regression, it could be important to include the correlations of the parameters.

4. If the structure of the second study is exactly the same as the first study, then you are, essentially, just collecting more data in the same study, and you can just collapse the two data sets into one. But presumably there is something different about the second study, so you can't just use the posterior of the first study directly and exactly for the second study, because if you did that's tantamount to saying you're just collecting more data in the same study. Instead, to capture the fact that there's something different about the second study, you should "relax" the prior, that is, make it more uncertain, because you're not sure how much the prior research really should apply to the second study. How much should you relax the prior? Well, remember the first sentence about a lot of judgment calls?

Other respondents may, of course, have more specific recommendations and might even point to examples in the literature that use informed priors. I'd just say to keep the above in mind as you read those examples.