Getting started with pomp

Structure of a POMP

As its name implies pomp is focused on partially observed Markov process models. These are also known as state-space models, hidden Markov models, and stochastic dynamical systems models. Such models consist of an unobserved Markov state process, connected to the data via an explicit model of the observation process. We refer to the former as the latent process model (or process model for short) and the latter as the measurement model.

Mathematically, each model is a probability distribution. Let $Y_n$ denote the measurement process and $X_n$ the latent state process, then by definition, the process model is determined by the density $f_{X_n|X_{n-1}}$ and the initial density $f_{X_0}$. The measurement process is determined by the density $f_{Y_n|X_n}$. These two submodels determine the full POMP model, i.e., the joint distribution $f_{X_{0:N},Y_{1:N}}$. If we have a sequence of measurements, $y^*_{n}$, made at times $t_n$, $n=1,\dots,N$, then we think of these data, collectively, as a single realization of the $Y$ process.

The following schematic shows shows causal relations among the process model, the measurement model, and the data. From the statistical point of view, the key perspective is that the model is, at least hypothetically, the process that generated the data.

Structure of a POMP. Arrows show the direction of causality. The closed loop from the state process to itself indicates the dynamic nature of this Markovian process. Information flow runs in the opposite direction.

Mathematically, a POMP is characterized by two conditions.

1. The state process, $X_n$, is Markovian, i.e., $$\mathrm{Prob}[X_n|X_0,\dots,X_{n-1},Y_1,\dots,Y_{n-1}]=\mathrm{Prob}[X_n|X_{n-1}].$$
2. The measurements, $Y_n$, depend only on the state at that time: $$\mathrm{Prob}[Y_n|X_0,\dots,X_{n},Y_1,\dots,Y_{n-1}]=\mathrm{Prob}[Y_n|X_{n}],$$ for all $n=1,\dots,N$.

These conditions can be represented schematically by the following diagram, which indicates the direct dependencies among model variables.

Conditional independence graph of a POMP. The latent state $X_n$ at time $t_n$ is conditionally independent of its history given $X_{n-1}$. The observation $Y_n$ is conditionally independent of all other variables given $X_n$. The underlying time scale can be taken to be either discrete or continuous and the observation times need not be equally spaced.

Basic model components

To implement a POMP model in pomp, we have to specify the measurement and process distributions. Note however, that, for each of the process and the measurement models there are two distinct operations that we might desire to perform:

1. we might wish to simulate, i.e., to draw a random sample from the distribution, or
2. we might wish to evaluate the density itself at given values of $X_n$ and/or $Y_n$.

Following the R convention, we refer to the simulation of $f_{X_n|X_{n-1}}$ as the rprocess component of the POMP model and the evaluation of $f_{X_n|X_{n-1}}$ as the dprocess component. Similarly, the simulation of $f_{Y_n|X_n}$ is the rmeasure component while the evaluation of $f_{Y_n|X_n}$ is the dmeasure component. Finally, we’ll call a simulator of $f_{X_0}$ the rinit component. Collectively, we’ll refer to these, and other, similarly basic elements of the model, as the model’s basic components.

The plug-and-play property

A method that makes no use of the dprocess component is said to be “plug-and-play” or to have the “plug-and-play property”. At present, pomp is focused on such methods, so there is no reason to discuss the dprocess component further in this document. In the following, we will illustrate and explain how one specifies the rprocess, rmeasure, and dmeasure components of a model in pomp. We will illustrate this using some simple examples.

Installing the package

To get started, we must install pomp, if it is not already installed. This package cannot yet be downloaded from CRAN (though that will change in the near future). However, the latest version is always available at the package homepage on Github. See the package website for installation instructions.

Important information for Windows and Mac users.

In this document, we will ultimately learn to implement POMP models using the package’s “C snippet” facility. This allows the user to write model components using snippets of C code, which is then compiled and linked into a running R session. This typically affords a manyfold increase in computation time. It is possible to avoid C snippets entirely by writing model components as R functions, and we will begin by doing so, but the speed-ups afforded by C snippets are typically too good to pass up. To use C snippets, you must be able to compile C codes. Compilers are not by default installed on Windows or Mac systems, so users of such systems must do a bit more work to make use of pomp’s facilities. The installation instructions on the package website give details.

The latent state process

Let us see how to implement a very simple POMP model. In particular, let’s begin by focusing on the famous Ricker model (Ricker 1954), which posits a nonlinear relationship between the size, $N(t)$, of a population in year $t$ and its size, $N(t+1)$, one year later: $$N(t+1)=r\,N(t)\,\exp\left(1-\frac{N(t)}{K}\right).\tag{1}$$ Here, $r$ and $K$ are constant parameters, usually termed the intrinsic growth rate and the carrying capacity, respectively. As written, this is a deterministic model: it does not allow for any variability in the population dynamics. Let’s modify the Ricker equation by assuming that $r$ is not constant, but instead has random variation from year to year. If we assume that the intrinsic growth rate varies from year to year as a lognormal random variable, we obtain $$N(t+1)=r\,N(t)\,\exp\left(1-\frac{N(t)}{K}+\varepsilon(t)\right),\tag{2}$$ where $\varepsilon(t)\sim\mathrm{Normal}(0,\sigma)$. Note that we’ve introduced a new parameter, $\sigma$, which quantifies the intensity of the noise in the population dynamics. Ecologically speaking, Eq. 2 is a model with environmental stochasticity.

Typically, it is relatively straightforward to simulate a POMP model. To accomplish this in pomp, as we’ve already discussed, we specify the rprocess component of the model. We’ll also need to choose values for the model parameters, $r$, $K$, and $\sigma$. We’ll also need to make a choice regarding the initial condition, $N(0)$. The simplest choice is to treat $N(0)=N_0$ as a parameter.

