difference of means of non-normal data using Bayesian model I have two distributions, group 'c1' and 'c2' that are not normally distributed, but rather log-normal:

Now, I'd like to compare the means between 'c1' and 'c2' using bayesian methods, as described in John Kruschke's paper "Bayesian estimation supersedes the t test" ("BEST") see here
To model the distributions I used pymc3, simply using an uniform prior on the mu and sigma parameters

The trace plot looks fine (0 divergences) so at the end I look at the difference between the posteriors of mu and sigma parameters of the two models, but also at the difference between the means given by $e^{\mu + \sigma^2/2}$:

In all three plots it's clear that 0 is outside HDI, thus 'c1' are 'c2' different. So my questions is: Does this model make sense to analyze the difference in means between groups 'c1' and 'c2'. If yes, does it make sense to show the difference between the calculated posterior means or is showing the difference between $\mu$s and $\sigma$ sufficient?
 A: So let me give you a somewhat more straightforward alternative to the ROPE method in Kruschke’s paper.  You would need to solve the convolution of the two distributions to get to where you want to be.  There is an easier solution, though it is less comparable to a Frequentist procedure.
I would like to make an observation about Fisher’s “No Effect” hypothesis.  In this case, the null hypothesis would be that group membership has no effect on the population parameter’s location.  In other words, $\mu_1-\mu_2\equiv{0}$.  If this were a regression such as $z=\beta_xx+\beta_yy+\alpha$, and the belief was that the variable $y$ had no impact on z, then the null would be that $\beta_y=0$.  However, there is an alternative specification that would also work as they are both equivalents.  That would be to have the hypothesis that $z=\beta_xx+\alpha.$
The Bayesian equivalent is not to determine if a parameter is zero, but whether you can simply drop the variable as it is not probable that the variable is part of how nature generates the data.
The two Bayesian hypotheses would be: $$H_1:x\text{ is drawn from a single population}$$ versus $$H2:x_1\text{ and }x_2\text{ are drawn from different populations.}$$
In the first model, the group designation means nothing.  If you flipped a coin to divide kindergarteners’ classroom in half and measured their height, the heads versus tails groups would be artificial and were not part of how nature chooses sizes.
In the second model, the group designation is meaningful.  If you divided a group of same age adults from the same population by birth sex, you would find that they were different populations because nature does assign height by sex.
As it is a Bayesian model, you will have to assign a prior probability to model one and model two.  It can either be your prior probabilities or, if you have an opponent, then you could use their prior probabilities.
If $f$ is the log-normal likelihood function, then the joint likelihood for model one is $$\prod_{i=1}^If(x_i|\mu;\sigma^2).$$
If $c$ is a categorical variable denoting group membership, then the join likelihood for model two is $$\prod_{j=1}^Jf(x_j|\mu_1;\sigma^2_1;c_j=1)\prod_{k=1}^Kf(x_k|\mu_2;\sigma^2_2,c_k=2), J+K=I.$$
You still need to set your prior distributions.  I also assumed independence of draws.
I would point out a couple of warnings.  First, when you use uniform prior distributions, there is a chance, particularly as the dimensionality of variables becomes three or higher, that the posterior will not sum to one.  You are comparing two groups, so it isn’t an issue here.  Second, a uniform distribution in the raw space is not uniform in log space.  You would be adding information, possibly unintentionally.  Third, it no longer matters about finding some area around zero.  You will not be subtracting parameters anymore.  Fourth, you would not be comparing mean values anymore.
You can have two lognormal distributions with the same population mean, but with different values for $\mu$ and $\sigma^2$.  The Bayesian method would not determine if $$\mu_1+\frac{1}{2}\sigma_1^2=\mu_2+\frac{1}{2}\sigma_2^2.$$  Instead, it would directly determine if they are from the same population.  After all, variables could be drawn from different populations but would have the same population mean if $$\mu_1=\mu_2+\frac{1}{2}\sigma_2^2-\frac{1}{2}\sigma_1^2.$$
Remember that $\mu$ is the mean only after logarithmic transformation; otherwise, it is the mode.
Suppose you set $\mu_1=1$, $\mu_2=1.1$, $\sigma_1=1$ and $\sigma_2=0.8944272$. In that case, the t-test will say they are the same if performed in the raw data space, without the log transform. Simultaneously, the Bayesian method will say they are different, with just a modest sample size.
A: Given that the distributions are log-normal, it would seem more natural to test the ratio of the means (the difference of the logs).
A: If your goal is to estimate the magnitude of the difference between groups, then your procedure is fine (in principle). You can decide whether it's more meaningful to consider $\mu_1-\mu_2$ or $\mu_1/\mu_2$ based on whether the original scale or the log() scale is more sensible to deal with. The HDI's won't be exactly the same (because HDI's aren't invariant under nonlinear transformations) but will probably yield the same conclusion. Same for $\sigma_1 - \sigma_2$. If you want to make a decision based on the magnitude of difference, you need a decision procedure; see this article (disclosure: I wrote it).
If your goal is to test the null hypothesis that the difference is zero then you're probably wanting to frame the question as a Bayesian model comparison. This can be a very useful approach but it can be very sensitive to the choice of prior distribution, unlike parameter estimation. For caveats about this approach, and references for how to pursue it, see this article (which I also wrote).
