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

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)

• If you want to perform a Bayesian analysis, then, since the model (the likelihood) is the same in the two studies, the best way would be to use Bayesian updating, i.e., 1) choose prior, 2) compute posterior from Study 1, 2) use posterior from Study 1 as prior for Study 2, 3) compute posterior. Note that this will result in exactly the same posterior as if you would have combined the two studies together. – DeltaIV Apr 4 '17 at 20:02