library(pomp)

simulate(t0=0, times=1:20,
  params=c(r=1.2,K=200,sigma=0.1,N_0=50),
  rinit=function (N_0, ...) {
    c(N=N_0)
  },
  rprocess=discrete_time(
    function (N, r, K, sigma, ...) {
      eps <- rnorm(n=1,mean=0,sd=sigma)
      c(N=r*N*exp(1-N/K+eps))
    },
    delta.t=1
  )
) -> sim1

## Warning: 'rmeasure' unspecified: NAs generated.

Notice that we’ve specified the rinit and rprocess components of the model as functions. These functions take as arguments the relevant variables (whether these are state variables or parameters). Importantly, they return named numeric vectors. Names of variables and parameters are very important in pomp. Notice too that we’ve used the discrete_time function, which encodes the fact that our Ricker model is a discrete-time stochastic process (a Markov chain). The first argument of discrete_time is an R function encoding Eq. 2; the second argument specifies the discrete time-step.

Note also that the parameters are furnished in the form of a named vector, and that we’ve specified both t0 and times. The former is the initial time, $t_0$, i.e., the time at which the initial conditions apply. Since our initial condition is that $N(0)=N_0$, our initial time is $t_0=0$. The times argument specifies the observation times $t_1,\dots,t_N$.

Finally, note that we received a warning about NA values being generated. We will soon see what this is about.

What sort of an object is sim1? If we print it, we see

sim1

## <object of class 'pomp'>

sim1 is evidently an object of class ‘pomp’. We refer to these as ‘pomp’ objects.

To get more insight into the structure of sim1, we can use spy:

spy(sim1)

pomp provides methods for plotting ‘pomp’ objects. For example,

plot(sim1)

We can also recast a ‘pomp’ object as a data frame:

as(sim1,"data.frame")

Casting the ‘pomp’ object as a data frame allows us to use ggplot2 graphics:

ggplot(data=as.data.frame(sim1),aes(x=time,y=N))+
  geom_line()

Tp return to the warning we got when we ran simulate: it was telling us that simulate could not make a random draw from the measurement process because we had not supplied it with any information about this process. In particular, we had not supplied a measurement-model simulator. Let’s now see how to specify the measurement-model simulator, or rmeasure.

The measurement model

Let’s suppose that non-lethal traps are used to sample the population to determine its size. Each year, some number of traps are set out and $Y_t$ is the number of animals captured. We might posit $$Y_t\;\sim\;\mathrm{Poisson}(b\,N(t)),$$ where the parameter $b$ is proportional to sampling effort, e.g., the number of traps. This is a measurement model, and we can implement a simulator for it by specifying another function:

simulate(t0=0, times=1:20,
  params=c(r=1.2,K=200,sigma=0.1,N_0=50,b=0.05),
  rinit=function (N_0, ...) {
    c(N=N_0)
  },
  rprocess=discrete_time(
    function (N, r, K, sigma, ...) {
      eps <- rnorm(n=1,mean=0,sd=sigma)
      c(N=r*N*exp(1-N/K+eps))
    },
    delta.t=1
  ),
  rmeasure=function (N, b, ...) {
    c(Y=rpois(n=1,lambda=b*N))
  }
) -> sim2

Note that, again, the rmeasure function need take only the necessary arguments (and ...) and must return a named numeric vector.

Now, in the preceding chunk of code where we construct sim2, there was some redundancy with our earlier construction of sim1. In particular, we specified the same values of t0, times, rinit, and rprocess as before. Since these specifications were stored in sim1, however, we could have simply added in just the new pieces. For example,

simulate(
  sim1,
  params=c(r=1.2,K=200,sigma=0.1,N_0=50,b=0.05),
  rmeasure=function (N, b, ...) {
    c(Y=rpois(n=1,lambda=b*N))
  }
) -> sim2

As before, we can examine our handiwork:

spy(sim2)

as(sim2,"data.frame")

plot(sim2)

ggplot(data=pivot_longer(as(sim2,"data.frame"),-time),
  aes(x=time,y=value,color=name))+
  geom_line()

Notice that now our simulate call produces samples from both the $N$ and the $Y$ distributions. simulate will try to sample from the joint distribution of latent states and observables, but will sample from just the latent state process if the rmeasure component is undefined.

Multiple simulations

We can run multiple simulations using the same parameters:

simulate(sim2,nsim=20) -> sims

ggplot(data=pivot_longer(
         as.data.frame(sims),
         c(Y,N)
       ),
  aes(x=time,y=value,color=name,
    group=interaction(.L1,name)))+
  geom_line()+
  facet_grid(name~.,scales="free_y")+
  labs(y="",color="")

We can also run a single simulation from multiple parameters:

p <- parmat(coef(sim2),3)
p["sigma",] <- c(0.05,0.25,1)
colnames(p) <- LETTERS[1:3]

simulate(sim2,params=p,format="data.frame") -> sims
sims <- pivot_longer(sims,c(Y,N))
ggplot(data=sims,aes(x=time,y=value,color=name,
  group=interaction(.id,name)))+
  geom_line()+
  scale_y_log10()+
  expand_limits(y=1)+
  facet_grid(name~.id,scales="free_y")+
  labs(y="",color="")

In the above, coef extracts the parameter vector stored in sim2. parmat makes a matrix with columns that are copies this vector. The second line modifies this matrix so that the columns are a set of points in parameter space that lie along a line with increasing values of the $\sigma$ parameter.

We can even run multiple simulations at each one of a set of parameters:

