What methods can be used to specify priors from data?

Question

Background

I am generally interested in learning appropriate methods of using data to specify priors. A previous question asks how to elicit priors from experts and received some good recommendations. Here, I am interested in learning how to specify a prior using data. I plan to use these priors in a meta-analysis to synthesize additional data that I collect.

update Although John provides a 'correct' answer, in my case, it would require substantial modification of the original model to implement, so I would prefer to find a way to estimate the prior as a discrete step.

Questions

1. What is the best way to specify such a data informed prior?

2. If I am working with parameters for a particular species (monkeys), and this species belongs to a group of organisms(primates), and data are available for primates but not for monkeys, would it be appropriate to fit a distribution based on the primate data?

Example cases, first with proposed solution

1. I have 100 observations from 100 primate species of primate thumb length:

   set.seed(0)
   thumb <- rgamma(100, 4, 0.1)
   library(MASS)
   fitdistr(thumb, 'gamma')
   

   Indeed, when there is no apriori reason to select a particular distribution, the distribution can be chosen by maximum likelihood:

   for(dist in c('gamma', 'lognormal', 'weibull') {
   logLik(fitdistr(thumb, dist))
   }
   

2. I have collected 50 means, standard errors, and sample sizes from 50 different primate species, and 50 independent observations from another 50 species of eye diameter:

   eye <- data.frame( diameter = rgamma(100, 4, 0.1),
                      se       = c(rlnorm(50, 0.5,1), rep(NA, 50)),
                      n        = c(rep(1:5, 10), rep(1, 50)))      
   eye <- signif(eye, 3)
   

   How can I incorporate the sample statistics into my calculation of a prior?

Answer 1

If you have all this data, I think the best answer is to actually fit a single large model, using Hierarchical Modeling rather than do it in two steps (generating a prior then fitting a model). This is basically the answer I gave to this question. I explain this a little bit more there.

In a hierarchical model you model each of the parameters you are interested in (for example, the location and scale parameters of the thumb lengths for a species) are drawn from a common prior distribution. The hyperparameters of the hierarchical model parameterize the common prior distribution of the parameters for the species', and you estimate the hyperparameters at the same time as the parameters you're interested in. The hyperparameters of course need their own prior distribution, but these can be relatively diffuse.

Answer 2

A useful way to incorporate data into a prior distribution is the principle of maximum entropy. You basically provide constraints that the prior distribution is to satisfy (e.g. mean, variance, etc.,etc.) and then choose the distribution which is most "spread out" that satisfies these constraints.

The distribution generally has the form $p(x) \propto exp(...)$

Edwin Jaynes was the originator of this principle, so searching for his work is a good place to start.

See the wiki page (http://en.wikipedia.org/wiki/Principle_of_maximum_entropy) and links therein for a more detailed description.

Answer 3

I propose the following solution to 2), and would appreciate feedback:

1. Data include mean, $Y$, sample size $n$, and standard error $\sigma$; calculate precision ($\tau=\frac{1}{\sigma\sqrt{n}}$) because it is required for logN parameterization by BUGS
2. data $Y\sim \text{N}(\beta_0,\tau)$
3. precision $\tau\sim\text{Gamma}(\frac{n}{2},\frac{n}{2\tau})$
4. diffuse priors
5. use $N(\mu=\beta_0, \sigma=\frac{1}{\sqrt{\tau}}$) prior

6. Here is the code:

   library(rjags)
   data <- data.frame(Y = c(1.6, 2.5, 1.8, 1.8, 1.7, 2.5), 
                      n = c(4, 4, 4, 3, 4, 3), 
                      se = c(0.2, 0.41, 0.24, 0.27, 0.2, 0.14))
   # convert se to precision
   data <- transform(data, obs.prec = 1/se)[, colnames(data)!='se'] 
   # write a bugs model
   sink(file= 'model.bug') #put following in file 'model.bug'
                           #i don't think sink() actually works like this 
   model 
   {
       for (k in 1:length(Y)) {
           Y[k] ~ dnorm(beta.o, tau.y[k])
           tau.y[k] <- prec.y * n[k]
           u1[k] <- n[k]/2
           u2[k] <- n[k]/(2 * prec.y)
           obs.prec[k] ~ dgamma(u1[k], u2[k])
       }
           beta.o ~ dnorm(3, 0.0001)
       prec.y ~ dgamma(0.001, 0.001)
       sd.y  <- 1/sqrt(prec.y)
   }
   sink()
   
   
   model  <- jags.model(file = "model.bug", 
                        data = data, 
                        n.adapt = 500, 
                        n.chains = 4)
   
   mcmc.object <- coda.samples(model = model, 
                               variable.names = c( 'beta.o', 'sd.y'), 
                               n.iter = 10000, 
                               thin = 50)
   summary(mcmc.object)
   

Update

I have revised this approach to compute a posterior predictive distribution. It required some modifications, mostly computing a posterior predictive distribution for an unobserved sample.

Details here:

David S. LeBauer, Dan Wang, Katherine T. Richter, Carl C. Davidson, and Michael C. Dietze 2013. Facilitating feedbacks between field measurements and ecosystem models. Ecological Monographs 83:133–154. http://dx.doi.org/10.1890/12-0137.1 pdf

Examples of this and simpler approaches here: https://github.com/dlebauer/pecan-priors/blob/master/priors_demo.Rmd