pMCMC convergence diagnostics

[kingaa]

The following question was sent to me by @CorderRM. It raises some interesting issues about particle MCMC, so I include it here despite the fact that some of the context is lacking (he asked that I not post the model or the data).

We're starting to use pomp and its pmcmc function to calibrate a stochastic model we developed and I'd appreciate it if you could help me answer some simple questions. I attach here some preliminary results we obtained. Basically, we want to estimate two parameters and their respective credible intervals.
Both parameters are probabilities and should vary from 0 to 1. For this reason, they are logit transformed. If I understood correctly the output of the model, 2000 iterations (Nmcmc) were clearly not enough for convergence and I must increase it. Also related to the MCMC process, I could not find how to check the acceptance ratio in order to verify if the "proposal" input should be better adjusted. Regarding the number of particles (Np), I used 200. I was wondering how to check if the number of particles is reasonable in this case.
I used parallel computation for each chain (foreach function) and I was also wondering if it is not possible to parallelize other processes within the estimation in order to speed it up. It seems that only one of the chains has been initialized from the point specified by rinit. How were the other starting points selected by pomp?
Here are my questions:

1. I was just able to find how to parallelize the code with the foreach function for each chain. Because of this, the process is really slow. Is it possible to parallelize other processes within the pMCMC framework?
2. My preliminary results show that the number of iterations chosen is small once the trace of prob2 did not converge (Nmcmc = 2000). I was wondering how I can make sure that the number of particles is sufficient in this case (I used Np=200).
3. How can I verify the MCMC acceptance rate within the pmcmc framework in order to make sure that the proposal was reasonable for my application (proposal = mvn.diag.rw(c(prob1 = 0.2, prob2 = 0.2)))?

I replied:

I'll first try to answer your questions, then to offer some advice based on my observations about your problem.

1. You ask about how to accelerate the code, suggesting that perhaps one can parallelize it. I suspect that you will be able to achieve substantial speedup without the headache of parallelization by making your code more efficient. The bulk of the computational effort is spent in your rprocess function, which for runs of the size you describe is executed hundreds of thousands of times on each parallel processor. For example, at present, you create several data frames each time that function is called. This appears mostly to be for the sake of organizing your variables (i.e., making the code easier to read) but creation of data frames is expensive (this is not one of the proper uses of that structure) and it is questionable whether it actually improves the readability. I, for one, find it hard to read long arithmetic expressions in which the variables have long and redundant names. I would suggest instead using variables with elaborate (or not so elaborate) names to perform the calculations you do in rprocess. Thus for example states$HB_M[1] could be states_HB_M[1] or even HB_M[1], with no loss in readability and some gain in performance.
   Once you've done this, it will be relatively straightforward to take the next step and encode your process model using C snippets. Every pomp basic component can be furnished either as an R function (as you have) or as a C snippet. While the former is useful when people are learning about POMP models and inference algorithms, the C snippets are much faster (often 20-100x faster) so that they are usually necessary in actual scientific applications.
   This might seem like a lot of trouble when you could just parallelize the code. However, the particle filter---which is the basis for pmcmc---is not easily parallelized in an efficient way. This is because there is a need for communication among the particles quite frequently. Generally speaking, the less frequent the required communication, the more efficient is the parallelization. For this reason, pomp does not offer a parallel particle filter.

2. You ask how you can diagnose the performance of pmcmc. This algorithm is essentially Metropolis-Hastings (MH) using the particle filter to obtain the likelihood. Since the particle filter returns a noisy (though unbiased) estimate of the likelihood, it will mix more slowly than it would if the likelihood were given exactly. In particular, the algorithm can become stuck if the noise is too great. In general, one needs the same degree of precision on the likelihood estimates for the MH that one needs for statistical inference, i.e., no more than roughly 1 log unit of Monte Carlo noise is permissible. So diagnosis is not too difficult: simply ensure that you have sufficiently many particles that the log likelihood estimates of the particle filter are repeatable to within about 1 log unit.

3. You ask about convergence and mixing diagnostics for pmcmc. The traces function, operating on the output of a pmcmc command, gives an object of class mcmc or mcmc.list; these are defined in the popular coda package. All the diagnostic tools from coda are therefore at your disposal. In particular, one can examine acceptance rates and effective sample sizes, run convergence diagnostics, and plot and manipulate the MCMC traces in various ways using the tools of that package.

4. Earlier, you asked whether rinit was initializing all the chains in the same way. It does (or should, at any rate). If you have reason to believe that it is not, please consider filing a bug report.

Now for some more general remarks. I realize there is a risk in offering unsolicited advice: I hope that you can accept it in the spirit in which it is offered.

1. Your rmeasure and dmeasure components do not match. You may have reasons for this, but I point it out because unless they do match, simulations from the fitted model will not in fact be valid draws from that model.

2. You seem to be using the default priors, which in this case are effectively improper uniform on the logit scale. This corresponds to a bathtub-shaped distribution---Beta(0,0)---on the probability scale. I strongly advise you to spend some time thinking about reasonable priors. In particular, it is unlikely that you have an abundance of information in the data (13 measurements after all), so your prior may be driving the results. My advice is to think of the prior as just another component of the model and to approach estimating it in the same way you would any other feature of the model.