simulate(sim2,params=p,
  times=seq(0,3),
  nsim=500,format="data.frame") -> sims
ggplot(data=separate(sims,.id,c("parset","rep")),
  aes(x=N,fill=parset,group=parset,color=parset))+
  geom_density(alpha=0.5)+
  # geom_histogram(aes(y=..density..),position="dodge")+
  facet_grid(time~.,labeller=label_both,scales="free_y")+
  lims(x=c(NA,1000))

Exercise

Modify the sim2 object constructed above to change the measurement model. In particular, assume that $$Y_t\;\sim\;\mathrm{NegBin}(b\,N(t),\theta),$$ where $b$ and $\theta$ are parameters. (By this notation, we understand that $Y_t$ is negative binomially distributed with mean $b\,N(t)$ and variance $b\,N(t)+(b\,N(t))^2/\theta$.) Make and plot some simulations for $t=20,\dots,50$ with different values of the $\theta$ parameter. [Hint: use the rnbinom function with the mu, size parameterization.]

Exercise

The following codes produce several simulations for parameters $r=0.5$, $K=2000$, $\sigma=0.1$ and $b=0.1$, and plot them on the same axes as the data. Notice that the format="data.frame" and include.data=TRUE options facilitate this. Vary the parameters to try to achieve a better fit to the data, as judged purely “by eye”.

vp |>
  simulate(
    params=c(r=0.5,K=2000,sigma=0.1,b=0.1,N_0=2000),
    format="data.frame", include.data=TRUE, nsim=5) |>
  mutate(ds=case_when(.id=="data"~"data",TRUE~"simulation")) |>
  ggplot(aes(x=year,y=pop,group=.id,color=ds))+
  geom_line()+
  labs(color="")

Exercise

Using the ‘pomp’ objects just constructed, explore model simulations at a variety of different parameters. As before, plot simulations and data on the same axes for comparison purposes.

Exercise

Write a C snippet implementing the negative binomial measurement model you explored previously. Fold it into the ‘pomp’ objects just constructed. Remember, there is no need to re-specify components you have already specified: by calling pomp on a ‘pomp’ object you can modify some or all of the basic model components. Plot simulations and data on the same axes, as in the immediately preceding exercise.

Exercise

Write the dmeasure component for your negative binomial model both as an R function and as a C snippet.

Exercise

Compute the likelihood for the parameters you found in your attempt to estimate parameters “by eye”.

A discrete-time deterministic skeleton

A deterministic skeleton of the stochastic Ricker model is the Ricker map, Eq. 1. We implement this for pomp, again either as an R function or a C snippet and pass it to pomp functions via the skeleton argument. For example:

rickC |>
  pomp(
    skeleton=map(
      Csnippet("DN = r*N*exp(1-N/K);"),
      delta.t=1
    ),
    paramnames=c("r","K"), statenames="N"
  ) -> rickC

Here, the left-hand side of Eq. 1 is indicated by the D prefix: in skeleton snippets, we don’t over-write the state variable N. We indicate that the skeleton is a discrete-time map using the map function.

A continuous-time deterministic skeleton

A deterministic skeleton of the Verhulst-Pearl model is the vectorfield (or ordinary differential equation), $$\frac{dN}{dt} = r\,N\,\left(1-\frac{N}{K}\right).$$ We fold this into the vpC ‘pomp’ object so:

vpC |>
  pomp(
    skeleton=vectorfield(Csnippet("DN = r*N*(1-N/K);")),
    paramnames=c("r","K"), statenames="N"
  ) -> vpC

Since the skeleton here is a vectorfield, in this C snippet, DN is filled with the value of the time-derivative of N. See the package help (?Csnippet) for a complete set of rules for writing C snippets.

Trajectories of the deterministic skeleton

With the deterministic skeleton in place we can generate trajectories of the skeleton using trajectory. For example:

p <- parmat(c(r=0.5,K=2000,sigma=0.1,b=0.1,N_0=2000),10)
p["N_0",] <- seq(10,3000,length=10)
vpC |>
  trajectory(params=p,format="data.frame") |>
  ggplot(mapping=aes(x=year,y=N,color=.id,group=.id))+
  guides(color="none")+
  geom_line()+
  theme_bw()

As with simulate, one can use trajectory to compute multiple trajectories at once, for varying values of the parameters.

Parameter estimation using trajectory matching

In pomp, the function traj_objfun constructs an objective function quantifying the mismatch between model predictions and data. For this purpose, it uses the dmeasure component of the model. This function can be given to any of the large variety of numerical optimizers available in R and R packages. These optimizers search parameter space to find parameters under which the likelihood of the data, given a trajectory of the deterministic skeleton, is maximized.

We’ll demonstrate using the Verhulst-Pearl model.

vpC |>
  traj_objfun(
    est=c("K","N_0"),
    params=c(r=0.5,K=2000,sigma=0.1,b=0.1,N_0=2000),
    dmeasure=dmeas, statenames="N", paramnames="b"
  ) -> ofun

This invocation of traj_objfun creates an objective function, ofun, that can be used to estimate the three parameters $K$ and $N_0$. It will hold the other three parameters, $r$, $\sigma$, and $b$, fixed at the values they are given in params.

Notice that, in this code chunk, we had to specify dmeasure once again. Why? What would have happened had we not done so?

What kind of object is ofun?

ofun

## <object of class 'traj_match_objfun'>

‘traj_match_objfun’ objects, like the objective functions created by the pomp functions spect_objfun, probe_objfun, and nlf_objfun, is an R function, but it is also more than an R function. It contains not only the vpC ‘pomp’ object, but it additionally saves information each time it is evaluated. Let’s see how such stateful objective functions make it easy to use a wide range of numerical optimization routines.

