Reposting the model formulation for reference.
So I can already see some possible non-identifiabilities (you should however verify those with simulations, e.g. using the deSolve
package, I am just guessing here):
- if \theta is close to 0, h can take almost any value (this should be visible when plotting posterior of \theta against h as some sort of “funnel” or other weird shape). In plain words: if the immune system doesn’t recognize the parasite at all we can’t learn how quickly it kills it.
- if h is large, \theta can take almost any value. In plain words: if the immune system takes forever to kill the parasite we can’t learn how quickly it recognizes it.
- \theta and h only ever influence the results via g(P,\theta,h) = \frac{\theta P}{1 + \theta P h}, so, especially when P ends up not changing much, \theta and h can vary a lot while giving the same g(P,\theta,h). In plain words - if we don’t observe the parasite population changing, we can’t learn whether it is because it is quickly recognised but killed slowly or slowly recognized but killed quickly (or something in between)
- when the observed H stays close to zero for most of the time, c, \theta, h, \mu end up not really influencing the likelihood (and can thus take basically any value). In plain words, if we observe very few immune cells, we can’t learn much about their dynamics.
- when the observed P stays close to zero most of the time, r, \theta, h, c end up not really influencing the likelihood. In plain words, if we observe very few parasites, we can’t learn about their dynamics.
If those are indeed the problems, it is still hard to figure out what to do about it. Reparametrizing via q = \theta h might help a bit. Restricting the priors for \theta and h to avoid the degenerate cases might also help (but only if you can be sure those values really don’t occur in your experiment)
Also avoiding data that don’t demonstrate all of the dynamics can help: e.g. if P is well fit by an exponential curve, it admits the solution with \theta = 0 and the experiment thus can’t inform the other parameters. Similarly when P is roughly constant, the experiment just can’t inform both \theta and h.
One way to approach this in a bit more systematic way would be to try a set of models of increasing complexity, something like:
a)
b)
and c) the original model. (more intermediate steps can also be conceived, e.g. b), but ignoring the c parameter or the full model but ignoring \mu)
If one of the simpler models fits the data well, it would mean that you probably can’t really learn much new information by fitting a more complex model, as the more complex models have many ways to emulate the same dynamics as the simpler models.
Does that make sense?