3. The structure of your rprocess function suggests that you are imagining a deterministic version of your process model. In such a case, there is no advantage to using more than Np=1 particle. A more standard, and probably safer, way of proceeding would be to encode the deterministic model as a skeleton, and then compute the likelihood using traj_objfun and wrap it in a custom MH..

4. There is some danger in providing default arguments for the functions you furnish as basic model components. It would be safer to avoid doing so, since a class of important errors in coding will show up as runtime errors in that case.

5. Finally, I am struck by the complexity of the model you are using to attempt to explain what are, at the end of the day, very few data. Do you have reason to believe that the data contain information about the parameters you are seeking to estimate? And what of the many parameters you are assuming are known? Does your certainty about those parameters truly dwarf that about the two parameters you are estimating?

[CorderRM]

I really appreciate the completeness of your reply.

1. You ask about how to accelerate the code, suggesting that perhaps one can parallelize it...
R: Yes! Thank you for pointing it out. Just by removing the data.frame structure, the code is already much faster.

2. You ask how you can diagnose the performance of pmcmc...
R: I would like to understand this better. You say: "no more than roughly 1 log unit of Monte Carlo noise is permissible". Let's assume a hypothetical model that has just one parameter (p1) in which the best loglikelihood is found when p1=x. Dou you mean that the MCMC noise (proposal) for p1 should not exceed log(x)? Then, regarding the pmcmc, you say "simply ensure that you have sufficiently many particles that the log likelihood estimates of the particle filter are repeatable to within about 1 log unit". If my proposal is somehow close to log(x) and my initial value for p1 is close enough to x, I should not worry about too many particles once the loglikelihood estimates of the particle filter will be repeatable to within about 1 log unit of the best value for p1. Is that correct? In that case, are 100 particles, say, enough?

3. You ask about convergence and mixing diagnostics for pmcmc…
R: Thanks for pointing out the coda package for this purpose.

4. Earlier, you asked whether rinit was initializing all the chains in the same way. It does (or should, at any rate). If you have reason to believe that it is not, please consider filing a bug report.
R: I was wrong about this. My point was not about the initial state values but about the initial parameter values. If I’m running more than 1 chain and defining just one initial condition for parameter initial values, how does pomp choose the other ones? I could notice that each chain starts with different parameter values.

Now for some more general remarks.
1. Your rmeasure and dmeasure components do not match…
3. The structure of your rprocess function suggests that you are imagining a deterministic version of your process model…
R: In my model I was assuming that the measurement noise is 0. (This assumption can be relaxed later on.)

        rmeasure = function (CFPp, ...) {
          c(CFPp_model = CFPp)
        },

However, instead of noise, I had many stochastic processes within the model (rprocess). For example, for each time step, the model is sampling (from a distribution defined within the rprocess) the number of eggs each mosquito will oviposit according to their genotype. It means that, for example, at time step 1, individual with a certain genotype might generate 100 or 110 eggs and, at time step 2, an individual with the same genotype might generate 90 or 105, say. Therefore, at each iteration the model output is changing. Other stochastic processes are also considered within model.

By assuming measurement noise equal to zero, I’m not assuming that the output of my model should be exactly the same as my data. Due to the counting data experiment (600 individuals are checked at each generation) and to the stochasticity of the dynamics, I assumed that the model's fit residuals follow a Poisson distribution.

        dmeasure = function (CFPp_data, CFPp_model, ..., log) {
          c(ll = dpois(round(600*CFPp_data),round(600*CFPp_model),log = T))
        },

Note that CFPp_model is the output from rmeasure and an input to dmeasure (to be honest, I’m not sure if this setup is working properly within pomp framework). Due to the inherent stochasticity of all accounted processes, I was not able to encode a satisfactory deterministic model as a skeleton (to follow a more standard way of calibrating it). I would really appreciate your thoughts about this point (and any others that you think would be relevant to this analysis).

2. You seem to be using the default priors, which in this case are effectively improper uniform on the logit scale…
R: Thank you for pointing it out. I bounded the prior distributions within a comfortable range for my estimates. It’s still not the ideal, but definitely improved the alg.

4. There is some danger in providing default arguments for the functions…
R: Thank you for pointing it out.

5. Finally, I am struck by the complexity of the model you are using to attempt to explain what are, at the end of the day, very few data. Do you have reason to believe that the data contain information about the parameters you are seeking to estimate? And what of the many parameters you are assuming are known? Does your certainty about those parameters truly dwarf that about the two parameters you are estimating?
R: I really appreciate your comment. I haven’t mentioned so far that I’ve designed the model before having the complete dataset from all experiments. These 13 points are mean values of 3 experiments replicates. I’ll have access to other experiments replicates starting from different points. Moreover, the known parameters I’m assuming in the model were all measured by the same laboratory for the same species. Because I have distributions for all of the parameters, I’m considering the stochastic model. In the end, I still have to account for some other important biological phenomenon, but I’m expecting to find the probability I'm searching for around 0.9. I can see that the model is describing the data experiment fairly enough by playing with the 2 parameters I’m trying to estimate (figure below).