We can evaluate ofun at any point in the 2-dimensional $K$-$N_0$ space. For example:

ofun(c(1000,3000))

## [1] 1161.543

The value returned by ofun is the negative log likelihood, as returned by the model’s dmeasure component. We can estimate the parameters using, for example, the subplex algorithm implemented in the subplex package:

library(subplex)

subplex(c(2000,1500),fn=ofun) -> fit
fit 

## $par
## [1] 1968.474 1895.245
## 
## $value
## [1] 276.2893
## 
## $counts
## [1] 318
## 
## $convergence
## [1] 0
## 
## $message
## [1] "success! tolerance satisfied"
## 
## $hessian
## NULL

Note that fit contains, among other things, the estimated parameters (element par) and the minimized value of the negative log likelihood (value).

To make absolutely certain that ofun remembers these estimates, we evaluate it once at fit$par:

ofun(fit$par)

## [1] 276.2893

coef(ofun)

##        r        K    sigma        b      N_0 
##    0.500 1968.474    0.100    0.100 1895.245

logLik(ofun)

## [1] -276.2893

Then we can, for example, extract the fitted trajectory thus:

ofun |>
  trajectory(format="data.frame") |>
  ggplot(mapping=aes(x=year,y=N,color=.id,group=.id))+
  guides(color="none")+
  geom_line()+
  theme_bw()

We can superimpose the model predictions on the data as follows.

ofun |>
  trajectory(format="data.frame") |>
  mutate(
    pop=coef(ofun,"b")*N,
    .id="prediction"
  ) |>
  select(-N) |>
  rbind(
    parus |>
      mutate(
        .id="data",
        pop=as.double(pop)
      )
  ) |>
  pivot_wider(names_from=.id,values_from=pop) |>
  ggplot(aes(x=year))+
  geom_line(aes(y=prediction))+
  geom_point(aes(y=data))+
  expand_limits(y=0)+
  labs(y="pop")

Exercise

Estimate the parameters for the Verhulst-Pearl model assuming negative binomial errors. Do not attempt, at first, to estimate all parameters simultaneously. Focus on estimating $K$, $N_0$ and the parameters of the measurement model. It will probably be helpful to make use of parameter transformations to enforce the model constraints.

Exercise

Estimate some or all of the parameters of the Ricker model using trajectory matching. It is probably a good idea to use parameter transformations.

Local vs global search

The search for the maximum likelihood parameter estimates (MLE) is inherently a global one: the likelihood surface can (and often enough does) possess multiple local maxima, any of which might represent a plausible explanation of the data in hand. IF2, however, is a local search algorithm: essentially, it starts at a user-provided “guess” and attempts to move “uphill” toward the nearest local maximum it can find. To solve the global problem, therefore, it is advisable to run IF2 multiple times, starting from different “guesses”. By studying the results of such a calculation, one hopes to obtain not only the global MLE but, more importantly, an understanding of the shape of the likelihood surface itself.

To illustrate this, let us attempt to determine at which values of the four parameters $r$, $K$, $\sigma$, and $N_0$ our Verhulst-Pearl model best explains the Parus major data. Because this is only a toy problem, we will restrict ourselves to the assumption that $b=1$, despite the fact that it might plausibly be otherwise. In a real scientific study, of course, we would be much less cavalier about fixing parameters arbitrarily in this way.

We begin by constructing a large number of “guesses”. Our design is to distribute these around some multidimensional “box” that, we hope, contains the MLE. Importantly, although our hope is that the box we construct here embraces the MLE, our eventual success does not necessarily depend on it doing so: the IF2 algorithm will be free to search parameter space outside the box. We use the sobol_design function to construct an overdispersed sequence of points within the selected box.

sobol_design(
  lower=c(r=0,K=100,sigma=0,N_0=150,b=1),
  upper=c(r=5,K=600,sigma=2,N_0=150,b=1),
  nseq=100
) -> guesses

plot(guesses,pch=16)

In the above call, lower and upper give the bounds of the box and nseq specifies the number of guesses we desire. Now, from each of these starting guesses, we will run the IF2 algorithm (implemented as mif2 in pomp). As its name suggests, this algorithm works by iteratively applying the particle filter: Crucially, IF2 adds an extra ingredient to the particle filter. As the filter is applied, IF2 adds random perturbations to the parameters. This serves several purposes, including smoothing the likelihood surface, combating particle depletion, and allowing IF2 to obtain information about the local structure of the surface. However, this augmented stochastic variability represents an alteration of the model. IF2 therefore gradually reduces the intensity of the random perturbations.

Local search

Let us examine a single call to mif2.

vpC |>
  mif2(
    params=guesses[1,],
    Np=1000,
    Nmif=20,
    dmeasure=dmeas,
    partrans=parameter_trans(log=c("r","K","sigma","N_0")),
    rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02,N_0=ivp(0.02)),
    cooling.fraction.50=0.5,
    paramnames=c("r","K","sigma","N_0","b"),
    statenames=c("N")
  ) -> vpM

In this call, we specify the starting point with params: here, we’ve just taken this to be the first of the guesses. As in pfilter, Np sets the number of particles to be used. Nmif fixes the number of IF2 iterations that will be performed. We furnish the necessary dmeasure component via the dmeasure argument, and provide parameter transformations via partrans to enforce positivity of the model parameters. The rw.sd argument specifies the intensity of the perturbations that will be applied at each observation time. The perturbations are normally distributed on the transformed (here, the log) scale, so these represent roughly 2% perturbations per observation. The ivp syntax indicates that N_0 is an initial value parameter and accordingly should be treated slightly differently than the regular parameters. The cooling.fraction.50 argument controls the speed at which the magnitude of the perturbations is drawn down. Finally, the paramnames and statenames arguments are needed since dmeas and the parameter transformations are built using C snippets.

