Licensed under the Creative Commons Attribution-NonCommercial license. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC This document has its origins in the SISMID short course on Simulation-based Inference given by Aaron King and Edward Ionides.

Produced in R version 3.6.3 using pomp version 2.8.

Important Note: These materials have been updated for use with version 2.8. As of version 2, pomp syntax has changed substantially. These changes are documented on the pomp website.

Introduction: ecological and epidemiological dynamics

Noisy clockwork: Time series analysis of population fluctuations in animals

Six problems of Bjørnstad and Grenfell (2001)

Obstacles for ecological modeling and inference via nonlinear mechanistic models:

  1. Combining measurement noise and process noise.
  2. Including covariates in mechanistically plausible ways.
  3. Using continuous-time models.
  4. Modeling and estimating interactions in coupled systems.
  5. Dealing with unobserved variables.
  6. Modeling spatial-temporal dynamics.

The same issues arise for epidemiological modeling and inference via nonlinear mechanistic models.

Objectives

  1. To show how stochastic dynamical systems models can be used as scientific instruments.
  2. To give students the ability to formulate models of their own.
  3. To teach efficient approaches for performing scientific inference using POMP models.
  4. To familiarize students with the pomp package.
  5. To give students opportunities to work with such inference methods.
  6. To provide documented examples for student re-use.

Partially observed Markov process (POMP) models

Schematic of the structure of a POMP

showing causal relations.

The key perspective to keep in mind is that the model is to be viewed as the process that generated the data.

The Markov assumption

  • \(\mathbb{P}\left[{X_n|X_0,\dots,X_{n-1}}\right]=\mathbb{P}\left[{X_n|X_{n-1}}\right]\).
  • Interpretation: knowledge of the system’s state at any point in time is sufficient to determine the distribution of possible futures.
  • Alternative interpretation: the system’s state is sufficiently rich so as to encompass all important features of the system’s history
  • Systems with delays can usually be rewritten as Markovian systems, at least approximately.
  • An important special case: any system of differential equations is Markovian.

Notation for partially observed Markov process models

  • Write \(X_n=X(t_n)\) and \(X_{0:N}=(X_0,\dots,X_N)\). Let \(Y_n\) be a random variable modeling the observation at time \(t_n\).

  • The one-step transition density, \(f_{X_n|X_{n-1}}(x_n|x_{n-1};\theta)\), together with the measurement density, \(f_{Y_n|X_n}(y_n|x_n;\theta)\) and the initial density, \(f_{X_0}(x_0;\theta)\), specify the entire joint density via

\[f_{X_{0:N},Y_{1:N}}(x_{0:N},y_{1:N};\theta) = f_{X_0}(x_0;\theta)\,\prod_{n=1}^N\!f_{X_n | X_{n-1}}(x_n|x_{n-1};\theta)\,f_{Y_n|X_n}(y_n|x_n;\theta).\]

  • The marginal density for sequence of measurements, \(Y_{1:N}\), evaluated at the data, \(y_{1:N}^*\), is

\[ f_{Y_{1:N}}(y^*_{1:N};\theta)=\int f_{X_{0:N},Y_{1:N}}(x_{0:N},y^*_{1:N};\theta)\, dx_{0:N}.\]

Another POMP model schematic

showing dependence among model variables:

The state process, \(X_n\), is Markovian, i.e., \[\mathbb{P}\left[{X_n|X_0,\dots,X_{n-1},Y_1,\dots,Y_{n-1}}\right]=\mathbb{P}\left[{X_n|X_{n-1}}\right].\] Moreover, the measurements, \(Y_n\), depend only on the state at that time: \[\mathbb{P}\left[{Y_n|X_0,\dots,X_{n},Y_1,\dots,Y_{n-1}}\right]=\mathbb{P}\left[{Y_n|X_{n}}\right],\] for all \(n=1,\dots,N\).

To think algorithmically, we define some function calls:

  • rprocess( ): a draw from \(f_{X_n|X_{n-1}}(x_n| x_{n-1};\theta)\)
  • dprocess( ): evaluation of \(f_{X_n|X_{n-1}}(x_n| x_{n-1};\theta)\)
  • rmeasure( ): a draw from \(f_{Y_n|X_n}(y_n| x_n;\theta)\)
  • dmeasure( ): evaluation of \(f_{Y_n|X_n}(y_n| x_n;\theta)\)
  • initializer( ): a draw from \(f_{X_0}(x_0;\theta)\)

