5
$\begingroup$

I am running Bayesian models to estimate the number of fruits on a plant, given the presence/absence of herbivores. I get a posterior distribution on each mean. I then run a separate model to estimate the number of seeds per fruit given the presence/absence of herbivores, again with posterior distributions on each mean. Is there a way to combine these posterior distributions to get an estimate of total seed production? Could this be as simple as multiplying the posteriors? Similar to calculating the average number of fruits per plant and the average number of seeds per fruit and multiplying those together to get a coarse estimate of seed production?

Or, is it OK to run both models simultaneously in one sampler and then work with the generated quantities. For example, using JAGS (or STAN):

for(n in 1:N){
    fruitNumber[n] ~ dnorm(fruitHat[n], sd_fruit)
    seedMass[n] ~ dnorm(seedHat[n], sd_seed)

    fruitHat[n] <- B0 + B1*Herb[n]
    seedHat[n] <- G0 + G1*Herb[n]
}

fruitHerb <- B0
fruitNoHerb <- B0+B1
seedHerb <- G0
seedNoHerb <- G0 + G1

totalSeedHerb <- fruitHerb * seedHerb
totalSeedNoHerb <- fruitNoHerb * seedNoHerb
$\endgroup$
  • $\begingroup$ Technically, it is fine to do both parts of the model simultaneously, but your likelihoods reflect the assumption that the number of seeds per fruit is conditionally independent of the number of fruits per plant, given herbivores. If that seems implausible, perhaps you should use a bivariate normal likelihood and estimate the correlation, instead of two univariate normal likelihoods which implicitly assumes the correlation is zero. $\endgroup$ – Ben Goodrich Oct 19 '14 at 4:13
  • $\begingroup$ If you're adding the random variables, and if you assume independence (though it may be a doubtful assumption) wouldn't the distribution you're after be a convolution of the two components of the sum? $\endgroup$ – Glen_b -Reinstate Monica Aug 27 '15 at 12:59
1
$\begingroup$

One approach to combine results from multiple models in a manner that reflects how well each model works on the dataset under analysis is Bayesian model averaging.

One simple way of approximating full Bayesian model averaging is to use weights proportional to $e^{-\text{BIC}_m/2}$ (or using the effective number of parameters à la DIC for a hierarchical model) for model $m=1,\ldots,M$. These weights approximate the posterior model probabilities, if a-priori all models are equally likely. Other weights are also possible, e.g. you might a-priori favor more complex models to get a AIC like behavior in the model averaging.

You can either average estimates (on a suitable scale) using these weights (and there are straighforward formulae published for combining uncertainty assuming a reasonably normal posteriors from each model) or perhaps in a more Bayesian fashion sample from the posteriors in proportion to these weights.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.