What kind of object is vpM?

vpM

## <object of class 'mif2d_pomp'>

pomp provides various methods for ‘mif2d_pomp’ objects. The plot method, for example, produces a diagnostic plot.

plot(vpM)

The first plot (“filter diagnostics”) shows the results of the last iteration of the particle filter, comparable to the plot produced by calling plot on a ‘pfilterd_pomp’ object. The second plot (“convergence diagnostics”) summarizes the results of the sequence of IF2 iterations. We observe that the log likelihood has increased, concomitant with movement in the $r$, $K$, and $\sigma$ parameters and, to a lesser extent, the $N_0$ parameter. We observe that, as intended, the $b$ parameter has remained fixed.

Estimating the log likelihood

Now, although we have observed the intended improvement in the log likelihood, we should be careful to note that the log likelihood displayed in this plot is the log likelihood of the perturbed model. This model differs from the one we are interested in. To compute the likelihood of our focal model at the parameter returned by mif2, we need to perform a few particle filter operations:

vpM |> pfilter() |> logLik() |> replicate(n=5) |> logmeanexp(se=TRUE,ess=TRUE)

##          est           se          ess 
## -143.9880041    0.1035849    4.7923226

Modern computers have multiple processors (cores). To take advantage of these, we can parallelize the above particle-filter computation using the circumstance and doFuture packages:

library(circumstance)
library(doFuture)
plan(multicore)

vpM |> pfilter(Nrep=5) |> logLik() |> logmeanexp(se=TRUE,ess=TRUE)

##          est           se          ess 
## -144.0370140    0.1467573    4.6030137

In the above, the pfilter function is provided by circumstance. It causes the particle filters to be run in parallel, each using a different core. Internally, circumstance uses the doFuture package to parallelize computations. doFuture provides a number of other ways of parallelizing computations, inclusing multisession (needed on Windows computers) and cluster (for parallelizing across multiple machines). These are set using the plan() function. doFuture also provides parallel random-number generators, so that we can consider each of the parallel computations to be independent.

Global search

Now we repeat this calculation from each of our 100 guesses. Because these computations are each independent of the others, it is both easy and worthwhile to parallelize them. The circumstance package is also useful for this purpose. In particular, it defines a mif2 method that takes a data frame of guesses of the kind we defined above. For each of these, it performs an iterated filtering computation.

vpM |> mif2(starts=guesses) -> mifs

Notice that, in the above, the settings of mif2 are not explicitly given here. They are inherited from vpM, which we computed previously. If we had wanted to modify any of them (e.g., Nmif, Np, rw.sd, etc.), we would have simply done so in the call to mif2. For example, to increase the number of particles and the number of iterations, we would do:

vpM |>
  mif2(
    starts=guesses,
    Nmif=30,
    Np=2000
  ) -> mifs

Let us examine the results of our mif2 computation. The following code chunk extracts the traces of the IF2 calculations. These are the movements of the parameters and the likelihood plotted against filter iteration.

mifs |>
  traces() |>
  melt() |>
  filter(name!="b") |>
  ggplot(aes(x=iteration,y=value,group=.L1,color=.L1))+
  geom_line()+
  facet_wrap(~name,scales="free_y")+
  guides(color="none")

More simply, executing plot(mifs), would give us a less fancy plot.

Next, we examine the diagnostics produced in the last iteration of the filter.

mifs |> 
  as("data.frame") |> 
  pivot_longer(-c(.L1,year)) |>
  ggplot(aes(x=year,y=value,group=.L1,color=.L1))+
  geom_line()+
  facet_wrap(~name,scales="free_y",ncol=1)+
  guides(color="none")

Next, to estimate the log likelihoods, we run a few independent particle filters on the results of our mif2 computations. We combine the resulting log likelihood estimates with the parameter estimates themselves (obtained using coef()).

mifs |>
  pfilter(Nrep=5) |>
  logLik() |>
  melt() |>
  separate(name,into=c(".id","rep")) |>
  group_by(.id) |>
  reframe(melt(logmeanexp(value,se=TRUE))) |>
  ungroup() |>
  bind_rows(
    mifs |> coef() |> melt()
  ) |>
  pivot_wider() |>
  rename(
    loglik=est,
    loglik.se=se
  ) -> estimates

To display the results of such a multi-start search pairwise scatterplot matrices are useful:

estimates |>
  bind_rows(guesses) |>
  filter(is.na(loglik) | loglik>max(loglik,na.rm=TRUE)-30) |>
  mutate(col=if_else(is.na(loglik),"#99999955","#ff0000ff")) |>
  {
    \(dat) pairs(~loglik+r+sigma+K+N_0,data=dat,col=dat$col,pch=16)
  }()

More search effort

At this point, there’s no particular reason to suspect that our IF2 searches have arrived at their destination. In general, it is hard to know a priori how much effort will be required to find the MLE. Let’s continue the search, starting with the best points we’ve uncovered so far.

estimates |>
  filter(!is.na(loglik)) |>
  filter(loglik > max(loglik)-30) |>
  select(-.id,-loglik,-loglik.se) -> starts

vpM |>
  mif2(starts=starts) |>
  mif2() -> mf

