# Combining Posterior Distributions of Separate Models

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

• 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. Commented Oct 19, 2014 at 4:13
• 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? Commented Aug 27, 2015 at 12:59

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.