Suppose I have 2 models which represent the same chemical system. Both models have the exact same reactions, however the second model has a slight variation. The variation is regarding an intervention, for model 1 we assume that we have 1 intervention in which $1000$ of molecule $X$ is dumped into our chemical system in one big go, whereas model 2 assumes that we dump $200$ molecules of $X$ in in 5 even time intervals (so total number of molecules dumped into each system is equal). I want to see what affect this variation has on one of my rate constants in my chemical system, let's call this rate constant $c$.
I use the gillespie algorithm to simulate both models, the models are stochastic hence each result is unique from the last. I assume to have perfect information so I know the complete sample path, i.e. know the time and type of each reaction that occurs. As I know the time and type of each reaction I am able to calculate how many times each reaction occurs and update my prior beliefs of my rate constants to my updated posterior
My question is: What is a good Bayesian test to perform to test difference in rate constants or something along those lines.
Further detail: Following on from the 1st paragraph where I outline the models. A very important assumption I have about my rate constant $c$ is that it has a gamma distribution. Also for this reaction with rate constant $c$ I assume my prior beliefs to be the same for both models.
What I have done: So after thinking about this for a while, I thought a good place to start off would be to do 50 simulated repetitions of this system for each model. This allows me to get 50 unique posterior distributions for each model for my rate constant $c$ based on the data which was simulated. Now I know, as I have a gamma distribution, that the expectation (ie mean) for each updated posterior is $\frac{a}{b}$ for a posterior $\sim Ga(a,b)$ and variance $\frac{a}{b^2}$. Maybe I could test between the ranges between the expectations/variances for my rate constant $c$ between the 2 models.
I hope this is clear, if there is any part that is not I will make sure to edit the question.
A second, perhaps better, approach: Maybe what I could do is see what the bayes factor is given a data set when comparing my 2 models. I was wondering how I would go about this, as every time I've attempted to do bayes factors it seems to go horribly wrong
Edit: I was looking to do something along the lines of these plots (https://stackoverflow.com/questions/11546256/two-way-density-plot-combined-with-one-way-density-plot-with-selected-regions-in?rq=1) but not entirely sure how to go about this. Should I take the 50 variances for each model and do a density plot between the 2?