mf |>
  pfilter(Nrep=5) |>
  logLik() |>
  melt() |>
  separate(name,into=c(".id","rep")) |>
  group_by(.id) |>
  reframe(melt(logmeanexp(value,se=TRUE))) |>
  ungroup() |>
  bind_rows(
    mf |> coef() |> melt()
  ) |>
  pivot_wider() |>
  rename(
    loglik=est,
    loglik.se=se
  ) -> ests1

Note that, in the above, we’ve performed a total of 100 mif2 iterations per starting point. The code above also does the post-mif2 likelihood estimation. It returns just the parameter and likelihood estimates.

We can combine the new estimates with the old ones into a general database:

estimates |> bind_rows(ests1) |> select(-.id) -> estimates

estimates |> arrange(-loglik) |> head()

estimates |>
  filter(loglik>max(loglik,na.rm=TRUE)-4) |>
  select(loglik,r,sigma,K,N_0) |>
  plot_matrix()

In this plot, we begin to see the emergence of structure in the likelihood surface. In particular, what looks like a ridge of high likelihood is visible in the $r$-$\sigma$ projection. Such structures are very interesting in that they contain clues as to the manner in which the model is fitting the data. They can also pose challenges to efficient estimation, since climbing up to a ridge is harder than traversing it.

Likelihood profile

We can improve the quality of our estimates and obtain likelihood-ratio-test-based confidence intervals by constructing profile likelihoods. In a likelihood profile, one varies the focal parameter (or parameters) across some range, maximizing the likelihood over the remaining parameters at each value of the focal parameter. The following codes construct a likelihood profile over $r$.

estimates |>
  filter(loglik>max(loglik)-10) |>
  select(r,K,sigma,N_0,b) |>
  apply(2,range) -> ranges
ranges

##                r        K    sigma      N_0 b
## [1,]  0.05091814 202.5526 0.291154 139.8762 1
## [2,] 15.24609275 623.1927 1.300350 162.2693 1

profile_design(
  r=10^seq(
    from=log10(ranges[1,1]),
    to=log10(ranges[2,1]),
    length=20
  ),
  lower=ranges[1,-1],
  upper=ranges[2,-1],
  nprof=50
) -> starts

dim(starts)

## [1] 1000    5

starts |>
  select(r,sigma,K,N_0,b) |>
  plot_matrix()

vpM |>
  mif2(
    starts=starts,
    partrans=parameter_trans(log=c("K","sigma","N_0")),
    rw.sd=rw_sd(K=0.02,sigma=0.02,N_0=ivp(0.02)),
    paramnames=c("K","sigma","N_0","b")
  ) |>
  mif2() -> mf

mf |>
  pfilter(Nrep=5) |>
  logLik() |>
  melt() |>
  separate(name,into=c(".id","rep")) |>
  group_by(.id) |>
  reframe(melt(logmeanexp(value,se=TRUE))) |>
  ungroup() |>
  bind_rows(
    mf |> coef() |> melt()
  ) |>
  pivot_wider() |>
  select(-.id) |>
  rename(
    loglik=est,
    loglik.se=se
  ) -> r_prof

Notice that we’ve changed the mif2 perturbations (rw.sd): we’ve removed the perturbation on the $r$ parameter, since we want to hold this parameter fixed.

We add these points to our database:

estimates |>
  bind_rows(r_prof) -> estimates

Next, we plot the likelihood profile. The following plot shows the top two estimates for each value of $r$ with error bars showing ±2 s.e. and a loess smooth.

r_prof |>
  group_by(r) |>
  filter(rank(-loglik)<=2) |>
  ungroup() |>
  ggplot(aes(x=r,y=loglik,
    ymin=loglik-2*loglik.se,ymax=loglik+2*loglik.se))+
  geom_point()+
  geom_errorbar()+
  geom_smooth(method="loess",span=0.2)+
  scale_x_log10()

We see that the 95% likelihood ratio test confidence interval (CI) appears to be one-sided: the CI is roughly $r>0.73$.

If we plot the profile trace of $\sigma$, we see that, as we increase $r$, we have to increase the intensity of the environmental stochasticity to maintain a good fit. Why is this?

r_prof |>
  group_by(r) |>
  filter(rank(-loglik)<=2) |>
  ungroup() |>
  ggplot(aes(x=r,y=sigma))+
  geom_point()+
  geom_smooth(method="loess",span=0.2)+
  scale_x_log10()+
  labs(y=expression(sigma))

Running multiple pmcmc chains

If we seek to do full-information Bayesian inference, we can use particle Markov chain Monte Carlo, implemented in pomp in the pmcmc function. The following codes cause parallel pMCMC chains to be run, each beginning at one of the estimates obtained using mif2. Note that we add a prior by furnishing a prior probability density evaluation function (dprior=...). In this case, we assume a product prior: marginally uniform in each of the parameters $r$, $\sigma$, $K$, and $N_0$.

r_prof |>
  group_by(r) |>
  filter(loglik==max(loglik)) |>
  ungroup() |>
  filter(
    r > 0.75,
    r < 6,
    sigma < 2,
    K > 100,
    K < 600,
    N_0 > 100,
    N_0 < 600
  ) |>
  select(-loglik,-loglik.se) -> starts
starts

vpM |>
  pmcmc(
    starts=starts,
    Nmcmc=2000,Np=200,
    dprior=Csnippet(r"{
      lik = dunif(r,0,10,1)+dunif(sigma,0,2,1)+
            dunif(K,0,600,1)+dunif(N_0,0,600,1);
      lik = (give_log) ? lik : exp(lik);}"),
    paramnames=c("K","N_0","r","sigma"),
    proposal=mvn_rw_adaptive(
      rw.sd=c(r=0.02,sigma=0.02,K=50,N_0=50),
      scale.start=100,shape.start=100
    )
  ) -> chains

