# Can I perform typical Bayesian contrast analysis, if my data takes the form of quantiles?

I'm working with some data from a cell analyzer, which takes millions of measurements and returns mean, median and 5% / 95% quantiles for the distribution. Preliminary model comparisons have indicated good fit to the data for a gamma distribution.

The experimental design is a straightforward ANOVA-type hierarchical model:

$$y_{ij} \sim Gamma(Sh, Sh/\mu_{ij})$$ $$log(\mu_{ij}) = \beta_{0} + \beta^{i}_{1}, i \in \{ 1, \dots 3 \} +$$ $$\beta^{j}_{2}, j \in \{ 1, \dots 3 \}, + \beta^{i}_{1} * \beta^{j}_{2}$$

In the Bayesian context the $\beta$ and Sh variables would have some sort of uninformative hyperparameters.

I want to investigate the contrasts between groups. I would prefer to do this in a Bayesian manner as I prefer my answer in the form of a probability mass.

What I would like to do is

• Fit each group's mean and quantiles to a gamma distribution using, for example, get.gamma.par
• Draw samples from the fitted distribution for each group
• Evaluate contrasts by subtracting samples of contrast groups as described in, for example, the Kruschke textbook

I am not sure if this procedure is legit because I am not generating posterior distributions using MCMC, but rather have essentially been given the posterior parametrization by this cell analyzer.

Can I still run a Bayesian contrast analysis given data in this form?

• So summarizing: instead of raw data you are given summary statistics (mean, median, 5%, 95% quantiles) and want to perform Bayesian analysis of such data, do I understand you correctly? What kind of model do you assume for this data (ignoring the fact that you have only summary statistics)? – Tim Nov 28 '17 at 10:20
• Yes you are correct. I'm not really assuming, but choosing a model via model comparison. It is all-positive data with some right skew (putting mean well higher than median). A gamma is quite a good fit and makes sense for all-positive measurements. – barnhillec Nov 28 '17 at 10:31
• Still, can you describe in greater details what is exactly the model you want to estimate? You've written that it is hierarchical and ANOVA-like, but could you give us more details? – Tim Nov 28 '17 at 10:33
• Post edited to contain a quantitative model specification, apologies for not including before. – barnhillec Nov 28 '17 at 10:53

Since same summary statistics may be produced from samples taken from different distributions, and estimating the parameters for gamma distribution may differ when using moments as compared to maximum likelihood, you need to take into consideration that what you have is just the approximate description of your data.

The first thing that comes to my mind in here is using approximate Bayesian computation, i.e.

1. sample the parameters from the priors you assumed,
2. simulate your data using the generative model,
3. calculate the summary statistics from the simulated data,
4. compare the simulated summary statistics $\hat S_\text{sim}$ with real summary statistics $S$ using some distance measure $\rho$, and when $\rho(S, \hat S_\text{sim}) < \varepsilon$ then accept the parameters, otherwise reject them.

The naive implementation should be very easy to code by hand in your favorite programming language, however this may be very inefficient, so comming up with more clever implementation may need some more research on your side.

You can find friendly introduction to ABC on Rasmus Bååth's blog (see also the accompanying video).