Estimate particle filter fit

[jandraor]

Dear Professor King,

First of all thank you so much for this package & the related documentation. I've been following your course "Simulation-based Inference for Epidemiological Dynamics" & I've found it quite valuable to my research.

Recently, I've been working with an SEIR model (ODE) whose effective contact rate is modelled as Geometric Brownian Motion with no drift. Essentially, a non-negative random walk. I used iterated filtering to estimate a few parameters from that model. So far, the results are encouraging, i.e. quadratic shapes for the MLE. From this MLE, I want produce a "fit" & compare it to the real data. To do so, I plug the MLE into the particle filter:

pomp_mdl %>% 
  pfilter(params = mle_estimate, Np = 5e5, pred.mean = TRUE, 
          save.states = TRUE) -> fit

C represents the weekly incidence in my model & it's a state. At each time step (week 1, week 2, week 3), I compute from the 500,000 particles the mean & 2.5 % & 97.5% quantiles of that state, which I assume will be the confidence (credible?) intervals of the fit. If I plot those results against the data, it's almost perfect. Too good to be true.

However, that mean is different from the "pred.mean". I think I've got a misunderstanding in my concepts.

In short, what would be the appropriate way to estimate the mean & uncertainty intervals of a particle filter fit?

Thanks.

[kingaa]

Thanks @jandraor, for the kind words, and thanks for bringing up a good question.

If $X_k$ and $Y_k$ represent the latent state and observable processes at time $t_k$, respectively, then one can ask about the distribution of $X_k$ given $Y_1,...,Y_j$, for any $j$. For $j = k$, and $j=k-1$, pfilter computes some aspects of these distributions that can be used to diagnose the filtering. In particular, with the save.states option set, pfilter returns an estimate of $X_k$ given $Y_1,...,Y_k$ for each $k=1,...,K$. This is called the filter distribution. With the filter.mean option set, it returns an estimate of the expectation of this distribution. With the pred.mean option set, it returns an estimate of the expectation of the prediction distribution, i.e., $X_k$ given $Y_1,...,Y_{k-1}$.

Your intuition that the filter distribution is "too good to be true" is a valuable one and it is fortunate that you respect your intuition. You haven't made it entirely clear why you are doing what you are doing, but if I guess correctly, you are asking the extremely pertinent question, "How does the model explain the data?" and perhaps the closely related, "How well does the model explain the data?". The royal road to answering these questions lies through gaining an understanding of the distribution of the data under the model. To be precise, the fitted model is nothing other than a probability distribution for $Y_1,...,Y_K$; we think of the data as a single draw from this distribution. To explore it, the most straightforward path is through simulation, unconditional on the data. You can obtain additional draws from this distribution via the simulate function in pomp. Thus, for example, were you to simulate 1000 realizations using simulate(...,nsim=1000), you could plot the quantiles of this distribution and compare with the data. (By the way, while doing so might be revealing, there is nothing that replaces simply plotting a small number, say 10–20, realizations and overlaying the data.) More generally, for any statistic you can compute on the data, it can be computed on a collection of simulations, and the resulting distributions compared; this is facilitated by the probe function provided by pomp.

Now, you may (also) be interested in asking the very different question, "Given the model, what was the path of the latent process actually taken in the single realization that was observed?" For example, you wish to know what the value of your variable C was on a particular day or set of days. This is a question about $X$ conditional on the data. For this purpose, the filter distribution returned by pfilter may be of some use. In many applications, however, one wants to know the distribution of $X_k$ given $Y_1,...,Y_K$. This is called the smoothed distribution. To obtain this, more computation is required: for example, one can use pmcmc and extract trajectories using filter.traj to obtain samples from the smoothed distribution.

[jandraor]

Dear Professor King (@kingaa),

I appreciate your thorough response & the time you devoted to this. It has answered some of my doubts but also raised new ones. Here are the model equations to clarify the context of this work.

Csnippet("
    y = rpois(C);  
 ") -> rmeas

Csnippet("
  lik = dpois(y,C,give_log);
") -> dmeas

Csnippet("
    S = 5e6 - I_0;
    E = 0;
    I = I_0;
    R = 0;
    B = B_0;
    C = 0;
  ") -> rinit

 Csnippet("
    double dW     = rnorm(0, sqrt(dt));
    double lambda = B * I / N;
    S-= (S * lambda)*dt;
    E+= (S * lambda - sigma*E)*dt;
    I+= (sigma * E - gamma*I)*dt;
    R+= (gamma*I)*dt;
    B+= alpha*B*dW;
    C+= (rho*sigma*E)*dt;
  ") -> SEIR_GBM_step

The calibration of this model is in the context of a COVID-19 outbreak. As such, the assumption that the effective contact rate (B) is constant no longer holds, even for a short period given the quick measures imposed since first cases reported. For this reason, my purpose is to approximate the trend of a time-varying B, which translates into a time-varying basic reproduction number. I propose Geometric Brownian Motion as a vehicle to unravel such a trend. By no means, I conceive B as a random walk. For this reason, I don't think that this model is able to answer "How well does the model explain the data?". Considering that a random walk can produce a broad arrange of trajectories & incidences, the observed data has a tiny probability of being drawn. To complement this point, I think this also the reason the ESS goes relatively low. To 10% of the ensemble size.

Consequently, the focus of this research is on the question "Given the model, what was the path of the latent process actually taken in the single realization that was observed?". My approach is that the particle filter can answer that question. However, to run a particle filter, I need a complete set of initial values. In this case, I'm assuming I don't know I(0), B(0), alpha & rho. I followed your tutorial to estimate such values via Iterated Filtering. Fortunately, there's concentration of measure around a single region & I can compute confidence intervals using the profile likelihood approach. Subsequently, I plug the MLE into the particle filter to obtain an estimate of the time-varying B. With this, I can compute the time-varying basic & effective reproduction number.

Based on the above, I've got the following questions:

- Is this approach enough to answer the research question? When I said "too good to be true", it was because the filtering distribution (mean) fitted very well the data with narrow uncertainty intervals. I estimated the mean from the 500.000 particles returned by saved.states & calculated intervals via the quantile function in R.

• 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. Should I estimate the smoothing distribution to answer my research question?

• Finally, if the particle filter approach is enough, the variance in the fit is due to the randomness in the structure of the process model. However, I'm fixing parameters to the estimated MLE's. Should I run the particle filter from the complete distribution of the confidence intervals to account for the uncertainty in parameter estimates? Or should I commit to the MLE's?

Thanks.

[kingaa]

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.

[jandraor]

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

[kingaa]

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

[jandraor]

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

[kingaa]

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.

[jandraor]

Hi Prof King,

It's fascinating to see how disinterested guidance can contribute to the career of a researcher. In my case, this discussion shaped to a large extent what later would become a paper.

Thank you.

Here is the result:

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010206