chains |>
  pmcmc(
    Nmcmc=20000,
    proposal=mvn_rw(covmat(chains))
  ) -> chains

In the above, we first subset out those mif2 estimates that are consistent with our prior. At each of these, we then perform a number of MCMC moves, using an adaptive multivariate-normal random-walk proposal distribution mvn.rw.adaptive. To complete the inference, we do more MCMC iterations using a multivariate random-walk proposal (mvn.rw) with covariance matrix derived from the preceding computation (covmat).

Convergence diagnostics

Now we investigate the chains, to try and determine whether they have converged.

library(coda)
chains |> traces() -> traces
rejectionRate(traces[,c("r","sigma","K","N_0")])

##         r     sigma         K       N_0 
## 0.8169429 0.8169429 0.8169429 0.8169429

We see that the rejection rate is very good. Let us examine the autocorrelation in the chains.

traces |> autocorr.diag(lags=c(1,5,10,50,100))

##              loglik log.prior           r          K       sigma          N_0
## Lag 1   0.926929299       NaN 0.915944167 0.91002475  0.91642728  0.917472089
## Lag 5   0.708801155       NaN 0.662951832 0.64140619  0.66051346  0.657948453
## Lag 10  0.521030101       NaN 0.461219574 0.43428713  0.45264129  0.436384889
## Lag 50  0.066046657       NaN 0.046472319 0.06245856  0.02921802  0.027336265
## Lag 100 0.004874857       NaN 0.009990407 0.02272330 -0.01255195 -0.001893778
##           b
## Lag 1   NaN
## Lag 5   NaN
## Lag 10  NaN
## Lag 50  NaN
## Lag 100 NaN

traces |> effectiveSize()

##    loglik log.prior         r         K     sigma       N_0         b 
##  4399.400     0.000  5427.223  5810.303  5597.652  5865.711     0.000

traces <- window(traces,thin=100,start=2000)
traces |> effectiveSize()

##    loglik log.prior         r         K     sigma       N_0         b 
##  1251.555     0.000  1420.696  1512.574  1328.394  1203.467     0.000

The autocorrelation is strong, but drops to small values by lag 100. Accordingly, we thin the chains by a factor of 100. We discard 2000 iterations as a burn-in.

Now let us examine the traces.

traces |>
  lapply(as.data.frame) |>
  lapply(rowid_to_column,"iter") |>
  bind_rows(.id="chain") |>
  select(chain,iter,loglik,r,sigma,K,N_0) |>
  pivot_longer(c(-chain,-iter)) |>
  ggplot(aes(x=iter,group=chain,color=chain,y=value))+
  guides(color="none")+
  labs(x="iteration",y="")+
  geom_line(alpha=0.3)+geom_smooth(method="loess",se=FALSE)+
  facet_wrap(name~.,scales="free_y",strip.position="left",ncol=2)+
  theme(
    strip.placement="outside",
    strip.background=element_rect(fill=NA,color=NA)
  )

gelman.diag(traces[,c("r","sigma","K","N_0")])

## Potential scale reduction factors:
## 
##       Point est. Upper C.I.
## r              1       1.01
## sigma          1       1.01
## K              1       1.00
## N_0            1       1.01
## 
## Multivariate psrf
## 
## 1.01

The trace-plots show good mixing and the Gelman-Rubin statistic confirms that the chains are mixing among themselves: a good indicator of convergence. Thus, insofar as we can tell by these diagnostics, the MCMC iterations appear to have converged to their stationary distribution and we are sampling it well after thinning.

Examining the posterior density

The plots of the marginal posterior densities are shown below.

traces |>
  lapply(as.data.frame) |>
  lapply(rowid_to_column,"iter") |>
  bind_rows(.id="chain") |>
  select(chain,iter,loglik,r,sigma,K,N_0) |>
  pivot_longer(c(-chain,-iter)) |>
  ggplot(aes(x=value))+
  geom_density()+
  geom_rug()+
  labs(x="")+
  facet_wrap(~name,scales="free",strip.position="bottom")+
  theme(
    strip.placement="outside",
    strip.background=element_rect(fill=NA,color=NA)
  )

traces |> summary()

## 
## Iterations = 2000:20000
## Thinning interval = 100 
## Number of chains = 7 
## Sample size per chain = 181 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##               Mean     SD Naive SE Time-series SE
## loglik    -143.036  1.596 0.044848       0.045575
## log.prior  -15.790  0.000 0.000000       0.000000
## r            6.117  2.557 0.071823       0.068898
## K          214.133 13.897 0.390417       0.370018
## sigma        0.866  0.235 0.006603       0.006569
## N_0        149.054 11.985 0.336693       0.345759
## b            1.000  0.000 0.000000       0.000000
## 
## 2. Quantiles for each variable:
## 
##                2.5%       25%       50%      75%    97.5%
## loglik    -147.0794 -143.7932 -142.7763 -141.943 -140.789
## log.prior  -15.7896  -15.7896  -15.7896  -15.790  -15.790
## r            1.3148    4.0021    6.4429    8.351    9.810
## K          191.7318  204.4396  213.0764  221.698  245.346
## sigma        0.4222    0.7029    0.8738    1.027    1.293
## N_0        126.2790  140.7598  148.9529  157.217  172.633
## b            1.0000    1.0000    1.0000    1.000    1.000

traces |>
  lapply(as.data.frame) |>
  lapply(rowid_to_column,"iter") |>
  bind_rows(.id="chain") |>
  select(loglik,r,sigma,K,N_0) |>
  plot_matrix()

Hindcast and smoothing with parametric uncertainty