What does it mean for methodology to be simulation-based?

  • Simulating random processes is often much easier than evaluating their transition probabilities.
  • In other words, we may be able to write rprocess() but not dprocess().
  • Simulation-based methods require the user to specify rprocess() but not dprocess().
  • Plug-and-play, likelihood-free and equation-free are alternative terms for “simulation-based” methods.
  • Much development of simulation-based statistical methodology has occurred in the past decade.

The pomp package for POMP models

The following diagrams show the structure of a POMP model schematically.

Example: the Ricker model

The deterministic Ricker map

The Ricker map describes the deterministic dynamics of a simple population: \[N_{t+1} = r\,N_{t}\,\exp(-N_{t})\] Here, \(N_t\) is the population density at time \(t\) and \(r\) is a fixed value (a parameter), related to the population’s intrinsic capacity to increase. \(N\) is a state variable, \(r\) is a parameter. If we know \(r\) and the initial condition \(N_0\), the deterministic Ricker equation predicts the future population density at all times \(t=1,2,\dots\). We can view the initial condition, \(N_0\) as a special kind of parameter, an initial-value parameter.

Process noise

We can model process noise in this system by making the growth rate \(r\) into a random variable. For example, if we assume that the intrinsic growth rate is log-normally distributed, \(N\) becomes a stochastic process governed by \[N_{t+1} = r\,N_{t}\,\exp(-N_{t}+\varepsilon_{t}), \qquad \varepsilon_{t}\;\sim\;\mathrm{Normal}\left(0,\sigma\right),\] where the new parameter \(\sigma\) is the standard deviation of the noise process \(\varepsilon\).

Measurement error

Let’s suppose that the Ricker model is our model for the dynamics of a real population. However, we cannot know the exact population density at any time, but only estimate it through sampling.

Let’s model measurement error by assuming the measurements, \(y_t\), are Poisson with mean \(\phi\,N_t\): \[y_{t}\;\sim\;\mathrm{Poisson}\left(\phi\,N_{t}\right)\]

In this equation,

  1. \(N_t\) is the true population density at time \(t\),
  2. \(y_t\) is the number of individuals sampled at time \(t\),
  3. the choice of units for \(N\) is peculiar and depends on the parameters (e.g., \(N=\log(r)\) is the equilibrium of the deterministic model),
  4. the parameter \(\phi\) is proportional to our sampling effort, and also has peculiar units.

Working with the Ricker model in pomp.

The R package pomp provides facilities for modeling POMPs, a toolbox of statistical inference methods for analyzing data using POMPs, and a development platform for implmenting new POMP inference methods. The basic data-structure provided by pomp is the object of class pomp, alternatively known as a pomp object. It is a container that holds real or simulated data and a POMP model, possibly together with other information such as model parameters, that may be needed to do things with the model and data.

Let’s see what can be done with a pomp object. First, we’ll load some packages, including pomp.

library(ggplot2)
library(plyr)
library(reshape2)
library(pomp)
stopifnot(packageVersion("pomp")>="2.8")

A pre-built pomp object encoding the Ricker model comes included with the package. Construct it by executing

ricker() -> rick

This has the effect of creating a pomp object named rick in your workspace. We can plot the data by doing

plot(rick)

We can simulate by doing

x <- simulate(rick)

What kind of object have we created?

class(x)
## [1] "pomp"
## attr(,"package")
## [1] "pomp"
plot(x)

Why do we see more time series in the simulated pomp object?

We can turn a pomp object into a data frame:

y <- as.data.frame(rick)
head(y)
##   time   y
## 1    0  68
## 2    1   2
## 3    2  87
## 4    3   0
## 5    4  12
## 6    5 174
head(simulate(rick,format="data.frame"))
##   time .id           N            e   y
## 1    0   1  7.00000000  0.000000000  78
## 2    1   1  0.19117895 -0.400455545   3
## 3    2   1 10.00176418  0.348485839 116
## 4    3   1  0.02194212  0.079655415   0
## 5    4   1  0.84780449 -0.123815830   7
## 6    5   1 16.25242037  0.001151556 145

We can also run multiple simulations simultaneously:

x <- simulate(rick,nsim=10)
class(x)
## [1] "pompList"
## attr(,"package")
## [1] "pomp"
sapply(x,class)
##  [1] "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp" "pomp"
x <- simulate(rick,nsim=10,format="data.frame")
head(x)
##   time .id N e  y
## 1    0   1 7 0 73
## 2    0   2 7 0 59
## 3    0   3 7 0 61
## 4    0   4 7 0 74
## 5    0   5 7 0 62
## 6    0   6 7 0 76
str(x)
## 'data.frame':    510 obs. of  5 variables:
##  $ time: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ .id : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ N   : num  7 7 7 7 7 7 7 7 7 7 ...
##  $ e   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ y   : num  73 59 61 74 62 76 63 58 62 61 ...

Also,

x <- simulate(rick,nsim=9,format="d",include.data=TRUE)
ggplot(data=x,aes(x=time,y=y,group=.id,color=(.id=="data")))+
  geom_line()+guides(color=FALSE)+
  facet_wrap(~.id,ncol=2)

We refer to the deterministic map as the “skeleton” of the stochastic map. We can compute a trajectory of the the deterministic skeleton using trajectory:

y <- trajectory(rick)
dim(y)
## [1]  2  1 51
dimnames(y)
## $variable
## [1] "N" "e"
## 
## $rep
## NULL
## 
## $time
## NULL
plot(time(rick),y["N",1,],type="l")

Notice that rick has parameters associated with it:

coef(rick)
##        r    sigma      phi        c      N_0      e_0 
## 44.70118  0.30000 10.00000  1.00000  7.00000  0.00000

These are the parameters at which the simulations and deterministic trajectory computations above were done. We can run these at different parameters:

theta <- coef(rick)
theta[c("r","N.0")] <- c(5,3)
y <- trajectory(rick,params=theta)
plot(time(rick),y["N",1,],type="l")

x <- simulate(rick,params=theta)
plot(x,var="y")

We can also change the parameters stored inside of rick:

coef(rick,c("r","N.0","sigma")) <- c(39,0.5,1)
coef(rick)
##     r sigma   phi     c   N_0   e_0   N.0 
##  39.0   1.0  10.0   1.0   7.0   0.0   0.5
plot(simulate(rick),var="y")

In all of the above, it’s possible to work with more than one set of parameters at a time. For example:

p <- parmat(coef(rick),500)
dim(p); dimnames(p)
## [1]   7 500
## $variable
## [1] "r"     "sigma" "phi"   "c"     "N_0"   "e_0"   "N.0"  
## 
## $rep
## NULL
p["r",] <- seq(from=2,to=40,length=500)
y <- trajectory(rick,params=p,times=200:1000)
matplot(p["r",],y["N",,],pch=".",col='black',xlab='r',ylab='N',log='x')

How do you interpret the above plot? This is called a *one-parameter bifurcation diagram".

More information on manipulating and extracting information from pomp objects can be viewed in the help pages (methods?pomp).

There are a number of other examples included with the package. Do pompExample() to see a list of these. More examples can be found in the pompExamples package:

library(pompExamples)
pompExample()

Inference algorithms in pomp

pomp provides a wide range of inference algorithms. We’ll learn about these in detail soon, but for now, let’s just look at some of their general features.

The pfilter function runs a simple particle filter. It can be used to evaluate the likelihood at a particular set of parameters. One uses the Np argument to specify the number of particles to use:

pf <- pfilter(rick,Np=1000)
class(pf)
## [1] "pfilterd_pomp"
## attr(,"package")
## [1] "pomp"
plot(pf)

logLik(pf)
## [1] -157.8544

Note that pfilter returns an object of class pfilterd.pomp. This is the general rule: inference algorithms return objects that are pomp objects with additional information. The package provides tools to extract this information. We can run the particle filter again by doing

pf <- pfilter(pf)
logLik(pf)
## [1] -158.2566

which has the result of running the same computation again. Note that, because the particle filter is a Monte Carlo algorithm, we get a slightly different estimate of the log likelihood.

Note that, by default, running pfilter on a pfilterd.pomp object causes the computation to be re-run with the same parameters as before. Any additional arguments we add override these defaults. This is the general rule in pomp. For example,

pf <- pfilter(pf,Np=100)
logLik(pf)
## [1] -158.822

Here, the particle filtering has been performed with only 100 particles.

