Ah, very interesting. I have been intrigued by these "semi-parametric" methods that constrain some aspects of model form while leaving others only weakly constrained. I have done some work with them, using both spline-based and diffusion-based approaches. The latter is very closely related to your approach; see Wale et al. (PNAS 2019).

In short, I do think the approach is capable of yielding answers to your question. In particular, you extract the filter distribution to examine the sample paths of your B variable that, under the model, generated the data. This is a "hindcast". B is constrained to be continuous, with a degree of smoothness controlled by the intensity parameter alpha. Your intuition that this is too good to be true is a bit off the mark in this case: the data do constrain the path of B quite strongly. Your specific questions concern the two important points left to consider.

First, as you note, there are sources of uncertainty that are as yet unaccounted for. In particular, uncertainty in the constraints you are applying, i.e., the form of the model. It is useful to think of this as having a component associated with parameterization and another component associated with structure. If you were, for example, to consider the B-paths consistent with the data with somewhat different but equally or nearly equally well supported parameters, you would see larger differences than you do. If, moreover, you were to consider alternative models that represent modified constraints, you would see further differences. In our 2015 paper we present an approach to the inclusion of parametric uncertainty in both hindcasts and forecasts. In essence, this approach samples the neighborhood of the MLE, weighting parameters by their likelihood. This can be viewed as a kind of empirical Bayes approach. The code for this is in a datadryad archive.

The second point, of which again you are aware, is that, as I described in my first reply in this discussion, the filter distribution applies only a portion of the data to each portion of the hindcast. It would be better to use all the data, i.e., to obtain the smoothing distribution. Note that this will not necessarily reduce the uncertainty in the hindcasts: increases can be interpreted as indications of model misspecification. In the 2019 PNAS paper, and in full detail in the data supplement for that paper, we describe how we used PMCMC to obtain the smoothing distribution. This is just an implementation in pomp of the recommendation in the original PMCMC paper.

You mentioned pMCMC as a method to estimate the smoothing distribution. I presumed (perhaps incorrectly) that pMCMC was the Bayesian counterpart of the frequentist Iterated Filtering for the estimation of constant parameters.

Yes, PMCMC can be used to estimate parameters and in that sense is an alternative to IF. As such, it can be much less efficient in many applications, since one has to get the likelihood estimate with high precision (i.e., many particles) in order for the MCMC to mix. On the other hand, it is perhaps the best method for obtaining samples from a posterior distribution for POMP models. I typically recommend that users who have a good reason for performing a Bayesian analysis first identify the region(s) of high likelihood using IF before commencing PMCMC.

Thanks Professor King. Certainly your responses have been invaluable to keep going with my research. I will go through the recommended material & run some more experiments. I'll post again in case of any further questions. Hope to keep them to a minimum.

Thanks for your time.

Jair

Sounds good. Please consider sharing the results of your experiments, however briefly, with the community here!

Dear Professor King,

Thanks for the guidance & for sharing the code of your results. Very much appreciated.

I've got a question regarding uncertainty (I'm agnostic here) interval estimation. I'm estimating I(0), B(0) & alpha. I followed your code to estimate parameter uncertainty using the profile likelihood method for each parameter. Here is Alpha's profile. geom_smooth shows the quadratic trend in the profile. The red line indicates chi-square cutoff. However, the points by themselves don't present a smooth graph. Also, I get different MLEs & cutoffs for every profile.

On the contrary, when I collate all profiles in a single data frame, I get smoother marginal distributions (see graph below) and a single MLE.

1.1 Could I use a single cutoff (profile likelihood) on the collated dataframe rather than on single profiles with different MLEs & cutoffs? Would that be methodologically sound?

In my previous experience with MCMC methods, the procedure for estimating credible intervals is to use the quantile function on the samples?

In this case, the collated data frame is prof_est_df.

1.2 Could I do the following?

quantile(prof_est_df$alpha, c(0.025, 0.975))

Regarding the uncertainty in the filtering distribution, I followed your 2015 paper. Specifically, I implemented the code below for producing a model fit. Recall that prof_est_df is a data frame which contains the profile of each estimated parameter. run_particle_fit is a function that runs pfilter with save.states = TRUE for every row of params_df & returns a data frame of the hidden states particles over time & its respective likelihood (loglik).

mle_boundary <- prof_est_df %>% filter(loglik > max(loglik, na.rm = T) - 6)