Let us revisit the matter of hindcasting. In a preceding section, using an ensemble of particle filters, we were able to estimate the likely trajectory that the latent state process followed in generating the data we see, conditional on a given parameter vector. Naturally, since our estimate of the parameters is itself uncertain, our actual uncertainty regarding the hindcast is somewhat larger. We can use particle Markov chain Monte Carlo to fold in the parametric uncertainty. The following codes select one randomly-chosen particle trajectory from each iteration of our pMCMC chains. The chains are thinned and the burn-in period is discarded, of course. Note that it is not necessary to weight the samples by the likelihood: if the pMCMC chains have converged, the samples can be taken to be equally weighted samples from the smoothing distribution.

chains |>
  filter_traj() |>
  melt() |>
  filter(rep > 1000, rep %% 100 == 0) |>
  mutate(year=time(vpM)[time]) |>
  pivot_wider(names_from=name) |>
  group_by(year) |>
  reframe(
    label=c("lo","med","hi"),
    p=c(0.025,0.5,0.975),
    q=wquant(N,probs=p)
  ) |>
  ungroup() |>
  select(-p) |>
  pivot_wider(names_from=label,values_from=q) -> quants2

bind_rows(
  with=quants2,
  without=quants1,
  .id="uncert"
) |>
  ggplot()+
  geom_line(aes(x=year,y=med,color=uncert))+
  geom_ribbon(aes(x=year,ymin=lo,ymax=hi,fill=uncert),color=NA,alpha=0.4)+
  geom_point(data=parus,aes(x=year,y=pop))+
  labs(y="N",fill="parametric\nuncertainty",color="parametric\nuncertainty")

Simulation of the fitted model

Ultimately, since the model is viewed, at least hypothetically, as the process that generated the data, simulation of the fitted model is a central tool we have for model criticism. Let’s plot the data and several simulated realizations of the model process on the same axes.

r_prof |>
  filter(loglik==max(loglik)) -> mle

mlepomp <- as(mifs[[1]],"pomp")
coef(mlepomp) <- mle

mlepomp |>
  simulate(nsim=8,format="data.frame",include.data=TRUE) |>
  ggplot(mapping=aes(x=year,y=pop,group=.id,alpha=(.id=="data")))+
  scale_alpha_manual(values=c(`TRUE`=1,`FALSE`=0.2),
    labels=c(`FALSE`="simulation",`TRUE`="data"))+
  labs(alpha="")+
  geom_line()+
  theme_bw()

The first lines above simply extract the maximum likelihood estimates (mle) from our profle computation. The next pair of lines plug these MLE parameters into a ‘pomp’ object (mlepomp) containing the model and the data. The last set of lines do the simulation and the plotting.

Although it is clear from these plots that the estimated model has more variability and is thus able to explain the data better, it can be hard to read much from spaghetti plots such as this. It’s almost always a good idea to plot the data together with several simulated realizations in order to help assess how similar the two are.

mlepomp |>
  simulate(nsim=11,format="data.frame",include.data=TRUE) |>
  ggplot(mapping=aes(x=year,y=pop,group=.id,color=(.id=="data")))+
  scale_color_manual(values=c(`TRUE`="black",`FALSE`="grey50"),
    labels=c(`FALSE`="simulation",`TRUE`="data"))+
  labs(color="")+
  geom_line()+
  facet_wrap(~.id)+
  theme_bw()

Model checking with probes

Visual comparison of simulations and data is always a good idea. An indication that the data are not a plausible realization of the model is evidence for lack of fit. In particular, if we have any set of summary statistics, or probes, we can apply them to both simulated and actual data. The probe function facilitates this comparison. Let’s perform this operation using several of the summary statistics provided with pomp: we’ll use the mean, several quantiles, and the autocorrelation at lags 1 and 3.

mlepomp |>
  probe(nsim=200,probes=list(
    mean=probe_mean("pop"),
    q=probe_quantile("pop",probs=c(0.05,0.25,0.5,0.75,0.95)),
    probe_acf("pop",lags=c(1,3),type="corr",transform=log)
  )) -> vp_probe

vp_probe

## <object of class 'probed_pomp'>

For ‘probed_pomp’ object, there are summary and plot methods. There is also an as.data.frame method.

summary(vp_probe)

## $coef
##       loglik    loglik.se            r            K        sigma          N_0 
## -141.3529509    0.2174540    4.5902378  209.3591890    0.7347252  148.6349594 
##            b 
##    1.0000000 
## 
## $nsim
## [1] 200
## 
## $quantiles
##   mean   q.5%  q.25%  q.50%  q.75%  q.95% acf[1] acf[3] 
##  0.485  0.555  0.725  0.240  0.310  0.770  0.620  0.800 
## 
## $pvals
##      mean      q.5%     q.25%     q.50%     q.75%     q.95%    acf[1]    acf[3] 
## 0.9751244 0.8955224 0.5572139 0.4875622 0.6268657 0.4676617 0.7661692 0.4079602 
## 
## $synth.loglik
## [1] -19.36771

Evidently, summary returns a list with several elements. The quantiles element contains, for each probe, what fraction of the nsim simulations had probe values below the value of the probe applied to the data. The pvals element contains $P$-values associated with the two-sided test of the hypothesis that the data were generated by the model.

plot(vp_probe)

The plot depicts the multivariate distribution of the probes under the model, with the data-values superimposed. On the diagonal, we see the marginal distributions of the individual probes, represented as histograms, with the vertical line signifying the value of the corresponding probe on the data. Above the diagonal, the scatterplots show the pairwise distributions of probes and the crosshairs, the corresponding data-values. Below the diagonal, the panels contains the pairwise correlations among the simulated probes.