Building a custom pomp object

A real pomp data analysis begins with constructing one or more pomp objects to hold the data and the model or models under consideration. We’ll illustrate this process a dataset of Parus major abundance in Wytham Wood, near Oxford (McCleery and Perrins 1991).

Download and plot the data:

dat <- read.csv("http://kingaa.github.io/clim-dis/intro/parus.csv")
head(dat)
##   year pop
## 1 1960 148
## 2 1961 258
## 3 1962 185
## 4 1963 170
## 5 1964 267
## 6 1965 239
plot(pop~year,data=dat,type='o')

Let’s suppose that we want to fit the stochastic Ricker model discussed above to these data.

The call to construct a pomp object is, naturally enough, pomp. Documentation on this function can be had by doing ?pomp. Learn about the various things you can do once you have a pomp object by doing methods?pomp and following the links therein. Read an overview of the package as a whole with links to its main features by doing package?pomp. A complete index of the functions in pomp is returned by the command library(help=pomp). Finally, the home page for the pomp project is http://kingaa.github.io/pomp; there you have access to the complete source code, tutorials, manuals, issues page, news blog, etc.

Now, to construct our pomp object:

library(pomp)
parus <- pomp(dat,times="year",t0=1959)

The times argument specifies that the column of dat labelled “year” gives the measurement times; t0 is the “zero-time”, the time at which the state process will be initialized. We’ve set it to one year prior to the beginning of the data. Plot the new pomp object:

plot(parus)

Adding in the process model simulator

We can add the stochastic Ricker model to parus by writing a Csnippet that simulates one realization of the stochastic process, from an arbitary time \(t\) to \(t+1\), given arbitrary states and parameters. We provide this to pomp in the form of a Csnippet, a little snippet of C code that performs the computation. The following does this.