# Here I select the unknown variables B_0, I_0, alpha
unk_par_MLE <- mle_boundary %>% select(unk_names)

sapply(unk_par_MLE, function(col) {
  min_col <- min(col)
}) -> lower_lims

sapply(unk_par_MLE, function(col) {
  max_col <- max(col)
}) -> upper_lims

set.seed(813292)

params_df <- sobol_design(lower = lower_lims,
                      upper = upper_lims,
                      nseq  = 200)

hindcast_df <- run_particle_filter(pomp_mdl, params_df) %>%
  mutate(weight = exp(loglik - mean(loglik)))

summary_hindcast <- hindcast_df %>%  group_by(time, var) %>%
  summarise(value = wquant(value, weight),
                    type = c("lower_lim", "median", "upper_lim")) 

Regarding this, I've got the following questions:

2.1 Why do you cut the log likelihood off at exactly 6? What's the rationale for that number?

2.2 In the Sobol design, we assume the log likelihood surface as a hypercube. However, this is not entirely true given the correlations among parameters (see graph above). I presume wquant accounts for parameter correlations based on the penalisations given by the likelihood weights. In my previous experience with MCMC methods, I simply sample rows (in this case, it would be from params_df). By doing so, I don't have to use weights later. Would that procedure (sampling entire rows) be valid here or is it better the weight penalisation?

I followed your procedure of weight penalisation & here are the results:

Thanks in advance for the time you've dedicated to this. I've found the POMP package & your online lectures quite useful to my research.

Regards,

Jair

I've got a question regarding uncertainty (I'm agnostic here) interval estimation. I'm estimating I(0), B(0) & alpha. I followed your code to estimate parameter uncertainty using the profile likelihood method for each parameter. Here is Alpha's profile. geom_smooth shows the quadratic trend in the profile. The red line indicates chi-square cutoff. However, the points by themselves don't present a smooth graph. Also, I get different MLEs & cutoffs for every profile.

Yes, the algorithms for estimating the parameters and for computing the likelihood are stochastic algorithms. Therefore, there is some Monte Carlo variability associated with them. To make the Monte Carlo error sufficiently small, one has to put in a sufficient amount of effort. For example, since differences of 1 log unit of likelihood are statistically meaningful, one typically has to make sure that the Monte Carlo error on log likelihood estimates is less than 1 unit. In your first figure, you are getting close to this, but would need to put in a bit more effort to achieve the requisite degree of precision. When using pfilter for likelihood estimation, one does this by increasing the number of particles Np.

Not all of the error in your first figure is Monte Carlo likelihood-estimation error, however. In addition, there is estimation error: some of the variability in the points shown is due to their sitting at random points. To reduce this, one can either spend more effort on each estimation attempt, or do more estimation attempts. In practice, the latter seems to work best. With more estimation attempts, the consistency (or lack thereof) among attempts gives confidence (or lack thereof) in the accuracy of the estimates. In your first figure, there are really too few points to be sure of the precise location of the MLE to within, say, 0.5 units of α.

This explains why, as you describe below, you have more confidence in the shape of the profile when you put everything into a single data frame. I do point out, however, that your second figure does not show how the log likelihood varies.

In passing, I mention that we have recently developed an approach to improve the robustness of profile estimation for POMP models by making adjustments for the Monte Carlo error. See Ionides et al., "Monte Carlo profile confidence intervals for dynamic systems" (J. R. Soc. Interface, 2017).

On the contrary, when I collate all profiles in a single data frame, I get smoother marginal distributions (see graph below) and a single MLE.