stochStep <- Csnippet("
  e = rnorm(0,sigma);
  N = r*N*exp(-N+e);
")
pomp(parus,rprocess=discrete_time(step.fun=stochStep,delta.t=1),
     paramnames=c("r","sigma"),statenames=c("N","e")) -> parus

Note that in the above, we use the exp and rnorm functions from the R API. In general any C function provided by R is available to you. pomp also provides a number of C functions that are documented in the header file, pomp.h, that is installed with the package. See the Csnippet documentation (?Csnippet) to read more about how to write them. Note too that we use discrete.time.sim here because the model is a stochastic map. We specify that the time step of the discrete-time process is delta.t, here, 1 yr.

At this point, we have what we need to simulate the state process:

sim <- simulate(parus,params=c(N.0=1,e.0=0,r=12,sigma=0.5),format="d")
plot(N~year,data=sim,type='o')

Adding in the measurement model and parameters

We complete the specification of the POMP by specifying the measurement model. To obtain the Poisson measurement model described above, we write two Csnippets. The first simulates:

rmeas <- Csnippet("pop = rpois(phi*N);")

The second computes the likelihood of observing pop birds given a true density of N:

dmeas <- Csnippet("lik = dpois(pop,phi*N,give_log);")

[Note the give_log argument. When this code is evaluated, give_log will be set to 1 if the log likelihood is desired, and 0 else.] We add these into the pomp object:

pomp(parus,rmeasure=rmeas,dmeasure=dmeas,statenames=c("N"),paramnames=c("phi")) -> parus

Now we can simulate the whole POMP. First, let’s add some parameters:

coef(parus) <- c(N.0=1,e.0=0,r=20,sigma=0.1,phi=200)
library(ggplot2)
sims <- simulate(parus,nsim=3,format="d",include.data=TRUE)
ggplot(data=sims,mapping=aes(x=year,y=pop))+geom_line()+
  facet_wrap(~.id,ncol=1,scales="free_y")

Adding in the deterministic skeleton

We can add the Ricker model deterministic skeleton to the parus pomp object. Since the Ricker model is a discrete-time model, its skeleton is a map that takes \(N_t\) to \(N_{t+1}\) according to the Ricker model equation \[N_{t+1} = r\,N_{t}\,\exp(-N_{t}).\] The following implements this.

skel <- Csnippet("DN = r*N*exp(-N);")

We then add this to the pomp object:

parus <- pomp(parus,skeleton=map(skel),paramnames=c("r"),statenames=c("N"))

Note that we have to inform pomp as to which of the variables we’ve referred to in skel is a state variable (statenames) and which is a parameter (paramnames). In writing a Csnippet for the deterministic skeleton, we use D to designate the map’s value. The map call tells pomp that the skeleton is a discrete-time dynamical system (a map) rather than a continuous-time system (a vectorfield).

With just the skeleton defined, we are in a position to compute the trajectories of the deterministic skeleton at any point in parameter space. For example, here we compute the trajectory and superimpose it on a plot of one simulation:

traj <- trajectory(parus,params=c(N.0=1,r=12),format="d")
plot(N~year,data=sim,type='o')
lines(N~year,data=traj,type='l',col='red')

A note on terminology

If we know the state, \(x(t_0)\), of the system at time \(t_0\), it makes sense to speak about the entire trajectory of the system for all \(t>t_0\). This is true whether we are thinking of the system as deterministic or stochastic. Of course, in the former case, the trajectory is uniquely determined by \(x(t_0)\), while in the stochastic case, only the probability distribution of \(x(t)\), \(t>t_0\) is determined. To avoid confusion, we use the term “trajectory” exclusively to refer to trajectories of a deterministic process. Thus, the trajectory command iterates or integrates the deterministic skeleton forward in time, returning the unique trajectory determined by the specified parameters. When we want to speak about sample paths of a stochastic process, we use the term simulation. Accordingly, the simulate command always returns individual sample paths from the POMP. In particular, we avoid “simulating a set of differential equations”, preferring instead to speak of “integrating” the equations, or “computing trajectories”.

Exercises

Exercise: Ricker model parameters

Fiddle with the parameters to try and make the simulations look more like the data. This will help you build some intuition for what the various parameters do.

Exercise: Reformulating the Ricker model

Reparameterize the Ricker model so that the scaling of \(N\) is explicit: \[N_{t+1} = r\,N_{t}\,\exp\left(-\frac{N_{t}}{K}+\varepsilon_t\right).\]

Modify the pomp object we created above to reflect this reparameterization.

Modify the measurement model so that \[\mathrm{pop}_t \sim \mathrm{Negbin}\left(\phi\,N_t,k\right),\] i.e., \(\mathrm{pop}_t\) is negative-binomially distributed with mean \(\phi\,N_t\) and clumping parameter \(k\). See ?NegBinomial for documentation on the negative binomial distribution and the R Extensions Manual section on distribution functions for information on how to access these in C.

Exercise: Beverton-Holt

Construct a pomp object for the Parus major data and the stochastic Beverton-Holt model \[N_{t+1} = \frac{a\,N_t}{1+b\,N_t}\,\varepsilon_t,\] where \(a\) and \(b\) are parameters and \[\varepsilon_t \sim \mathrm{Lognormal}\left(-\tfrac{1}{2}\sigma^2,\sigma\right).\] Assume the same measurement model as before.

Back to course homepage

R codes for this document

References

Bjørnstad, O. N., and B. T. Grenfell. 2001. Noisy clockwork: Time series analysis of population fluctuations in animals. Science 293:638–643.

He, D., E. L. Ionides, and A. A. King. 2010. Plug-and-play inference for disease dynamics: Measles in large and small populations as a case study. Journal of the Royal Society, Interface 7:271–283.

King, A. A., M. Domenech de Cellès, F. M. G. Magpantay, and P. Rohani. 2015. Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to ebola. Proceedings of the Royal Society of London. Series B 282:20150347.

King, A. A., E. L. Ionides, M. Pascual, and M. J. Bouma. 2008. Inapparent infections and cholera dynamics. Nature 454:877–880.

Laneri, K., A. Bhadra, E. L. Ionides, M. Bouma, R. C. Dhiman, R. S. Yadav, and M. Pascual. 2010. Forcing versus feedback: Epidemic malaria and monsoon rains in northwest India. PLoS Computational Biology 6:e1000898.

McCleery, R. H., and C. M. Perrins. 1991. Effects of predation on the numbers of great tits Parus major. Pages 129–147 in Bird population studies. Relevence to conservation and management. Oxford University Press, Oxford.

Romero-Severson, E. O., E. Volz, J. S. Koopman, T. Leitner, and E. L. Ionides. 2015. Dynamic variation in sexual contact rates in a cohort of hiv-negative gay men. American Journal of Epidemiology 182:255–262.