Yes, and this is our usual recommendation (e.g., the pimp my pomp wiki and "Simulation-based Inference for Epidemiological Dynamics", Lesson 4. Making the plots as you have in your second figure is what we call a "poor man's profile", or would be if you included the log likelihood as one of the variables.

1.1 Could I use a single cutoff (profile likelihood) on the collated dataframe rather than on single profiles with different MLEs & cutoffs? Would that be methodologically sound?

This is not perfectly sound (hence the "poor man's" qualifier in our phraseology). As we explain in "Simulation-based Inference for Epidemiological Dynamics", Lesson 3, to compute a profile likelihood for one parameter, one has to maximize over the other parameters. However, one could do a round of estimation, using each of the points (or perhaps only points near the upper envelope) as starting points, leaving the value of the profiled parameter fixed. To the extent that the maximization is successful (i.e., keeping the remarks above about Monte Carlo error in mind), the result would then be the basis of a valid profile likelihood.

The poor man's profile is very useful, however. Suppose we project all these parameter/likelihood estimates onto the θ-ℓ plane, where θ is any parameter and ℓ is the log likelihood. The upper envelope of the points is a lower bound on the profile likelihood for θ. It follows that, as you add to the database progressively better-refined estimates, the poor man's profiles approach the true profiles.

In my previous experience with MCMC methods, the procedure for estimating credible intervals is to use the quantile function on the samples?

In this case, the collated data frame is prof_est_df.

1.2 Could I do the following?

quantile(prof_est_df$alpha, c(0.025, 0.975))

You haven't said exactly where the points in your data frame come from. I'm guessing that they are the results of many parallel and independent IF2 runs. Is that right? If so, then you are not justified in assuming they are samples from a posterior distribution. In particular, IF2 is designed to maximize the likelihood, not to approximate a posterior distribution. Accordingly, if everything is sufficiently regular and one had a flat prior, the IF2 estimates would congregate near the maximum of the posterior, but they could be made to lie arbitrarily close to the maximum and their spread would bear no necessary relation to the width of the posterior.

On the other hand, if these points represent samples from the posterior distribution, obtained for example, by using pmcmc, then one could construct one sort of credible interval just as you describe.

Regarding the uncertainty in the filtering distribution, I followed your 2015 paper. Specifically, I implemented the code below for producing a model fit. Recall that prof_est_df is a data frame which contains the profile of each estimated parameter. run_particle_fit is a function that runs pfilter with save.states = TRUE for every row of params_df & returns a data frame of the hidden states particles over time & its respective likelihood (loglik).

mle_boundary <- prof_est_df %>% filter(loglik > max(loglik, na.rm = T) - 6)

# Here I select the unknown variables B_0, I_0, alpha
unk_par_MLE <- mle_boundary %>% select(unk_names)

sapply(unk_par_MLE, function(col) {
  min_col <- min(col)
}) -> lower_lims

sapply(unk_par_MLE, function(col) {
  max_col <- max(col)
}) -> upper_lims

set.seed(813292)

params_df <- sobol_design(lower = lower_lims,
                      upper = upper_lims,
                      nseq  = 200)

hindcast_df <- run_particle_filter(pomp_mdl, params_df) %>%
  mutate(weight = exp(loglik - mean(loglik)))

summary_hindcast <- hindcast_df %>%  group_by(time, var) %>%
  summarise(value = wquant(value, weight),
                    type = c("lower_lim", "median", "upper_lim")) 

Regarding this, I've got the following questions:

2.1 Why do you cut the log likelihood off at exactly 6? What's the rationale for that number?

I'm afraid I don't recall, specifically. I should point out that the goal in this piece of code is to sample the neighborhood of the MLE. Lowering the threshold will enlarge the neighborhood, but the weights of the additional points will be lower, since the likelihoods are lower. With a threshold set at e^-6 below the maximum, the smallest weights are about 1/400 of the highest weights. One can always do more work, with diminishing returns.

2.2 In the Sobol design, we assume the log likelihood surface as a hypercube. However, this is not entirely true given the correlations among parameters (see graph above). I presume wquant accounts for parameter correlations based on the penalisations given by the likelihood weights. In my previous experience with MCMC methods, I simply sample rows (in this case, it would be from params_df). By doing so, I don't have to use weights later. Would that procedure (sampling entire rows) be valid here or is it better the weight penalisation?

It does not assume anything about the shape of the likelihood surface. Two assumptions are made about the Sobol' design: (1) that the box is sufficiently large that it embraces the entire region of high likelihood and (2) that the sampling is sufficiently dense to represent the full range of model behaviors consistent with the data. Item (1) can be checked easily, using the results of the profile calculation, though if there is another mode in the likelihood surface that we have not seen, we of course will be ignorant of it. Item (2) is more subtle in general, though in our case, since the model was so simple, it was not difficult to convince ourselves that we were not missing anything.

An alternative approach would be to sample from the posterior, using pmcmc for example. Then, assuming your sampling has converged and you have thinned it sufficiently, one could simply plug in the posterior samples to the simulator to generate simulations and sample one trajectory per sample to construct the smoothed trajectories.